Monte Carlo Conditional Inference

When Conditional Enumeration Is Too Large

A conditional exact test calculates a tail probability under the null distribution given specified sufficient statistics or margins. If the conditional sample space is small, this probability can be found by summing over every table or outcome in that space. When the space is enormous, the same probability can instead be estimated by simulation. The result is then an approximation, and its Monte Carlo error must be reported or assessed separately from the uncertainty represented by the test itself.

Let \(\Omega_s\) be the set of possible outcomes consistent with the observed conditioning statistic \(S=s\). Under a specified null hypothesis, write \(q(y)=P_0(Y=y\mid S=s)\) for the conditional probability of \(y\in\Omega_s\). For a test statistic \(T\), define the inclusive upper-tail set \(A=\{y\in\Omega_s:T(y)\ge T(y_{\mathrm{obs}})\}\). The conditional exact p-value is the probability of this set under \(q\).

This definition concerns the probability under the conditional null law; it does not depend on whether that probability is computed by enumeration or simulation. The choice of tail, including the definition of \(T\), must be made as for the corresponding exact test. Simulation changes the accuracy of the calculation, not the test statistic or its null distribution.

Independent Conditional Simulation

Suppose \(Y_1,\ldots,Y_B\) are independent draws from \(q\), where \(B\) is the number of simulations. For each draw, record whether it falls in the tail set \(A\). Let \(I_b=1\) when \(Y_b\in A\), and \(I_b=0\) otherwise. The number of tail hits \(K=\sum_{b=1}^B I_b\) has a binomial distribution with parameters \(B\) and \(p\).

Proof. Each \(I_b\) is Bernoulli with \(P(I_b=1)=p\), so \(E[I_b]=p\) and \(\operatorname{Var}(I_b)=p(1-p)\). Independence gives \(E[K/B]=B^{-1}\sum_b E[I_b]=p\). It also gives \(\operatorname{Var}(K/B)=B^{-2}\sum_b\operatorname{Var}(I_b)=p(1-p)/B\), since the covariances between distinct indicators are zero. Replacing the unknown \(p\) in this variance formula with \(\hat p\) gives the stated estimated standard error. \(\square\)

For large \(B\), the central limit theorem for independent Bernoulli variables gives an approximate normal sampling distribution for \(\hat p\) when \(0<p<1\). This approximation can be poor when the expected number of hits \(Bp\) is small. In particular, the estimated standard error is zero when \(K=0\), but that does not establish that \(p=0\). A zero-hit run is evidence that the simulation has not resolved the tail at its chosen size.

Choosing a Simulation Size and Interpreting Zero Hits

A useful planning calculation follows directly from the variance formula. The standard error is largest over \(0\le p\le1\) at \(p=1/2\), where it is \(1/(2\sqrt B)\). Thus, \(B=10{,}000\) gives a worst-case standard error of \(0.005\). For a small tail probability, relative error is often more relevant: when \(p\) is small, the relative standard error is approximately \(1/\sqrt{Bp}\). To aim for relative standard error near 10% when \(p\) is about \(0.001\), one needs roughly \(Bp=100\), or \(B=100{,}000\) simulations.

If no tail events occur, a one-sided binomial upper bound gives more information than a zero plug-in standard error. With \(K=0\), the probability of this outcome at a candidate tail probability \(p\) is \((1-p)^B\). Setting this equal to \(0.05\) gives the 95% upper bound \(p=1-0.05^{1/B}\). For \(B=1000\), this is approximately \(0.00299\). Thus, zero hits in 1000 independent draws is compatible with a tail probability of nearly 0.003; it does not justify reporting an exact p-value of zero.

Monte Carlo Tests and the Plus-One Estimate

There is a useful distinction between estimating a fixed tail probability and constructing a simulation-based test. The proportion \(K/B\) estimates the conditional p-value for the observed data. For a test based on simulation, one can also count the observed outcome among the compared outcomes and use \((K+1)/(B+1)\). The added one is not a correction that makes \(K/B\) unbiased; rather, it yields a conservative simulation p-value when the observed outcome and simulations are exchangeable under the null.

Proof. Set \(m=B+1\). Because the \(m\) statistics are exchangeable, their ordering is equally likely to be any ordering after ties are broken uniformly at random. The resulting descending rank \(R^*\) of \(T(Y_0)\) is uniform on \(\{1,\ldots,m\}\). The inclusive rank \(R=1+\sum_{b=1}^B\mathbf{1}\{T(Y_b)\ge T(Y_0)\}\) is the largest rank that \(T(Y_0)\) could receive among its tied values. Therefore \(R\le k\) implies \(R^*\le k\). For \(k=\lfloor\alpha m\rfloor\), \(P(R/m\le\alpha)=P(R\le k)\le P(R^*\le k)=k/m\le\alpha\). Since \(p_{\mathrm{MC}}=R/m\), the claimed inequality follows. \(\square\)

This validity result assumes independent draws from the correct conditional null distribution, including the observed outcome in the exchangeability argument. It does not establish that a numerical estimate of a tiny p-value is precise. The value \(p_{\mathrm{MC}}\) also has a minimum of \(1/(B+1)\), so a simulation with limited \(B\) cannot report finer test resolution.

Rare Tails and Importance Sampling

When the tail is very rare, direct conditional simulation may produce few or no hits even with a substantial computational budget. One alternative is importance sampling: draw from a proposal distribution \(r\) that places more probability on the tail, then weight each draw by \(q(y)/r(y)\). The proposal must assign positive probability wherever \(q(y)\) contributes to the target tail. For independent draws \(Z_b\sim r\), the estimate \[ \hat p_{\mathrm{IS}}=\frac{1}{B}\sum_{b=1}^B \mathbf{1}\{Z_b\in A\}\frac{q(Z_b)}{r(Z_b)} \] is unbiased, because its expectation is \(\sum_{y\in A}r(y)q(y)/r(y)=\sum_{y\in A}q(y)=p\). Its variability depends on the proposal and weights; targeting the tail does not guarantee a precise estimate if a few draws receive enormous weights.

Markov-Chain Sampling and Mixing

Some conditional distributions are difficult to sample independently but can be explored by a Markov chain whose stationary distribution is \(q\). For example, a chain may move between tables in \(\Omega_s\) while preserving the required margins. If the state space is finite and the chain is irreducible and aperiodic, it has a unique stationary distribution; when initialized from that distribution, each state has marginal distribution \(q\). The draws are nevertheless dependent. After a suitable initial transient, averages can converge to the desired conditional probability, but their uncertainty is not given by the independent-binomial formula.

For a stationary chain \(Y_1,\ldots,Y_B\), set \(I_b=\mathbf{1}\{Y_b\in A\}\). If the chain has summable autocorrelations \(\rho_k=\operatorname{Corr}(I_b,I_{b+k})\), then \[ \operatorname{Var}\left(\frac{1}{B}\sum_{b=1}^B I_b\right) =\frac{\operatorname{Var}(I_1)}{B} \left[1+2\sum_{k=1}^{B-1}\left(1-\frac{k}{B}\right)\rho_k\right]. \] This follows by expanding the variance of the sum: there are \(B\) individual variances and \(2(B-k)\) covariance terms at lag \(k\), and each such covariance is \(\operatorname{Var}(I_1)\rho_k\). When the chain satisfies a Markov-chain central limit theorem, the large-\(B\) variance is approximately \(\operatorname{Var}(I_1)\tau/B\), where \(\tau=1+2\sum_{k=1}^{\infty}\rho_k\). The effective sample size is approximately \(B/\tau\); positive serial correlation makes it smaller than \(B\).

Practical Checks and Common Pitfalls

Simulation output should include enough information to judge how well the conditional tail has been resolved. For independent sampling, report \(B\), \(K\), the estimate, and an appropriate uncertainty assessment; if \(K\) is small, avoid relying on a normal approximation alone. For a Markov chain, report the chain length and assess dependence in the tail indicator, not only in the raw table entries. Multiple chains and trace plots can reveal some problems, but apparent stability does not prove that all relevant regions of the conditional space have been explored.

• Mistake: treating zero hits as a zero p-value. A finite run with no hits only shows that the event was not seen. Use a binomial upper bound or increase the simulation effort.
• Mistake: using \(K/B\) as if it had no simulation error. It is an estimate of a conditional probability, not the probability itself. State its Monte Carlo uncertainty.
• Mistake: applying an independent-sampling standard error to MCMC output. Correlated draws contain less information. Account for autocorrelation or effective sample size.
• Mistake: assuming more simulations automatically fix poor mixing. A chain that repeatedly explores only part of the relevant state space can remain misleading even when it produces many iterations.
• Mistake: treating importance weights as optional. Proposal samples estimate a different distribution unless corrected by the likelihood ratio.