Monte Carlo Error in Estimated Allocation Probabilities

Why Monte Carlo Error Matters

A posterior best-arm probability is a probability under the current posterior distribution, conditional on the trial information observed so far. When it is evaluated by simulation, the reported value is only an estimate of that probability. This introduces Monte Carlo error: variation caused by the finite number of posterior draws, not by uncertainty about the treatment responses themselves.

That distinction matters in Thompson sampling. The allocation probability used for an arm may be approximated by the fraction of simulated posterior draws in which that arm is best. If this fraction is imprecise, the implemented randomization probability can differ from the intended posterior probability. The right number of draws depends on the precision needed for the decision; there is no single universally adequate simulation size.

Definition and Sampling Model

Let $D$ denote the trial information available at the allocation decision, and let $\theta_j$ be the response parameter for arm $j$, for $j=1,\ldots,K$. For a response where larger values are better, define the posterior best-arm probability by

q_j=P\left(\theta_j=\max_{1\leq \ell\leq K}\theta_\ell\mid D\right).

Assume the posterior distribution is continuous, so ties have probability zero. For each of $M$ independent draws from the joint posterior given $D$, record whether arm $j$ is best. The resulting indicator is $I_{mj}=1$ if arm $j$ is best on draw $m$, and $I_{mj}=0$ otherwise. The simulated estimate is

\hat q_j=\frac{1}{M}\sum_{m=1}^{M}I_{mj}.

Conditional on $D$, the indicators for a particular arm are independent Bernoulli random variables with probability $q_j$, provided the posterior draws are independent. Consequently, $\hat q_j$ is an unbiased estimate of $q_j$, with variance $q_j(1-q_j)/M$. The estimated Thompson sampling allocation probability is therefore subject to simulation error even though $D$ is being held fixed.

Definition: The Monte Carlo standard error of $\hat q_j$, under independent posterior draws, is the standard deviation of its simulation-to-simulation variation conditional on $D$. Its plug-in estimate is $$ \widehat{SE}_{MC}(\hat q_j)=\sqrt{\frac{\hat q_j(1-\hat q_j)}{M}}. $$

The notation $SE$ is often used for this standard error; the subscript $MC$ emphasizes its source. It measures numerical precision of the simulated probability, not posterior uncertainty in $\theta_j$, uncertainty about a treatment effect, or variation in patients’ outcomes.

Deriving the Standard Error

For a fixed arm, each $I_{mj}$ has mean $q_j$ and variance $q_j(1-q_j)$, because it takes the value 1 with probability $q_j$ and 0 otherwise. Independence across draws means the variance of their sum is the sum of their variances. Dividing by $M^2$, as required for the average, gives

\begin{aligned} E(\hat q_j\mid D)&=\frac{1}{M}\sum_{m=1}^{M}E(I_{mj}\mid D)=q_j,\\ \operatorname{Var}(\hat q_j\mid D) &=\frac{1}{M^2}\sum_{m=1}^{M}\operatorname{Var}(I_{mj}\mid D) =\frac{q_j(1-q_j)}{M},\\ SE_{MC}(\hat q_j\mid D)&=\sqrt{\frac{q_j(1-q_j)}{M}}. \end{aligned}

Since $q_j$ is the unknown quantity being estimated, replace it with $\hat q_j$ to estimate the standard error. For sufficiently large $M$, a normal approximation gives an approximate 95% Monte Carlo interval $\hat q_j\pm1.96\widehat{SE}_{MC}$. This is an interval for the numerical approximation to $q_j$, conditional on $D$; it is not a 95% posterior interval for a treatment response parameter.

The indicators across arms on the same draw are not independent: exactly one arm is best. Thus, $\sum_{j=1}^K\hat q_j=1$, and an increase in one estimated probability must be balanced by decreases in others. Under the no-tie assumption, for distinct arms $j$ and $\ell$, the covariance of the estimates is

\operatorname{Cov}(\hat q_j,\hat q_\ell\mid D)=-\frac{q_jq_\ell}{M}.

This dependence does not change the single-arm standard error formula, but it matters when assessing contrasts or making a simultaneous statement about several arms.

Worked Example: Standard Error for a Simulated Probability

In a hypothetical three-arm trial, an analyst draws $M=4{,}000$ independent sets of response parameters from the posterior. Arm A has the largest response parameter in 2,380 draws. Estimate its best-arm probability and Monte Carlo standard error.

The estimated probability is the fraction of draws in which A is best: $$ \hat q_A=\frac{2{,}380}{4{,}000}=0.595. $$

Substituting this estimate and the draw count into the plug-in formula gives $$ \widehat{SE}_{MC}(\hat q_A) =\sqrt{\frac{0.595(1-0.595)}{4{,}000}} =\sqrt{0.00006024} =0.007762. $$

The Monte Carlo standard error is about 0.0078, or 0.78 percentage points. A normal-approximation 95% Monte Carlo interval is $$ 0.595\pm1.96(0.007762) =0.595\pm0.01521 =(0.5798,\ 0.6102). $$

Conditional on the trial information, this interval describes approximate simulation precision for $q_A$. It does not say there is a 95% probability that the fixed posterior probability $q_A$ lies in the interval, and it does not quantify uncertainty about A’s response rate.

Choosing the Number of Draws

Because the standard error decreases in proportion to $1/\sqrt{M}$, reducing it by a factor of two requires approximately four times as many draws. For a desired approximate 95% margin of error $e$, a conservative single-arm planning calculation uses the largest possible Bernoulli variance, $q_j(1-q_j)\leq1/4$. It follows that

1.96\sqrt{\frac{1}{4M}}\leq e \quad\Longrightarrow\quad M\geq\frac{1.96^2}{4e^2}.

This normal-approximation calculation plans precision for one probability. If all $K$ probabilities must be accurate at once, a simple conservative bound can control the chance that any estimate differs from its target by more than $e$. Hoeffding’s inequality gives $P(|\hat q_j-q_j|\geq e\mid D)\leq2\exp(-2Me^2)$ for each arm. Applying the union bound across arms yields

P\left(\max_{1\leq j\leq K}|\hat q_j-q_j|\geq e\mid D\right) \leq 2K\exp(-2Me^2).

Therefore, it is sufficient to choose $M\geq\log(2K/\delta)/(2e^2)$ to make this upper bound no greater than $\delta$. This guarantee does not require knowing the $q_j$, but it can be substantially more conservative than necessary.

Worked Example: Planning for a Target Precision

Suppose a hypothetical analysis needs each of three estimated best-arm probabilities to be within 0.01 of its target with at least 95% probability, according to a conservative bound. Here $e=0.01$, $K=3$, and $\delta=0.05$. The required draw count satisfies $$ M\geq\frac{\log(2K/\delta)}{2e^2} =\frac{\log(120)}{2(0.01)^2} =\frac{4.7875}{0.0002} =23{,}937.5. $$

Rounding up, choose at least $23{,}938$ independent posterior draws. The result is a sufficient, not necessary, draw count for the stated simultaneous absolute-error guarantee. If only one arm’s probability needed approximate 95% precision of 0.01, the worst-case normal calculation would instead give $M\geq1.96^2/[4(0.01)^2]=9{,}604$. The difference reflects both the simultaneous requirement and the conservatism of the bound.

Assessing Rare and Unobserved Best Arms

For a small $q_j$, absolute Monte Carlo error may look modest while relative error remains large. The relative standard error is approximately $\sqrt{(1-q_j)/(Mq_j)}$, so it can be substantial unless the expected number of times the arm is best, $Mq_j$, is sufficiently large. In particular, a zero simulated count does not establish that the true posterior best-arm probability is zero.

Worked Example: A Rare Best-Arm Event

In a hypothetical simulation with $M=1{,}000$ independent posterior draws, arm C is best 18 times. Its estimated probability and standard error are $$ \hat q_C=\frac{18}{1{,}000}=0.018,\qquad \widehat{SE}_{MC}=\sqrt{\frac{0.018(0.982)}{1{,}000}}=0.004204. $$

The relative standard error is $0.004204/0.018=0.234$, or about 23.4%. Thus the estimate is 1.8%, but its Monte Carlo uncertainty is substantial relative to that estimate. A normal-approximation interval is particularly questionable for such a small count. More draws, or an interval method better suited to a binomial proportion, is appropriate if this probability is decision-relevant.

Now suppose that in a separate hypothetical run, arm C is best zero times in $M=500$ draws. The plug-in standard error is zero because $\hat q_C=0$, but that is not evidence of zero simulation uncertainty. Under independent draws, a one-sided exact 95% binomial upper limit after zero successes is $1-0.05^{1/500}=0.00597$. The simulation has not ruled out a best-arm probability of approximately 0.006.

Diagnostics, Alternatives, and Limitations

The standard error formula relies on independent posterior draws from the target posterior conditional on the current $D$. With direct independent draws, check that each draw uses the intended posterior and that best-arm indicators are recorded consistently. Verify that exactly one arm is counted as best on each draw; then the estimated probabilities should sum to one. Small departures from one can indicate a coding or bookkeeping error, not Monte Carlo noise.

If draws come from a Markov chain rather than independent sampling, successive indicators may be correlated. The binomial standard error above will generally understate uncertainty in that case. Use an effective sample size or a Monte Carlo standard error based on the autocorrelation of the indicator sequence, and inspect chain mixing. Similarly, repeatedly checking a conventional interval and stopping at the first time it is narrow enough can invalidate its nominal coverage. A fixed draw count avoids that particular repeated-checking issue; sequential stopping requires methods designed for sequential monitoring.

For moderate counts, the plug-in standard error and a normal approximation are convenient. When counts are small or zero, a Wilson or exact binomial interval is more informative than relying on a symmetric normal interval. These methods assess simulation uncertainty for each probability; a collection of separate 95% intervals does not automatically provide 95% simultaneous coverage across arms. The Hoeffding bound above is one conservative alternative when simultaneous absolute accuracy is required.

Increasing $M$ addresses simulation error only. It cannot correct a misspecified likelihood, an unsuitable prior, an incorrect response direction, or a bug in the posterior sampler. It also cannot remove posterior uncertainty arising from limited trial data. Computationally cheaper approximations may be useful, but they must be assessed against the desired precision and the allocation decision at hand. If two estimated allocation probabilities are close, compare their Monte Carlo uncertainty before interpreting their ordering as meaningful.

Common Mistake

Common mistake: “The simulated fraction is 0.60, so the allocation probability is exactly 0.60; increasing the number of draws is only needed to make the posterior more certain.” The correction is that 0.60 is a finite-simulation estimate of a fixed posterior probability conditional on $D$. More draws reduce its Monte Carlo error at a rate proportional to $1/\sqrt{M}$, but they do not reduce the posterior uncertainty about treatment responses. Those are distinct sources of uncertainty.

Key takeaway: For independent posterior draws, estimate a best-arm probability by the fraction of draws on which the arm is best, and estimate its Monte Carlo standard error by $\sqrt{\hat q_j(1-\hat q_j)/M}$. Choose $M$ for the required precision, and treat rare or zero counts with care.

Check Your Understanding

Answer each question, then compare with the solutions below.

In 2,500 independent draws, an arm is best 1,000 times. Calculate its estimated probability and Monte Carlo standard error.
By what factor must $M$ increase to halve the Monte Carlo standard error, approximately?
Why can the plug-in standard error be misleading when an arm is never best in the simulated draws?
What changes in the standard error calculation if posterior draws are autocorrelated?

Solutions. (1) $\hat q=1{,}000/2{,}500=0.400$, and $\widehat{SE}_{MC}=\sqrt{0.4(0.6)/2{,}500}=0.009798$, or about 0.98 percentage points. (2) Since the standard error is proportional to $1/\sqrt M$, $M$ must increase by a factor of four. (3) A zero estimate makes the plug-in formula return zero, although a finite run with no observed successes does not establish that the true probability is zero; use more draws or a binomial interval appropriate to a zero count. (4) Correlation means the independent-draw variance formula is generally too small; account for autocorrelation, for example through an effective sample size or an indicator-based Monte Carlo standard error.