# Bret Larget
# Simulation for p-value for example
# Compare to chisq-test statistic p-value

# Observed data
o = matrix(c(3,4,47,46),2,2)

# Compute expected counts
x = apply(o,1,sum)
y = apply(o,2,sum)
n = sum(o)
e = (x %o% y) / n

# Probabilities for likelihood ratio test
theta.null = e/n
theta.full = o/n

# Function to find the test-statistic for the chi-square test
test.stat = function(m) {
  n = sum(m)
  x = apply(m,1,sum)
  y = apply(m,2,sum)
  e = (x %o% y) / n
  if(!any(e==0))
    ts = sum((m-e)^2/e)
  else
    ts = 0
  return(ts)
}

# Function to return a sample drawn from the null MLE
sample.m = function(m) {
  s = sample(4,sum(m),replace=T,prob=(apply(m,1,sum) %o% apply(m,2,sum))/sum(m)^2)
  counts = c(sum(s==1),sum(s==2),sum(s==3),sum(s==4))
  return(matrix(counts,2,2))
}

# Compute the test-statistic for 10,000 simulated data sets.
# Is this true distribution similar to chi-square(1)?
out = rep(NA,10000)
for(i in 1:10000)
  out[i] = test.stat(sample.m(o))
p.sim = sum(out >= test.stat(o))/10000
p.chisq = 1 - pchisq(test.stat(o),1)
cat(paste("Simulation-based p-value =",round(p.sim,4),"\n"))
cat(paste("Chi-square test  p-value =",round(p.chisq,4),"\n"))