Lepage検定は時系列データなどを扱う時に使うことがあるのだけど、そのままノンパラの検定として使うことは少なくて、統計量を見て不連続的変化の検出に使う用途が多い（と思っている）。その用途に使う場合、とにかく繰り返し検定を行うことになるのでいちいち手書きで全組み合わせを書くのも面倒で、量が増えてくれば現実的でなくなってくるため、まとめてやってくれるコードを書いた。
まずはLepage検定を使えるようにする。
https://rdrr.io/github/tpepler/nonpar/
のInstallationにある
install.packages("remotes")
remotes::install_github("tpepler/nonpar")
をRで実行しましょう。
ネットワークの都合などでうまくできない場合は以下のページから関数のコードをコピーして貼り付ける。
https://github.com/tpepler/nonpar/blob/master/R/lepage.test.R
パッケージ NSM3 が必要。パッケージをインストールするときに関連するパッケージも全部インストールするのチェックボックスが出る場合はそれにチェックを入れないと正しくインストールされない。install.packages関数でインストールする場合はdependencies=Tとする。Rでは基本的なことらしいがここで一回つまづいたので念のため。
以下の.txtファイルにBFlepage.testの関数でラページ検定をデータに対して総当たりで行えるようにしたコードがある。BFlepage.txtの中身をRにコピペすればBFlepageという関数ができる。
BFlepage.test(x, n = NA, method = T/F)で、xはデータ（6個以上）、nは前後nずつで区切る数（3以上）、method = Tの時はデータを前後nずつに区切って検定を行った統計量のみを出す。method = F、あるいは指定なしの場合は端から3つとn、4つとnというふうに末端に近い部分を丁寧に見ていく。xの数に対して大きなnを指定した時もエラーは返さずmethod = Fとして処理を行う。統計量はhklistに記録される。BFlepage.testを実行後に
hklist
と入力すると計算された統計量の一覧が出る。
時系列データでLepage検定を用いた研究を見ているとnを適当に決めているが、最適なnというのはデータの構造次第なところがある。nが小さければ小さな不連続点も発見しやすいけれど検出力が低くノイジーで、nが大きいとノイズは減るが細かな変化が見えなくなることがある。そのため多くの場合はさまざまなnについて見てやることが重要である。
そこで、nを3からとりうる最大の数まで網羅的にBFlepage.testを実行して結果のhklistをデータフレームに保存、それをさらにグラフにして視覚化するには以下。このあたりから計算量が指数関数的に増加するので最初は少なめのデータから試して計算にかかる時間を確認すること。###(xに入れるデータの個数-3)をここにいれる###と###ここにデータを入れる###に数字とデータを入れる。
noxmat <- NULL
xmat <- NULL
rownamesxmat <- NULL
colnamesxmat <- NULL
for (i in 3:(###(xに入れるデータの個数-3)をここにいれる###)){
BFlepage.test(x = ###ここにデータを入れる###, n = i)
assign(sprintf("hklist%d", i), hklist)
xmat <- cbind(xmat, get(sprintf("hklist%d", i)))
noxmat <- append(noxmat, i)
rownamesxmat <- append(rownamesxmat, i)
colnamesxmat <- append(colnamesxmat, sprintf("hklist%d", i))
}
rownames(xmat) <- rownamesxmat # 行の名前
colnames(xmat) <- colnamesxmat # 列の名前
matplot(noxmat, xmat, type = 'l', xlab="x", ylab="HK")
最後のmatplotでグラフを書いてLepage検定の統計量を視覚的に示してくれる。
データフレームは
xmat
と入れると出てくるので全体像を確認できる。
たとえばx=61の時に最高の統計量HKを出しているnを知りたいときは
which.max(xmat["61",])
とするとhklistXXのように出力してくれる。このXXがnの値。""で囲まないと61行目のデータについて返すのでずれが発生するので注意。
統計量の最大値は
max(xmat["61",])
で求められる。
ある程度プログラミングができる人が上記のコードを見れば、処理が遅くなるのでよくないとされているfor文が多用されていることに気づくでしょう。計算速度を上げるならfor文を速い処理に書き直して、並列化をして利用するCPUのコア数を増やしてやれば数百〜数万倍の速度を得られるはずです。改良する気（と能力）は少なくとも今回のコードを書いた段階ではありませんので、時間がかかって困る方は自力で速くなるように組み直してください。