読者です 読者をやめる 読者になる 読者になる

RでLepage Test

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


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


Lepage Test

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

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

詳しくは統計科学研究所のサイトにわかりやすい資料があります.
学習参考資料の統計データ解析入門にノンパラメトリック検定の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) )
  } 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が破線を超えていると不連続が有意ですねって話です. 見た目は良い感じに評価できてそうだけど, どこか間違っている. けどなにが間違っているかわからぬ. ぐぬぬ.

教えてエロい人!!


f:id:aaaazzzz036:20121128001740j:plain

追記

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