axjack's blog

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.

連の検定(Runs Test for Detecting Non-randomness)

はじめに

連の検定についてよく分からなかったので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の大きさを持つ系列は、

  • もしA・Bがランダムに現れるなら(←帰無仮説)連の数はどのような分布になるのか → 正規分布
  • 連の数の平均値・分散が以下の式*3どおりとなるか

f:id:axjack:20190926230103p:plain

を確かめた。

ソース

# 連の検定####
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

プロット

f:id:axjack:20190926230756p:plain

結果の考察

  • 計算した平均値・分散がシミュレーションの結果とだいたい合っていることがわかった。
  • ヒストグラムとカーブプロットが大まかにあっている。

おわりに

帰無仮説の下での平均値・分散の式が奇妙で覚えられないし式を見てもあまり実感が湧かなかったものの、シミュレーションしてみると確かに平均値・分散がそれっぽく合うのでこの式の妥当性を確認することができた。疑わしい系列があってもこれで確かめることができそうだ。

参考

連の検定関連

R言語のplot関連

*1:[連の検定]より引用

*2:[連の検定]より引用

*3:式は[連の検定]より引用