8月28日（木）午後 - 2014年度マイコースプログラム

山田先生からの課題の再考

原因を確かめようと調べてみたところ、

x<-c()
for(i in 1:10000)x<-c(x,myfunc01())
x.t<-table(x)
chisq.test(x.t)

> chisq.test(x.t)

        Chi-squared test for given probabilities

data:  x.t
X-squared = 1500.347, df = 3, p-value < 2.2e-16

y<-c()
for(i in 1:10000)y<-c(y,myfunc02())
y.t<-table(y)
chisq.test(y.t)

> chisq.test(y.t)

        Chi-squared test for given probabilities

data:  y.t
X-squared = 710.7872, df = 2, p-value < 2.2e-16

どうやらmyfunc01,myfunc02の時点でおかしかったようである。

よく見ると、

myfunc01<-function(){
if (runif(1)>3/4)return(4)
else if (runif(1)>2/4)return(3)
else if (runif(1)>1/4)return(2)
else return(1)
 }

ここで、runif(1)のたびに乱数を生成していたようである。それだと1～4の出る確率はそれぞれ、3/32、9/32、3/8、1/4となってしまう。

とりあえずここを修正。

myfunc01<-function(){
A<-runif(1)
if (A>3/4)return(4)
else if (A>2/4)return(3)
else if (A>1/4)return(2)
else return(1)
 }#1～4を同様に確からしく抽出する関数
 
myfunc02<-function(){
A<-runif(1)
if (A>2/3)return(3)
else if (A>1/3)return(2)
else return(1)
 }#1~3を同様に確からしく抽出する関数

また、

myfunc04<-function(){
z<-c()
for(i in 1:2)z<-c(z,10^(2-i)*myfunc03()[i])
return(sum(z))
}

ここも、for-loop処理のたびに別の乱数を生成してしまうので、問題あり。

正しくは

w<-function(){
myfunc03()
}
myfunc04<-function(x){
v<-x()
s<-c()
for(i in 1:2)s<-c(s,10^(2-i)*v[i])
return(sum(s))
}

u<-c()
for(i in 1:10000)u<-c(u,myfunc04(w))

こんなかんじにすべきだったようである。乱数をうまい形でループ処理するのに非常に苦労した。

まとめると、

myfunc01<-function(){
A<-runif(1)
if (A>3/4)return(4)
else if (A>2/4)return(3)
else if (A>1/4)return(2)
else return(1)
 }#1～4を同様に確からしく抽出する関数
 
myfunc02<-function(){
A<-runif(1)
if (A>2/3)return(3)
else if (A>1/3)return(2)
else return(1)
 }#1~3を同様に確からしく抽出する関数

myfunc03<-function(){
x<-c(1,2,3,4)
y<-x[c(-myfunc01())]
z<-y[c(-myfunc02())]
return(z)
}

myfunc04<-function(x){
v<-x()
s<-c()
for(i in 1:2)s<-c(s,10^(2-i)*v[i])
return(sum(s))
}

u<-c()
for(i in 1:10000)u<-c(u,myfunc04(w))

u.t<-table(u)

chisq.test(u.t)

> chisq.test(u.t)

        Chi-squared test for given probabilities

data:  u.t
X-squared = 2.6396, df = 5, p-value = 0.7553

というわけで、今度こそ一様に取り出す関数が出来た、と言える。

chisq.test()を自作する

chisq.test()を使って検定したはいいものの、そもそもchisq.test()って一体なんなのかよくわかっていなかった。

というわけで、chisq.test()の自作を試みることにした。

上の場合は、12,13,14,23,24,34の6種類で10000回試行した場合のchisq.test()を行ったわけであるが、まずは3種類で9回行う場合から考えてみる。

すなわち、

s <- sample(1:3,9,replace=TRUE)
s.t <- table(s)

このp値を求めてみる。

とりあえず各事象が起こる確率を求めたいというわけで、

> for (i in 0:9) 
+ for (j in 0:(9-i)) 
+ {if (i<0)next; if (j<0)next;cat(i, j, 9-i-j,"\n")}
0 0 9 
0 1 8 
0 2 7 
0 3 6 
0 4 5 
0 5 4 
0 6 3 
0 7 2 
0 8 1 
0 9 0 
1 0 8 
1 1 7 
1 2 6 
1 3 5 
1 4 4 
1 5 3 
1 6 2 
1 7 1 
1 8 0 
2 0 7 
2 1 6 
…（以下略）

これをなんとか行列の形で表したいのだが、いきなりここで行き詰まる。

まだ行列を使ったことがないので、習熟が必要なようだ。

明日フレッシュな頭で考えることにして、今日は帰って他の勉強でもすることにする。

以上