r-statistics-fanの日記

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

ロケみつのサイコロの統計その1(イカサマサイコロなのか?)

はじめは、遊び記事から。

 

筆者はロケみつが好きである。

海外の学会で発表する機会があるのだが、学生時代から

英語は大の苦手である。仕事柄読むのはそれなりに出来る

のだが、英会話が全然ダメなのだ。センター試験でヒアリング

なんて無かったような気がする。外国人英語教師なんて大学

まで一度もお目にかからなかった。LとRなんで聞き分け不可能。

 

ロケみつ以前の海外学会では、

1:頭のなかでいうことを考える

2:英作文する

3:文法が正しいか吟味する

4:やっと発声する

というような感じであった。当然、会話のテンポがものすごく

悪く、会話にならない。3の時点で、あきらめることもしばしば。

 

ところが、ロケみつの早希ちゃん、ゆりやんレトリィバァさん

の雄姿をみて、何かが弾けました。

 

文法なんて細かいことが考えずに一生懸命話せば通じるやん!

 

ロケみつ後の海外では

1:とにかく思ったことをしゃべりまくる

2:通じなかったら、更にしゃべりまくる

 

ただ、これだけ。観光も仕事も捗りまくり。めっちゃ楽しかった。

 

と、いうことで、ロケみつに感謝はすれども、野暮な物言いを

する意図はありませんので。

 

http://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B1%E3%81%BF%E3%81%A4_%E3%82%B6%E3%83%BB%E3%83%AF%E3%83%BC%E3%83%AB%E3%83%89

ここから、出目を収得する

チャラD第弐百拾八話以降の出目は

サイコロ1~6が、それぞれ0,1,3,1,2,6回

 

これをRで解析する

deme <- c(0,1,3,1,2,6)  #サイコロの出目1~6

chisq.test(deme)

##出力

Chi-squared test for given probabilities

 

data:  deme

X-squared = 10.5385, df = 5, p-value = 0.06134

 

Warning message:

In chisq.test(deme) :  カイ自乗近似は不正確かもしれません 

 

これではいけないので青木先生のコードを使ってみる

http://aoki2.si.gunma-u.ac.jp/R/gft.html

 

source("http://aoki2.si.gunma-u.ac.jp/R/src/gft.R", encoding="euc-jp")

gft(deme, p=rep(1, 6))

 

カイ二乗値は 10.5385,自由度は 5,P値は 0.061339

正確なP値は 0.0683837

 

自作コードだと(我ながらメモリを贅沢に使って非効率的で遅いと思う)

deme <- c(0,1,3,1,2,6)  #サイコロの出目1~6

m <- sum(deme)  

M <- choose(6 + m - 1, 5) 

x <- matrix(0, nrow=M, ncol=6)    

n <- 1  

y <- numeric(M)

p <- dmultinom(deme, prob = c(rep(1/6,6))) 

 

for (n1 in 0:m){

  for (n2 in 0:m){

    if (n1 + n2 > m) {next}

    else {

      for (n3 in 0:m){

        if (n1 + n2 + n3 > m) {next} 

        else {  

          for (n4 in 0:m){

            if (n1 + n2 + n3 + n4 > m) {next} 

            else  {

              for (n5 in 0:m){

                if (n1 + n2 + n3 + n4 + n5 > m) {next} 

                else{

                  x[n,1] <- n1

                  x[n,2] <- n2

                  x[n,3] <- n3  

                  x[n,4] <- n4

                  x[n,5] <- n5

                  x[n,6] <- m - n1 - n2 - n3 - n4 - n5

                  n <- n + 1

                  next}

              }

            }

          }

        }

      }}}}

 

y <- apply(x, 1, function(x) dmultinom(x, prob = c(rep(1/6,6)))) #それぞれの組み合わせの確率

sum(y[y <= p])  #正確なP値 #計算誤差を考慮しても結果は不変のようだ

 

 

###出力

> sum(y[y <= p])  #正確なP値

[1] 0.08783549

 

あれ、なんか違う

 

カイ二乗:P=0.06134

aoki先生正確検定:P=0.0683837

自作:P=0.08783549

 

いずれにせよ、P>0.05なので、サイコロは怪しくない範囲とい

分析結果になる。

なんで結果に差がでたかというと、aoki先生のコードはカイ二乗値で、

より極端な表を同定しているのに対して、自作コードは多項分布の

生成確率で極端な表を同定しているためである。

 

たとえば、5,4,2,2,0,0という出目は、カイ二乗基準だと極端ではない

という判断になるが、多項分布だと極端という判定になる。

カイ二乗5,4,2,2,0,0: 9.615385 #出目よりありふれていると判断

カイ二乗0,1,3,1,2,6:10.53846 

多項分布確率5,4,2,2,0,0:4.138677e-05 #出目より極端と判断

多項分布確率0,1,3,1,2,6:5.518237e-05 

このように異なる判断が生じる組み合わせがあるからである。

 

いずれにせよ、サイコロの出目は問題なしと結論。よかった。

 

多項分布での正確確率の記事は検索した限り無い?ので、面白い

結果が出て自己満足。色々値を変えてみた範囲で、自作コードの

方がP値は大きくなる傾向があるようだ(全く同じ値のこともある)。

有意差を出したい向きには自作コードは需要はないな。