Monte Carlo Error In Posterior Quantiles

Why Quantile Precision Needs Its Own Check

Posterior quantiles are often reported as interval endpoints for a treatment effect. But a large total number of retained draws does not, by itself, show that an endpoint is accurately computed. MCMC draws may be serially dependent, and the dependence relevant to a particular quantile is the dependence of whether draws fall below that quantile. The lower and upper endpoints can therefore have different Monte Carlo precision, even when they come from the same simulation.

The earlier material on effective sample size and Monte Carlo standard errors introduced the effect of serial dependence on simulation accuracy. Here the target is specifically a posterior quantile. Its Monte Carlo error is the repeated-simulation variation in the estimated quantile when the same posterior is sampled again with a new MCMC run. It is distinct from posterior uncertainty: the posterior distribution describes uncertainty about the parameter given the data, while Monte Carlo error describes numerical uncertainty in approximating a feature of that distribution.

Definition: Let $\theta$ have posterior cumulative distribution function $F$. For $0<p<1$, the posterior $p$-quantile is $q_p=\inf\{x:F(x)\ge p\}$. Given $N$ retained draws $\theta_1,\ldots,\theta_N$, the estimated quantile $\hat q_p$ is the corresponding empirical quantile. Its Monte Carlo standard error (MCSE) estimates the standard deviation of $\hat q_p$ over repeated simulations targeting the same posterior.

For a continuous posterior with positive density near $q_p$, the precision of $\hat q_p$ depends on two things: how well the simulation estimates the posterior mass below $q_p$, and how much posterior density lies near the quantile. A high local density means a small change in probability corresponds to a small change in the quantile. A low local density means that the same probability error can move the quantile farther.

Quantile-Specific Effective Sample Size

Define a threshold indicator for each draw, $I_t=1\{\theta_t\le q_p\}$. At the true quantile, its posterior mean is $p$, and its variance is $p(1-p)$, provided the posterior is continuous at $q_p$. Serial dependence in this indicator determines how precisely the simulation estimates the fraction of draws below the quantile.

Let $\gamma_k$ be the lag-$k$ autocovariance of the stationary indicator sequence, and let $\sigma_p^2=\gamma_0+2\sum_{k=1}^{\infty}\gamma_k$ be its long-run variance. When a central limit theorem applies, the variance of the estimated fraction below $q_p$ is approximately $\sigma_p^2/N$. This motivates the quantile-specific effective sample size:

\operatorname{ESS}_p = \frac{N p(1-p)}{\sigma_p^2}.

Equivalently, if the indicator autocorrelations are $\rho_k$, then $\sigma_p^2=p(1-p)\{1+2\sum_{k=1}^{\infty}\rho_k\}$, so $\operatorname{ESS}_p$ is approximately $N/\{1+2\sum_{k=1}^{\infty}\rho_k\}$. In practice, $q_p$ is unknown and the indicator is formed using $\hat q_p$; autocorrelation and long-run variance must also be estimated from finite output. Thus the reported ESS is an estimate, not a known count.

Result: If the posterior has a continuous density $f$ that is positive and regular near $q_p$, and the MCMC indicator sequence satisfies an appropriate central limit theorem, then $$ \operatorname{MCSE}(\hat q_p) \approx \frac{\sqrt{\sigma_p^2/N}}{f(q_p)} = \sqrt{\frac{p(1-p)}{\operatorname{ESS}_p f(q_p)^2}}. $$ The density $f(q_p)$ is measured in inverse units of $\theta$, so the MCSE is in the same units as $\theta$.

This expression requires an estimate of the posterior density near the quantile, as well as an estimate of the indicator’s long-run variance. Common practical approaches estimate the density locally from posterior draws and estimate the long-run variance using a spectral or batch-based method. The estimate should be examined for sensitivity to reasonable choices of bandwidth or truncation: a noisy density estimate can dominate the resulting MCSE.

Why the Formula Works

The key approximation is that a small error in the empirical cumulative distribution function (CDF) near $q_p$ translates into a quantile error inversely proportional to the local density. Write $F_N$ for the empirical CDF. Under local smoothness, a first-order expansion around $q_p$ gives $F(\hat q_p)-F(q_p)\approx f(q_p)(\hat q_p-q_p)$. Since $F(q_p)=p$, the empirical quantile satisfies $F_N(\hat q_p)\approx p$, and therefore:

\hat q_p-q_p \approx -\frac{F_N(q_p)-p}{f(q_p)}.

The numerator is the simulation error in the estimated fraction of draws at or below $q_p$. For an autocorrelated stationary sequence, its variance is approximately the long-run variance divided by $N$. Dividing that variance by $f(q_p)^2$ yields the quantile variance formula above. This justification relies on more than a large $N$: it needs quantile consistency, a positive and sufficiently smooth local density, and a valid central limit approximation for the indicator sequence.

Worked Example: Quantile ESS for a Trial Treatment Effect

Suppose a hypothetical Bayesian trial model produces $N=20{,}000$ retained draws of a treatment log-odds effect $\theta$. The estimated 95th percentile is $\hat q_{0.95}=1.20$ log-odds units. For the indicator of draws below this quantile, suppose estimated autocorrelations are $0.40$ at lag 1 and $0.20$ at lag 2, with later autocorrelations negligible. A local density estimate is $0.25$ per log-odds unit. Use these approximations to estimate the quantile ESS and MCSE.

The integrated autocorrelation factor for the indicator is $1+2(0.40+0.20)=2.20$. Thus:

\operatorname{ESS}_{0.95} \approx \frac{20{,}000}{2.20} = 9{,}091.

Using $p(1-p)=0.95(0.05)=0.0475$, the estimated MCSE is:

\operatorname{MCSE}(\hat q_{0.95}) \approx \sqrt{\frac{0.0475}{9{,}091(0.25)^2}} = \sqrt{\frac{0.0475}{568.19}} \approx 0.00914

The quantile’s estimated Monte Carlo standard error is about $0.00914$ log-odds units. Under a roughly normal repeated-simulation approximation, a 1.96-MCSE range for simulation error is about $1.20\pm1.96(0.00914)$, or $1.20\pm0.0179$ log-odds units. This describes numerical uncertainty in the estimated 95th percentile; it is not a posterior credible interval for the treatment effect.

Repeated-Batch Quantile Estimates

A second practical approach divides a sufficiently long, approximately stationary simulation into equal-length batches, calculates the target quantile separately in each batch, and examines the variation among those estimates. If batch estimates are approximately independent and have similar simulation variance, their sample standard deviation divided by the square root of the number of batches estimates the Monte Carlo standard error of their mean. With long batches, the mean of the batch quantiles and the quantile from all draws target the same population quantile and are asymptotically close. In finite samples, they need not be exactly equal.

Worked Example: Checking a Quantile With Eight Batches

In a hypothetical trial analysis, eight equal-length batches of retained draws produce 90th-percentile estimates for a treatment log-odds effect of $0.31,\ 0.34,\ 0.33,\ 0.36,\ 0.32,\ 0.35,\ 0.34,$ and $0.33$. Estimate the standard error across the batches’ quantile estimates.

Their mean is $2.68/8=0.335$ log-odds units. Expressing deviations in hundredths of a log-odds unit gives $-2.5,\ 0.5,\ -0.5,\ 2.5,\ -1.5,\ 1.5,\ 0.5,$ and $-0.5$. These deviations sum to zero, and their squared sum is $18$ hundredths-squared, or $0.0018$ squared log-odds units. The sample variance and standard deviation are:

s_{\mathrm{batch}}^2 = \frac{0.0018}{8-1} \approx 0.00025714, \qquad s_{\mathrm{batch}} \approx 0.01604.

If the eight batch estimates are sufficiently independent and have comparable variance, the estimated standard error of their mean is:

\frac{s_{\mathrm{batch}}}{\sqrt{8}} = \frac{0.01604}{\sqrt{8}} \approx 0.00567

This gives an estimated Monte Carlo standard error of about $0.00567$ log-odds units for the average of the batch quantile estimates. The calculation depends on the batches being long enough to reduce within-batch dependence and the run being stable; eight batch estimates alone do not guarantee either condition.

Interpreting the Estimates and Checking for Failure

Quantile-specific ESS and repeated batches answer related but not identical practical questions. Quantile ESS isolates the serial dependence of the threshold indicator and, together with a local density estimate, produces an MCSE directly in parameter units. Batch quantiles provide a relatively direct view of between-batch variation, but their precision depends on the number and length of batches. Agreement between approaches can be reassuring; disagreement calls for investigation rather than choosing the smaller MCSE.

Before interpreting either estimate, check that the retained chains adequately explore the posterior and that the quantile estimate is not being driven by an unsettled chain or a poorly explored mode. Inspect trace plots and chain-to-chain behavior, and examine batch quantiles in their order through the simulation. A gradual drift or clusters of batch estimates can signal that stationarity or adequate mixing is doubtful. In that case, an MCSE computed under a stationary-chain approximation may understate the actual numerical uncertainty.

Other warning signs include a quantile ESS much smaller than the total number of draws, unstable indicator autocorrelation estimates, a local density estimate that changes substantially with smoothing choices, and batch MCSEs that change markedly when batch length changes. Tail quantiles can be particularly difficult: there is less posterior mass beyond a tail endpoint, and the local density may be small. More draws can help, but only if the additional simulation explores the relevant posterior region.

Worked Example: Independent Draws and a Median MCSE

Suppose a hypothetical trial analysis uses $N=10{,}000$ independent posterior draws of a treatment-effect parameter. At its median, the density is estimated as $0.10$ per parameter unit. For the indicator below the median, suppose the estimated long-run variance is $0.50$. Calculate the quantile ESS and MCSE.

For a median, $p=0.50$, so $p(1-p)=0.25$. The quantile-specific ESS is:

\operatorname{ESS}_{0.50} = \frac{10{,}000(0.25)}{0.50} = 5{,}000.

The corresponding MCSE is:

\operatorname{MCSE}(\hat q_{0.50}) \approx \sqrt{\frac{0.50}{10{,}000}}\frac{1}{0.10} = \frac{\sqrt{0.00005}}{0.10} \approx \frac{0.007071}{0.10} = 0.07071

The estimated MCSE is $0.07071$ parameter units. The stated independence is not consistent with a long-run variance of $0.50$, because independent median-indicator draws have variance $0.25$ and zero nonzero-lag autocovariances. Thus this numerical calculation is an illustration of the formula with the supplied long-run variance, not a coherent independent-draw scenario. For independent draws in this example, the long-run variance would be $0.25$, yielding $\operatorname{ESS}_{0.50}=10{,}000$ and MCSE $\sqrt{0.25/10{,}000}/0.10=0.0500$ parameter units.

Assumptions, Limitations, and a Common Mistake

The standard MCSE approximation is most defensible when the retained simulation is stationary and representative of the target posterior, the relevant indicator sequence has a stable long-run variance, and a central limit approximation is reasonable. The posterior CDF must also be locally smooth at the target quantile with $f(q_p)>0$. For discrete or mixed posteriors, ties and jumps in the CDF can invalidate the density-based approximation; quantile behavior then needs methods appropriate to that distribution rather than a mechanically applied continuous formula.

Density estimation is a separate source of approximation error. A reported quantile MCSE often treats the estimated density as if it were accurate enough, but a noisy or bandwidth-sensitive density estimate can make the MCSE itself unreliable. Repeated batches avoid explicit density estimation, yet still require an adequately long run and enough batches to estimate between-batch variation. With few batches, the sample standard deviation is itself uncertain. Neither approach compensates for modes the sampler failed to visit or nonstationary output.

Mistake: using total draws as the precision argument. “There are 50,000 draws, so the 95th percentile must be precise” ignores dependence and local posterior density. Report a quantile MCSE or assess batch quantile stability instead.
Mistake: using an ESS for the parameter mean as the quantile ESS. Dependence in $\theta_t$ and in $1\{\theta_t\le q_p\}$ need not be the same. Estimate dependence for the threshold indicator associated with the quantile of interest.
Mistake: interpreting MCSE as posterior uncertainty. An MCSE describes computational variability in an estimate across repeated simulations; it is not the width of a credible interval or a probability statement about the fixed quantile.

Key takeaway: For a posterior quantile, total draws are not a precision diagnostic. Estimate the Monte Carlo error for the quantile itself—using the threshold indicator’s long-run variance and local posterior density, or the variation among suitable repeated batches—and verify that the chain and regularity assumptions are credible.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why can a 99th-percentile estimate have a larger MCSE than a median estimate from the same simulation?
What is the quantile-specific ESS for independent draws, and why?
Using $p=0.50$, $N=10{,}000$, $\sigma_p^2=0.50$, and $f(q_p)=0.10$ per unit, calculate the quantile MCSE.
Name two conditions needed for the density-based quantile MCSE approximation, and give one reason it may fail for a discrete posterior.

Solutions. (1) The 99th-percentile indicator may have stronger serial dependence, and the posterior density near the 99th percentile may be lower. A larger long-run variance or smaller density increases the MCSE. (2) It is $N$: for independent draws all nonzero-lag autocovariances vanish, so the indicator long-run variance is $p(1-p)$, giving $\operatorname{ESS}_p=N$. (3) $\sqrt{0.50/10{,}000}/0.10=\sqrt{0.00005}/0.10\approx0.007071/0.10=0.07071$ units. (4) Conditions include a positive, sufficiently smooth density near the quantile and a valid central limit approximation for the indicator sequence. A discrete posterior can have a jump in its CDF, ties, or no positive continuous density at the quantile, so the local density approximation need not apply.