RでLepage Test

血尿でるまで研究頑張ろうとおもったけど先にア●ルから血がでてきました.

あっ, どーも僕です. みなさんは元気ですか.

Lepage Test

レポートには課されていないのですが, マンケンドール検定とともにせっかく教わったので(・・・というか研究がうまくいかない現実逃避です.)ラページ検定をRで組んでみました.

ラページ検定ってなんだというと, ウィルコンスク検定とアンサリー・ブラッドレイ検定を同時におこなってみれば自由度２のカイ二乗分布に従うので, ノンパラメトリックに不連続的変化の検出できるよねって話です.

詳しくは統計科学研究所のサイトにわかりやすい資料があります.
学習参考資料の統計データ解析入門にノンパラメトリック検定のPDFがあります.
統計科学研究所

で, これと授業ノートにしたがってやってみてそれっぽい図がかけたのですが先生がくれた資料と結果が違う・・・. もう生きててごめんなさいって感じです.

Rのスクリプトはこちら. データは先生が見本として見せくれた金沢市の月平均気温です.

#****************************************************
# 金沢市における
# 月平均気温を
# ラページ検定する
#
# データ源 気象庁HP　2012/11/16
# http://www.data.jma.go.jp/obd/stats/etrn/index.php?prec_no=57&block_no=47616
#***************************************************
# 動くけど先生と同じにならない・・・・・詰んだわ

# *********************** temprature ****************************
# ラページ検定の検定統計量
compute_lepageHK <- function ( d, k, m)
{
  # d  : 一次元データ
  # k: 境界よりも前のデータ数
  # m: 境界よりも後のデータ数
  N <- length(d) 
  border_index <- which((1:N > k) & (1:N < (N - m + 1)))# 境界となりうるデータ位置
  n   <- k + m
  UI  <- sapply(border_index, compute_ui, k, m, d)
  
  # HK算出のための各項を計算
  
  # ウィルコクソン検定統計量
  W   <- apply( UI, 2, function(x) sum(1:n * x))
  EW  <- k * (n + 1) / 2
  VW  <- k * m * (n + 1) / 12
  
  # アンサリー・ブラッドレー検定統計量
  A   <- apply( UI, 2, compute_A, k, m)
  if (n %% 2 == 0)
  {
    EA <- k * (n + 2) / 4
    VA <- k * m * (n - 2) * (n + 2) / (48 * (n - 1) )[f:id:aaaazzzz036:20180120103143p:plain]
  } else
  {
    EA <- k * (n + 1)^2 / (4 * n)
    VA <- k * m * (n + 1) * (n^2 + 3) / (48 * n^2)
  }
  
  lepageHK <- (W - EW)^2 / VW + (A - EA)^2 / VA
  return(lepageHK)
}


# uを計算する関数
compute_ui <- function ( i, k, m, d)
{
  # i : 境界となるデータ位置
  x <- d[i - (k:1)]
  names(x) <- rep(1, k)
  y <- d[i + 1:m]
  names(y) <- rep(0, m)
   
  rtn <- as.numeric(names(sort(c(x,y))))
  return(rtn)
}

# Aを計算する関数
compute_A <- function ( d, k, m)
{
  tmp <- (k+1):(k+m)
  a1 <- sum(1:k * d[1:k])
  a2 <- sum((k + m - tmp + 1) * d[tmp])
  return(a1 + a2)
}

# Aを計算する関数
# (結果が同じだから使わないけどこっちが定義どおり)
compute_AA <- function ( x, k, m, n)
{
  a1 <- k * (n + 1) / 2
  a2 <- sum(abs(1:n - (n + 1)/2) * x)
  return(a1 - a2)
}


# データ
dat <- read.table("temperature_kanazawa.txt", 
                  header = T, sep = "\t",
                  row.names = 1)

# 結果
result <- sapply(dat, compute_lepageHK, 20, 20)

#  2月だけ書いてみる
def.par <- par(no.readonly = TRUE) # デフォルトの環境を保存
par(mai = rep(1.5, 4))
x.axis <- as.numeric(rownames(dat))
plot(x.axis, dat[,2], type = "n", 
     ylab = "Mean Temprature", xlab = "Year", 
     tcl = 0.5, ylim = c(0,1), yaxt = "n",
     main = "February")
axis(side = 2, at = seq(0, 1, length = 5), tcl = 0.5, 
     labels = round(seq(0, max(dat[,2]), length = 5), 1))
axis(side = 4, at = seq(0, 1, length = 5), tcl = 0.5,
     labels = round(seq(0, max(result[,2]), length = 5), 1))
points(x.axis, dat[,2]/max(dat[,2]), col = "#0068b7", 
       lwd = 3, type = "l")
points(x.axis[-c(1:20, (nrow(dat)-19):nrow(dat))], 
       result[,2]/max(result[,2]), lwd = 3,
       type = "l", col = "#f39800")
mtext(side = 4, text = "LepageHK", line = 3)
abline( h = qchisq(c(0.95, 0.99), df = 2)/max(result[,2]), lty = 5)
text(locator(1), labels = "1%", cex = 2)
text(locator(1), labels = "5%", cex = 2)
legend("bottomright",  bty = "n", legend = c("Temp", "HK"), 
       lty = 1, col = c("#0068b7", "#f39800"), lwd = 1, cex = 1)
par(def.par)　　　　　　　　　　　 # デフォルトの環境を復帰

2月だけグラフを書いてみた結果がこちら. HKが破線を超えていると不連続が有意ですねって話です. 見た目は良い感じに評価できてそうだけど, どこか間違っている. けどなにが間違っているかわからぬ. ぐぬぬ.

教えてエロい人！！

追記

おそらく, i番目を境界に20個ずつとするときにi-20~i-1の二十個とi+1~i+20の二十個のデータではなくi-20~i-1の二十個とi~i+19の二十個のデータでやれば良さげ

再実装（2018/1/20）

コメントで計算がおかしいとの指摘がありました。指摘いただき、ありがとうございます。そこで改めて実装してみることにしました。検証方法はざっくりとですが、論文の追試を行いました。

対象の論文は「Discontinuous Changes of Precipitation in Japan after 1900 Detected by the Lepage Test, Tsuneharu Yonetani, Journal of the Meteorological Society of Japan, Vol.70, No.1」です。
Googleスカラーに伺いを立てるとPDFを得ることができます。
この論文では気象庁のデータを使いラページ検定で年降水量の変化点検出をしています。使われているデータは気象庁HPから取得することができるので、誰でも追試が可能です。

ラページ検定の結果はFig2.です。この図を再現することとしました。
結果がこちらです。上がラページ検定の結果で、下がデータのグラフです。ここには載せられませんが論文の図と見比べてもらうと、再現できていることがわかります。

なお、論文と比べるといくつかデータが異なっているようですが、数点異なるくらいでは結果に大きな影響はないと思います。

再実装したプログラムは下記となります。lepageHKが統計量を計算する関数です。

library(tidyverse)
convert_to_ui <- function(x, y) {
    names(x)  <- rep("1", length(x))
    names(y)  <- rep("0", length(y))
    sorted_xy <- sort(c(x,y))
    parse_integer(names(sorted_xy))
}
compute_W <- function(ui) {
    sum( seq_along(ui) * ui )
}
compute_EW <- function(n1, n2) {
    .5 * n1 * (n1 + n2 + 1)
}
compute_VW <- function(n1, n2) {
    n1 * n2 * (n1 + n2 + 1) / 12
}
compute_A <- function(ui, N) {
    n  <- N / 2
    a1 <- sum(1:n * ui[1:n])
    a2 <- sum((N-(n+1):N + 1) * ui[(n+1):N])
    a1 + a2
}
compute_EA <- function(n1, n2) {
    N <- n1 + n2
    if(N %% 2 == 0)
        n1 * (N + 2) / 4
    else
        n1 * (N + 1) ** 2 / 4 / N
}
compute_VA <- function(n1, n2) {
    N <- n1 + n2
    if(N %% 2 == 0)
        n1 * n2 * (N - 2) * (N + 2) / 48 / (N - 1)
    else
        n1 * n2 * (N + 1)(N ** 2 + 3) / 48 /(N **2)
}
lepageHK <- function(x, y) {
    n1 <- length(x)
    n2 <- length(y)
    N  <- n1 + n2

    ui <- convert_to_ui(x, y)

    W  <- compute_W(ui)
    EW <- compute_EW(n1, n2)
    VW <- compute_VW(n1, n2)

    A  <- compute_A(ui, N)
    EA <- compute_EA(n1, n2)
    VA <- compute_VA(n1, n2)

    (W - EW) ** 2 / VW + (A - EA) ** 2 / VA
}

今回のように時系列データへの適用は次のようなラッパーで対応が可能です。
dat1dにベクトルでデータを与えて、左側のサイズと右側のサイズを与えてやります。

lepage_test_values <- function(dat1d, left_size, right_size) {
    N <- length(dat1d)
    test_values <- rep(NA_real_, N)
    i   <- left_size + 1
    lid <- 1:left_size
    rid <- left_size + 1:right_size
    while(max(rid) <= N) {
        test_values[i] <- lepageHK(dat1d[lid], dat1d[rid])
        i   <- i + 1
        lid <- lid + 1
        rid <- rid + 1
    }
    test_values
}