simtab = function(P, n)
{
  ## Convert to cumulative probabilities.
  CP = cumsum(P)
  ## Storage for the data being simulated.
  N = array(0, c(2,2))
  ## Simulate one contingency table.
  U = runif(n)
  N[1,1] = sum(U <= CP[1])
  N[1,2] = sum( (U > CP[1]) & (U <= CP[2]) )
  N[2,1] = sum( (U > CP[2]) & (U <= CP[3]) )
  N[2,2] = sum(U > CP[3])
  return(N)
}

n=1e4
P=(1:4)/sum(1:4)
N=simtab(P, n)
F=N/n
round(F, 3)

N=t(array(rmultinom(n=1,size=n,prob=P),c(2,2)))
F=N/n
round(F, 3)




reps=1000
##  Population Cell probs
P=(1:4)/sum(1:4)

POWER=NULL
for ( n in seq(from=10, to=5000, by=100)  )
{
  ## Generate reps tables with sample size n
  N=rmultinom(n=reps,size=n,prob=P)

  ## Compute the sample log-odds ratios, and SEs
  SLOR = log(N[1,]*N[4,]/N[2,]/N[3,])
  SE = sqrt(1/N[1,] + 1/N[2,] + 1/N[3,] + 1/N[4,])

  ## Compute the 95% CI
  LB = SLOR-2*SE
  UB = SLOR+2*SE

  ## Get a list of indicators of
  ## CI's not including zero
  REJECT = (LB > 0 ) | (UB < 0)

  ## Do not reject if a cell count is zero
  REJECT[is.na(REJECT)]=FALSE
  POWER = c(POWER,mean(REJECT))
}

plot(seq(from=10, to=5000, by=100) , POWER, t="l", 
     xlab="Sample Size,n", ylab="POWER", main="")



## The Same thing as above, but we use modified 
## code from lecture notes

reps=1000
##  Population Cell probs
P=(1:4)/sum(1:4)

POWER=NULL
for ( n in seq(from=10, to=5000, by=100) )
{
  ## Array of rejection indicators.
  REJECT = array(0, reps)
  for (k in 1:reps) 
  {
    ## Generate a single table
    N = simtab(P, n)
    ## The sample log-odds ratio and its standard error.
    LR = log(N[1,1]) + log(N[2,2]) - log(N[1,2]) - log(N[2,1])
    SE = sqrt(1/N[1,1] + 1/N[1,2] + 1/N[2,1] + 1/N[2,2])
    
    ##  The 95 Percent CI
    LB = LR - 2*SE
    UB = LR + 2*SE

    ## Check to see if Zero is outside the interval.
    if (!is.finite(SE)) { REJECT[k] = 0 }
    else { REJECT[k] = (LB > 0 ) | (UB < 0) }
  }
  POWER=c(POWER, mean(REJECT))
}


plot(seq(from=10, to=5000, by=100), POWER, t="l", 
     xlab="Sample Size,n", ylab="POWER", main="")