はじめに
連の検定についてよく分からなかったのでRで実装して確かめてみました。
連とは?
2値{A,B}を取る系列ABAABBAAABBBAB
があった時、これをA | B | AA | BB | AAA | BBB | A | B
と連続する同じ値ごとに分割できるように見える。この時同じ文字または記号のひと続きを連(run)
*1と言い、またひと続きの個数を連の長さ
*2と言う。この場合、Aの連が4つ(= {A,AA,AAA,A})とBの連が4つ(= {B,BB,BBB,B})あることが分かる。AAAAAAAAABBBBBBBBBBB
であれば、Aの連が1つでBの連が1つとなる。
確かめたこと
Aの数(= Na)とBの数(= Nb)が与えられた時、N = Na+Nbの大きさを持つ系列は、
を確かめた。
ソース
# 連の検定#### rm(list=ls()) # -------------------- # 関数: generateRun #### # -------------------- # 引数 # Na_: Aの個数 # Nb_: Bの個数 # 戻り値 # 連の数を返す関数 # -------------------- generateRun <- function(Na_, Nb_){ return( function(dummy){ N <- Na_ + Nb_ Arep <- rep(TRUE,Na_) Brep <- rep(FALSE,Nb_) AB <- c(Arep,Brep) xx <- sample(AB, N, replace = FALSE) yy <- xx[1:(N-1)] zz <- xx[2:N] return( sum( yy != zz) + 1 ) } ) } # シミュレーション #### ## 各種変数設定 Na <- 9 # Aの個数 Nb <- 11 # Bの個数 sim_N <- 10000 # シミュレーション回数 ## 平均値と分散 ( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値 ( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散 # シミュレーション実行 Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) ) # 要約 table(Result.Run) summary( Result.Run ) var( Result.Run ) # グラフ表示 hist(Result.Run,breaks = seq(min(Result.Run),max(Result.Run),1)) par(new=T) # 正規分布を重ね合わせ curve(dnorm(x,t_mean,sqrt(t_var)),min(Result.Run),max(Result.Run),xlab="",ylab="",xaxt="n", yaxt="n",col="red",bty="n")
結果の確認
数値の確認
> # シミュレーション #### > ## 各種変数設定 > Na <- 9 # Aの個数 > Nb <- 11 # Bの個数 > sim_N <- 10000 # シミュレーション回数 > ## 平均値と分散 > ( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値 [1] 10.9 > ( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散 [1] 4.637368 > # シミュレーション実行 > Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) ) > # 要約 > table(Result.Run) Result.Run 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1 3 13 28 167 312 852 1218 1725 1699 1699 1144 704 279 121 30 5 > summary( Result.Run ) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.00 9.00 11.00 10.89 12.00 18.00 > var( Result.Run ) [1] 4.666794
プロット
結果の考察
- 計算した平均値・分散がシミュレーションの結果とだいたい合っていることがわかった。
- ヒストグラムとカーブプロットが大まかにあっている。
おわりに
帰無仮説の下での平均値・分散の式が奇妙で覚えられないし式を見てもあまり実感が湧かなかったものの、シミュレーションしてみると確かに平均値・分散がそれっぽく合うのでこの式の妥当性を確認することができた。疑わしい系列があってもこれで確かめることができそうだ。
参考
連の検定関連
R言語のplot関連
- 49. 図の重ね合わせ
- 図の重ね合わせ方法を参考にしました。
- 正規分布のグラフをRで描く。【curve()の使い方】
- 正規分布をcurveで書く方法を参考にしました。
- R 使い方 軸・ラベルの調整(向き・サイズ・色など) グラフの描き方
- plotを重ねる時に軸やラベルを消す方法を参考にしました。