## Here, because n=6, we can actually enumerate all possible outcomes
## There are (n+c-1)! / [(c-1)!n!] different tables with c categories and given margin n
n = 6
luis = expand.grid(0:6,0:6,0:6) # all possible permutations of numbers 0:6
dim(luis)
6^3
luis1 = luis[rowSums(luis)==6,]  #filter out only those that have a sum of 6
dim(luis1)  #this should have (n+c-1)! / [(c-1)!n!] different rows, and it does!
prob0 = c(0.4,0.4,0.2) #null probabilities
table.probs = apply(luis1,1,function(ro) dmultinom(ro,size=6,prob=prob0)) #this gives the exact probability of every possible table under the null hypothesis
X2 = apply(luis1,1,function(ro) sum((ro-n*prob0)^2/(n*prob0)))   #compute X2 for eachpossible table
obs = c(1,3,2)
X2.obs = sum((obs-n*prob0)^2/(n*prob0))
p.exact = sum(table.probs[X2>=X2.obs])
p.exact
p.asympt = 1-pchisq(X2.obs,df=2)
p.asympt


## Generating the exact distribution by randomly simulating from the multinomial under H0:
n = 186
prob0 = c(0.4,0.4,0.2) #null probabilities
sims = 10000
tables = rmultinom(sims,n,prob0)
X2 = apply(tables,2,function(co) sum((co-n*prob0)^2/(n*prob0)))
obs = c(70,90,26)
X2.obs = sum((obs-n*prob0)^2/(n*prob0))
mean(X2>=X2.obs)
hist(X2)
abline(v=X2.obs, col="red")
chisq.test(obs, p=prob0, correct=FALSE) #this gives asymptotic P-value
chisq.test(obs, p=prob0, correct=FALSE, simulate.p.value=TRUE, B=10000) #this gives exact P-value


G2 = apply(tables,2,function(co) 2*sum(co*log(co/(n*prob0))))
G2.obs = 2*sum(obs*log(obs/(n*prob0)))
mean(G2>=G2.obs)
hist(G2)
abline(v=G2.obs, col="red")