適当に分割表を作る

• 総サンプル数を固定して、行数、列数を指定して適当に周辺度数を与えて、その期待度数の表を作る

# <desc> 総サンプル数がnで、行数がs、列数がtの分割表を適当に作る。自然数nをsに分けるとき、長さnの数列のn-1の分割可能箇所からs箇所を選ぶことで可能なので、そのような選択をsample()を用いた定義関数Partition()で作成してやる </desc>


IntegerTable<-function(n=10,s=4,t=3){
 S<-Partition(n,s)$partition
 T<-Partition(n,t)$partition
 ret<-matrix(rep(1,s*t),nrow=s)
 for(i in 1:s){
  for(j in 1:t){
   ret[i,j]<-S[i]*T[j]/n
  }
 }
 write.table(ret,file="TXTFILE1.txt",sep="\t",quote=F,row.names=F,col.names=F)
 return(ret)
}

Partition<-function(n=10,k=3){
 partpnt<-sample(1:(n-1),k,replace=FALSE)
 partpnt<-c(sort(partpnt),n)
 partpnt2<-c(0,partpnt[1:(length(partpnt)-1)])
 return(list(partition=partpnt-partpnt2))
}

• 期待度数がすべて自然数の表を作る。

# <desc> 行数がs、列数がtので期待度数がすべて自然数であるような分割表の期待度数表を適当に作る。行と列の周辺度数を自然数のベクトルで与え、セルの値をその積とする</desc>

MakeIntegerExpTable<-function(S=rpois(3,3)+1,T=rpois(4,3)+1){
 ret<-matrix(rep(0,length(S)*length(T)),nrow=length(S))
 mrow<-rep(0,length(S))
 mcol<-rep(0,length(T))
 for(i in 1:length(S)){
  for(j in 1:length(T)){
   ret[i,j]<-S[i]*T[j]
   mrow[i]<-mrow[i]+ret[i,j]
   mcol[j]<-mcol[j]+ret[i,j]
  }
 }
 
 return(list(ratioRow=S,ratioCol=T,table=ret,mrow=mrow,mcol=mcol,N=sum(ret)))
}

• "r2dtable()" 2d(二次元)のtableをr(ランダム発生)というコマンドがあると教えていただきました

n<-1000
r<-c(10,20,30)
c<-c(5,40,10,5)

ts<-r2dtable(n,r,c)
ps<-rep(0,n)
for(i in 1:n){
 ps[i]<-chisq.test(ts[[i]])$p.value
}

plot(sort(ps))