Why Validate MCMC Against an Exact Answer?
MCMC is useful when a posterior distribution cannot be calculated in closed form, but that flexibility makes implementation errors harder to spot: a chain can look stable while targeting the wrong distribution. A conjugate model supplies a controlled test. Its posterior is available exactly, so the simulated output can be compared with known means, variances, and quantiles. This checks whether the implementation reproduces a target distribution in a setting where the answer is not in doubt.
The comparison must account for Monte Carlo error. A simulated posterior mean will almost never equal its exact value digit for digit, even when the sampler is correct. Instead, ask whether the discrepancy is plausible given the effective amount of simulation. A very large discrepancy calls for investigation; a small one is evidence of agreement, not proof that every part of a more complex model is correct.
Set Up an Exact Reference
For a beta-binomial model, suppose \(Y\mid p\sim\operatorname{Binomial}(n,p)\), where \(Y\) is the number of successes in \(n\) trials and \(p\) is the success probability. Let the prior be \(p\sim\operatorname{Beta}(a,b)\), with \(a,b>0\). If the observed count is \(y\), the likelihood is proportional to \(p^y(1-p)^{n-y}\). Multiplying by the prior density gives a posterior proportional to \(p^{a+y-1}(1-p)^{b+n-y-1}\), so
For a \(\operatorname{Beta}(A,B)\) distribution, the exact posterior mean and variance are
where \(A=a+y\) and \(B=b+n-y\). These formulas assume the stated binomial sampling model and beta prior. For validation, use precisely the same observed data and prior in both the exact calculation and the code being tested.
A second useful test is a normal model with known sampling variance. Suppose \(X_i\mid\theta\sim N(\theta,\sigma^2)\) independently for \(i=1,\ldots,n\), and \(\theta\sim N(\mu_0,\tau_0^2)\), with known positive \(\sigma^2\) and \(\tau_0^2\). Completing the square in the product of prior and likelihood yields
The independence assumption and known sampling variance are part of this model. The normal-conjugate test is especially useful for checking that the implementation has the intended posterior center and spread, including on an unconstrained parameter scale.
Worked Example: Checking a Beta-Binomial Posterior
Suppose a hypothetical trial arm has 8 responses among 20 participants. Use a \(\operatorname{Beta}(2,2)\) prior for the response probability \(p\). Assume participants’ binary outcomes are conditionally independent with a common response probability, as required by the binomial model.
The exact posterior is \(\operatorname{Beta}(2+8,2+20-8)=\operatorname{Beta}(10,14)\). Its mean is
and its variance is
Now suppose a hypothetical chain, after warmup and with diagnostics otherwise satisfactory, has effective sample size 10,000 for \(p\). It gives a posterior mean of 0.4178 and variance of 0.00965. Using the exact variance as an approximation for the variance of posterior draws, the mean Monte Carlo standard error is
The mean discrepancy is \(0.4178-0.4167=0.0011\), about 1.1 mean Monte Carlo standard errors. That is a plausible simulation difference. The variance estimate is also close to the exact 0.009722; its comparison should use a variance-specific Monte Carlo error estimate, not the mean’s standard error. This check supports agreement for these summaries, but does not by itself verify every likelihood or decision calculation in a larger model.
Compare Each Summary on Its Own Error Scale
For \(M\) independent draws from a posterior with variance \(V\), the sample mean has standard error \(\sqrt{V/M}\). MCMC draws are generally autocorrelated, so the raw number of saved iterations is not the effective number of independent draws. Use an appropriate effective sample size or a chain-based Monte Carlo standard error (MCSE), as covered in the earlier material on effective sample size and Monte Carlo standard errors. If the MCSE is available directly, prefer it to treating all saved draws as independent.
Variance estimates also have simulation error. For independent draws from a distribution with variance \(V\) and fourth central moment \(\mu_4\), the large-sample variance of the sample variance is approximately \((\mu_4-V^2)/M\). Thus the corresponding standard error is approximately \(\sqrt{(\mu_4-V^2)/M}\). For dependent draws, this independent-draw expression does not account for autocorrelation; use a suitable chain-based estimate or an effective-sample approximation with care. A normal approximation \(\sqrt{2V^2/M}\) is available when the draws are normal, since then \(\mu_4=3V^2\). It should not be assumed for an arbitrary posterior.
For a posterior quantile \(q_p\) at probability level \(p\), the large-sample independent-draw approximation is
where \(f(q_p)\) is the posterior density at the exact quantile. The approximation requires a positive, sufficiently smooth density near the quantile. For an MCMC chain, dependence must also be accounted for; an effective sample size for the quantile or a chain-based quantile MCSE is more appropriate than automatically substituting the total iteration count.
Worked Example: Checking a Normal-Conjugate Posterior
Consider a hypothetical model with \(\theta\sim N(0,4)\), known sampling variance \(\sigma^2=1\), \(n=3\) independent observations, and observed sample mean \(\bar{x}=1\). The assumptions are a normal sampling model with common known variance and the stated normal prior. The posterior variance and mean are
Suppose a hypothetical MCMC run has effective sample size 5,000 for \(\theta\), with estimated mean 0.9250 and variance 0.3120. The mean MCSE, using the exact posterior variance, is
The mean differs from the exact value by approximately \(0.0019\), much less than one MCSE. Since the exact posterior is normal, the approximate standard error of the variance estimate is
The variance discrepancy, \(0.3120-0.3077=0.0043\), is also less than one such standard error. These comparisons are consistent with the exact target. They are approximate error checks, not formal acceptance tests with a universal pass threshold.
Worked Example: Checking a Posterior Quantile
For the \(\operatorname{Beta}(2,2)\) distribution, the density is \(f(p)=6p(1-p)\), and its cumulative distribution function is \(F(p)=3p^2-2p^3\). Solving \(F(q)=0.025\) gives a lower 2.5% quantile of approximately \(q=0.0943\). The density there is
Suppose the quantile estimate comes from 4,000 independent-equivalent draws, with dependence properly reflected in that effective amount of simulation. Its approximate MCSE is
If the chain estimates this quantile as 0.0980, the difference from the exact value is \(0.0980-0.0943=0.0037\), less than one estimated MCSE. This is compatible with ordinary simulation variation. The calculation depends on the density at the quantile: using a density that is too small would overstate the quantile’s MCSE, while a density that is too large would understate it.
A Practical Validation Sequence
A useful check is deliberately small and transparent. First calculate the exact posterior independently of the MCMC code. Then run the sampler on exactly that model and compare more than one feature of the output. A discrepancy in the mean with agreement in the variance can point to a location or coding issue; a correct mean with the wrong spread can reveal a scale, prior, or likelihood error. Quantiles probe distributional shape that matching only the first two moments can miss.
Write the likelihood, prior, data, and parameterization explicitly, and record the exact posterior distribution and summaries.
Check that the code uses the same data, prior hyperparameters, likelihood convention, and parameter scale as the analytic calculation.
Use Monte Carlo standard errors that reflect autocorrelation; do not treat every stored iteration as an independent draw.
Assess means, variances, and selected quantiles relative to their own simulation errors, and investigate large or systematic discrepancies.
Re-run the known-answer test after changing samplers, transformations, model code, or summary calculations.
This comparison complements, rather than replaces, chain diagnostics such as trace plots, rank-normalized \( \hat{R} \), effective sample size, and checks for divergent transitions. A chain can agree with an exact answer by chance on one summary while still mixing poorly; conversely, excellent-looking diagnostics do not establish that the coded target is the intended posterior.
Failure Modes and Alternatives
When a discrepancy is too large to explain by the estimated MCSE, first verify that the analytic calculation and code refer to the same target. Common causes include reversing success and failure counts, omitting a prior contribution, using a variance where a standard deviation is required, transforming draws but comparing them on the wrong scale, or summarizing a chain before warmup has been removed. If the target matches, inspect mixing and error estimates: poor exploration can make the apparent MCSE misleadingly small.
A conjugate model is a known-answer test, not a general guarantee. It may not expose bugs that occur only with multiple parameters, hierarchical structure, nonconjugate likelihoods, or a particular boundary region. Nor does matching a few summaries prove that two distributions are identical. More summaries, graphical comparisons, repeated runs, and tests across more than one simple model provide stronger evidence.
| Approach | Strength | Limitation |
|---|---|---|
| Exact conjugate comparison | Provides reference values without simulation uncertainty in the target summaries | Applies only to the specified test model; cannot certify a complex model on its own |
| Simulation-based checks without an exact answer | Can test broader model code and repeated-sampling behavior | Does not provide the same direct, known posterior values for comparison |
| Chain diagnostics alone | Assess aspects of mixing and sampling behavior | Do not establish that the implementation targets the intended posterior |
Worked Example: Detecting a Count-Coding Error
Suppose a hypothetical validation data set has 3 responses among 12 participants and uses a \(\operatorname{Beta}(1,1)\) prior. Under the stated binomial model, the exact posterior is \(\operatorname{Beta}(4,10)\), with mean \(4/14=0.2857\). Assume the MCMC implementation has effective sample size 5,000 and returns a mean of 0.3330. Its mean MCSE under the exact target is approximately
The discrepancy is \(0.3330-0.2857=0.0473\), or approximately \(0.0473/0.001650=28.7\) MCSEs. Ordinary simulation variation is not a plausible explanation. Check whether the code used the correct number of successes and failures, and whether the prior was included as specified. This calculation identifies a serious mismatch, though the exact source still requires debugging.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- With a \(\operatorname{Beta}(1,1)\) prior and 6 successes among 10 trials, what is the posterior distribution and its exact mean?
- Why should the raw number of saved MCMC iterations not automatically be used as \(M\) in an independent-draw MCSE formula?
- For a quantile, what role does the posterior density at that quantile play in its approximate Monte Carlo standard error?
- A simulated posterior mean is 0.020 from the exact value, with MCSE 0.010. Is the discrepancy conclusive evidence of an implementation error? Explain.
- What can a successful conjugate validation check establish, and what can it not establish?
Solutions. (1) The posterior is \(\operatorname{Beta}(1+6,1+10-6)=\operatorname{Beta}(7,5)\), with mean \(7/(7+5)=0.5833\). (2) MCMC draws are generally autocorrelated, so the stored iterations contain less information than the same number of independent draws; use an MCSE or effective sample size that accounts for dependence. (3) The quantile MCSE is approximately \(\sqrt{p(1-p)/M}/f(q_p)\); a lower density near the quantile means greater uncertainty in its estimated location. (4) The discrepancy is two MCSEs. It merits attention but is not by itself conclusive evidence of an error; consider the full pattern of summaries, chain diagnostics, and repeated-run behavior. (5) It provides evidence that the implementation reproduces the exact target in the tested model and parameterization, within simulation error. It does not prove correctness for more complex models, all summaries, or every possible coding path.