Quantifying Monte Carlo Error in Posterior Probabilities

Why quantify simulation error?

A posterior decision rule may depend on a probability such as the posterior chance that an arm’s benefit exceeds a clinically meaningful threshold. In practice, that probability is often estimated from simulated posterior draws. The estimate is itself subject to simulation error: a different run of the simulation would generally give a slightly different value, even if the data, model, and posterior distribution were unchanged.

This error matters when the estimated probability is close to a decision cutoff. An estimate just above the cutoff is not persuasive evidence that the underlying posterior probability is above it if simulation noise is large. Monte Carlo standard errors (MCSEs) and intervals quantify this numerical uncertainty. They do not quantify uncertainty about the treatment effect itself; that uncertainty is already represented by the posterior distribution.

Estimating a posterior decision probability

Let \(y\) denote the observed data, and let \(A\) be the event that defines a decision-relevant posterior probability. For example, \(A\) might be the event that an arm’s treatment effect exceeds a specified benefit threshold. The target is \(p=P(A\mid y)\). Given \(M\) posterior draws, define \(I_m=1\) if draw \(m\) satisfies \(A\), and \(I_m=0\) otherwise. The simulation estimate is the average of these indicators.

If the posterior draws are independent, the indicators are independent Bernoulli variables with success probability \(p\). Their variance is \(p(1-p)\), so the variance of their average is \(p(1-p)/M\). Replacing the unknown \(p\) by \(\hat p\) gives the usual estimated MCSE:

This formula applies to genuinely independent posterior draws. MCMC draws are generally dependent, so treating their count as if it came from \(M\) independent draws can understate the error.

Accounting for MCMC dependence

For MCMC, the indicator sequence \(I_1,\ldots,I_M\) can be autocorrelated even when the chain has reached its target distribution. Let \(\gamma_k\) be the lag-\(k\) autocovariance of the indicator sequence, and let \(\rho_k=\gamma_k/\gamma_0\) be its lag-\(k\) autocorrelation. For a stationary chain, the variance of the sample average is

The bracketed factor accounts for dependence. When positive autocorrelation persists across lags, it is greater than one, so the variance exceeds the independent-draw variance \(\gamma_0/M\). For an indicator with marginal probability \(p\), \(\gamma_0=p(1-p)\). In practice, the autocovariances are estimated from the simulated indicator sequence, using a suitable spectral-variance or batch-means estimator.

A useful summary is the effective sample size for the indicator sequence, denoted here by \(ESS_I\). It is the number of independent Bernoulli draws that would give approximately the same variance as the correlated simulation. The corresponding approximation is

In applications, the infinite sum is estimated with a method that controls the effect of noisy high-lag autocorrelations. Use the ESS and MCSE for the decision indicator, not automatically the ESS for a treatment-effect parameter: the chain can mix differently for a threshold event than for a continuous parameter.

Intervals and threshold-based decisions

When the estimate is not too close to zero or one and the simulation is sufficiently regular, a normal approximation gives an approximate \(100(1-\alpha)\%\) Monte Carlo interval:

where \(z_{1-\alpha/2}\) is the corresponding standard normal quantile. For a 95% interval, it is approximately 1.96. The interval is an approximation to repeated-simulation uncertainty in \(\hat p\), conditional on the data and posterior target. It is not a posterior interval for an effect, and it is not a probability statement that the fixed target \(p\) has a particular chance of lying in this realized interval.

For a rule that compares \(p\) with a cutoff \(c\), the interval gives a direct precision check. If its lower endpoint is greater than \(c\), the simulation supports \(p>c\) at the approximate interval level; if its upper endpoint is less than \(c\), it supports \(p<c\). If the interval overlaps \(c\), the simulation has not resolved the comparison at that precision. That does not prove \(p=c\): more simulation may narrow the interval, or the posterior probability may genuinely be close to the threshold.

Choosing a precision target and checking the calculation

A useful simulation plan specifies a precision target before inspecting the final result. For example, one may require the Monte Carlo interval’s half-width to be no more than \(h\), or require the interval to lie entirely on one side of a decision cutoff. Under the independent-draw approximation, a rough effective sample size requirement for half-width \(h\) is

Because the required size depends on \(p\), planning may use a conservative value such as \(p=0.5\), which maximizes \(p(1-p)\), or a plausible value from pilot simulation. For MCMC, the target is an effective sample size, not merely a raw number of saved iterations.

As a short justification, the interval half-width is \(h=z_{1-\alpha/2}\sqrt{\hat p(1-\hat p)/ESS_I}\). Requiring it to be at most a chosen \(h\), squaring both nonnegative sides and rearranging gives the displayed ESS requirement. For example, with \(\hat p=0.70\), a 95% half-width target of 0.01 requires approximately

so an effective sample size of about 8,068 or greater is a rough target. This is a planning calculation, not a guarantee: the realized ESS estimate and chain diagnostics still need inspection.

Diagnostics, limitations, and alternatives

Before trusting an MCSE, check whether the simulation is credible for the event being estimated. The convergence and mixing diagnostics considered in the earlier material on diagnosing MCMC near decision boundaries remain relevant. In particular, inspect the indicator trace, compare chains, and check for long stretches in which the event status barely changes. A chain can appear well behaved for a parameter while rarely crossing the decision boundary, leaving the indicator ESS small and its MCSE unreliable.

• Do not use the raw iteration count as an independent sample size. Autocorrelation can make the effective information much smaller than \(M\).
• Do not report an MCSE without identifying what was monitored. The relevant series is the zero-one indicator of the decision event.
• Be cautious near zero or one. A normal interval can extend outside \([0,1]\), and if every simulated indicator is zero or every one, the plug-in MCSE is zero even though the true probability need not be exactly zero or one. More simulation and appropriate boundary-aware methods may be needed.
• Do not confuse precision with correctness. A tiny MCSE does not remedy a misspecified model, an unsuitable prior, a coding error in the event, or poor chain exploration.

For independent Bernoulli draws, a Wilson interval is one alternative to the simple normal interval, especially when the estimated probability is near a boundary. For MCMC draws, ordinary binomial intervals that treat every iteration as independent are generally inappropriate; dependence-aware variance estimation is needed. Increasing the number of iterations can improve precision, but only if the chain explores the relevant posterior region. The central practical comparison is therefore not “many draws versus few draws,” but “adequate effective information and diagnostics versus inadequate effective information.”