8月29日（金）午後

chisq.test()の自作（3種類で9回試行する場合）

「関数をつくる」ことと「関数で計算する」ことをごちゃまぜにしていたので、ちゃんと分けて考えることにする。

あと、「何の関数なのか」のメモも基本的にするようにして、整理しながら進めようと思う。

f1<-function(N,i){
y<-c()
for (j in 0:(N-i)){
x<-c(i,j,N-i-j)
y<-c(y,x)}
;y
}# iを固定したベクトル

f2<-function(N){
z<-c()
for(i in 0:N)
z<-c(z,f1(N,i));z
}# f1を使って全部の出方を一列にしたベクトル

mat1<-matrix(f2(9),ncol=3,byrow=TRUE)
#　f2を行列にする

f3<-function(N,x,y,z){
log((factorial(N)*(1/3)^N)/
(factorial(x)*factorial(y)*factorial(z)))
}# 試行回数がNのとき、それぞれの出方がx,y,zとなるような確率

f4<-function(n,N){
factorial(N+n-1)/(factorial(N)*factorial(n-1))
}# 試行回数がN、n種類の出方があるときの出方の総数

f5<-function(N){
x<-c()
for(i in 1:f4(3,N))
x<-c(x,f3(N,mat1[i,1],mat1[i,2],mat1[i,3]))
;x
}# mat1からそれぞれの出方になる確率のベクトルをつくる

f6<-function(N,x,y,z){
w<-f5(N)[f5(N)<f3(N,x,y,z)]
return(sum(exp(w)))
}#f3の値よりも小さいf5の要素の和＝p値

mychisq.test<-function(N,c){
x<-c[1]
y<-c[2]
z<-c[3]
f6(N,x,y,z)
}#s.tの形で入れられるように少し調整

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

mychisq.test(9,s.t)

> mychisq.test(9,s.t)
[1] 0.5305594

それらしいものができた。

chisq.test()と同じといえるか確かめたいところではあるが、ちょっといい方法が思い浮かばない。

出てくる数字を眺めた限りでは、出てくる値は違うが、バラつき方は似たようなものに思える。

mychisq.test（）の拡張ができてからもう少し考えることにする。

以上