#Generate a table of N response with I response levels at each of T timepoints
tablesim=function(N,I,T){
ncells=I^T
p=matrix(0,T,I)
indx=NULL
for(i in 1:T){
#Generate the table indices
indx=cbind(indx,(rep(rep(1:I,each=I^(i-1)),I^(T-i))))
#Generate marginal probabilities of approx 1/I
px=rbeta(I,1,I-1)
px=px/sum(px)
p[i,]=px
}
#Generate a mutual independence model (use of outer product here could
#possibly take too much space
prob=p[1,]
for(i in 2:T){
prob=prob%o%p[i,]
}
prob=as.numeric(prob)
#Add a little noise to the mutual independence model
depadj=rbeta(ncells,1,1)
depadj=(depadj-mean(depadj))/10
prob=prob+depadj
if(min(prob)<0) prob=prob-min(prob)
prob=prob/sum(prob)
#Simulate N random draws
x=rmultinom(1,N,prob)
tab=cbind(indx,x)
return(tab)
}
#Save some data sets
tab=tablesim(1000,3,3)
write(t(tab),"z:/stat 770/3tothe3.txt",ncol=4)
tab=tablesim(1000,3,4)
write(t(tab),"z:/stat 770/3tothe4.txt",ncol=5)
tab=tablesim(10000,3,12)
write(t(tab),"z:/stat 770/3tothe12.txt",ncol=13)