## Simulate a 2x2 table with sample size n and cell probabilities 
## given in P.
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)
}

## The sample size.
n = 50

## Array of coverage indicators.
C = array(0, 1000)

for (k in 1:1000)
{
  ## Generate the four cell probabilities (p11, p12, p21, p22).
  P = 0.1 + 0.8*runif(4)
  P = P / sum(P)

  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 population log-odds ratio.
  PLR = log(P[1]) + log(P[4]) - log(P[2]) - log(P[3])

  ## Check for coverage.
  if (!is.finite(SE)) { C[k] = 1 }
  else { C[k] = (LR-2*SE < PLR) & (LR+2*SE > PLR) }
}