r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

ロケミツのサイコロの統計その2(連続して1が出ないのが20回続いたらイカサマ?)

ロケミツのサイコロの統計(連続して1が出ないのが20回続いたら)

#どうもグラフがうまくスマホで表示できない。グラフが見えなかったらごめんなさい。

 

1が連続20回出ていないらしい。 この確率を求める。

単純に1が連続20回でないのは、(5/6)20 =0.026 で2.6%しか 起こらない珍しい事象

しかし、実際はこれまでに、 http://gavangavan.hatenablog.com/entry/20131130/1385780461 を参考にすると、20+35+10+106+24+11+35=最低でも241回はサイコロ を振っているので、もっと起こりやすいと予想される。

これをRで計算してみよう。

m <- 241  #サイコロを振る回数
q <- 5/6  #たとえば'1が出ない'確率
r <- 20  #qの確率で起こる事象が連続r回起こる

prenzoku <- function(m, q, r) {
    pren <- numeric(m - r + 1)

    pren[1] <- q^r
    pren[2:(r + 1)] <- (1 - q) * q^r

    for (i in (r + 2):(m - r + 1)) {
        pren[i] <- (1 - q) * q^r * (1 - sum(pren[1:(i - r - 1)]))
    }

    return(sum(pren))  #答えの確率
}

prenzoku(m, q, r)  #1回以上”qの確率で起こる事象が連続r回起こる”確率
## [1] 0.6618

しかし、統計素人なので、間違っているかもしれないので検証する。

先人たちの業績を参考に検証する http://questionbox.jp.msn.com/qa4572311.html

prenzoku(5,½,2)=19/32 prenzoku(1000,½,15)=p[1000] / 21000 = 0.01495106644108 ならばOK

options(digits = 15)

prenzoku(1000, 1/2, 15)
## [1] 0.0149510664410799
0.01495106644108
## [1] 0.01495106644108

prenzoku(5, 1/2, 2)
## [1] 0.59375
19/32
## [1] 0.59375

うむ合っている

prenzoku(m, q, r)
## [1] 0.661756740764489

ということで、20回連続1が出ないのは66%くらいで 起こりうるありふれた事象ということになる。

今後何回連続1がでなかったら5%以下でしか無い珍しい事象になるか計算する

kk <- 20:100
kp <- numeric(length(kk) + 19)
for (k in kk) {
    kp[k] <- prenzoku(m + k - 20, q, k)
}
plot(kp[kk])

 

kp[36]
## [1] 0.0523819074409424
kp[37]
## [1] 0.0437959345660729

というわけで、37回連続して1が出なければ、やっと5%以下でしかおこらない

珍しい事象ということになる。(すでに起こっていることなどは今回考慮しない。

あくまでその時点で連続n回1が出ないことが最低1回は起こるだろう確率。)

すでに20回連続してでているので、あと17回か。先は長いな。

 

ということで、今回もロケみつのサイコロのイカサマ判定は白!