9月1日（月）午前
気を引き締めてかかりたいところである。
mychisq.test()の自作
とりあえず求めるものは出来上がった。
mat1<-function(n,N){
if (N==0) return (matrix(0,ncol=n))
else if(n==1)return(matrix(N))
else{
tmp.list<-list()
cnt<-1
for (i in 0:N){
tmp.out<-cbind(mat1(n-1,N-i),i)
tmp.list[[cnt]]<-tmp.out
cnt<-cnt+1
}
mat2<-numeric(0)
for(j in 1:(N+1)){
mat2<-rbind(mat2,tmp.list[[j]])
}
return(mat2)
}
}#n種でN回試行するときの出方の行列
f1<-function(c){
x<-0
for(i in 1:length(c))x<-x+lfactorial(c[i])
return((lfactorial(sum(c))+sum(c)*log(1/length(c)))-x)
}# それぞれの出方がベクトルcとなるような確率(の対数）
f2<-function(n,N){
x<-c()
for(i in 1:nrow(mat1(n,N)))x<-c(x,f1(mat1(n,N)[i,]))
return(x)
}#mat1からそれぞれの出方になる確率（の対数）のベクトルをつくる
mychisq.test<-function(c){
x<-f2(length(c),sum(c))[f2(length(c),sum(c))<f1(c)]
return(sum(exp(x)))
}#p値を計算する関数
s <- sample(1:3,9,replace=TRUE)
s.t <- table(s)
mychisq.test(s.t)
> mychisq.test(s.t)
[1] 0.5305594
s <- sample(1:5,10,replace=TRUE)
s.t <- table(s)
mychisq.test(s.t)
> mychisq.test(s.t)
[1] 0.5626086
というわけで、任意の種類・試行回数について適応度検定を行える関数ができた。
しかし、後者（5種類・10回）の検定の時点ですでに1分弱くらいかかっている。かなり負担がかかる関数なのだと思われる。
mychisq.test()の検証
ここで、一応もともとRに備えられているchisq.test()と結果がだいたい同じといえるかどうか、検定をしてみたい。
結果の母平均と母分散がだいたい同じなら同じといってよいのではないか、ということで、やってみる。
s <- function(){
sample(1:3,20,replace=TRUE)
}
s.t <- function(){
table(s())}#テーブル作成
x<-c()
for(i in 1:100)x<-c(x,mychisq.test(s.t()))
y<-c()
for(i in 1:100)y<-c(y,chisq.test(s.t())$p.value)
#データをベクトルにする
t.test(x,y,paired=FALSE)
> t.test(x,y,paired=FALSE)
Welch Two Sample t-test
data: x and y
t = -0.2052, df = 197.259, p-value = 0.8376
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.08615321 0.06991538
sample estimates:
mean of x mean of y
0.4586159 0.4667348 最初はN=50でやってみたのだが、計算に時間がかかっていつまでも終わらなかったため、N=20としたところ、昼休みの間に終わっていた。
平均に関しては同じと言ってよさそうである。
午後は分散についてやってみる。（時間がどれだけかかるか分からないが…。）
以上
