Monte Carlo, Importance Sampling, and MCMC

The sampling mechanism determines the estimator and its diagnostics. Each of the four branches below fails in its own way, so each is checked with its own number: a trial budget for direct simulation, a weight-tail parameter for reweighting, an autocorrelation summary for a dependent chain, an ensemble size for a population of particles. None of those numbers is a matter of taste. Each has been written down by a standards body, a journal or a review, and the rest of this lesson quotes them.

Many expectations and probabilities cannot be evaluated analytically. Monte Carlo methods approximate them using samples from a distribution that is easy enough to generate. The method is general because the error rate does not depend directly on geometric dimension in the same way as grid integration. It is also slow. Root-n convergence is slow enough that an international measurement standard had to legislate a number for it.

The metrology standard JCGM 101:2008 fixes the sample budget outright. Its own text reads: “A value of M = 10^6 can often be expected to deliver a 95 % coverage interval for the output quantity such that this length is correct to one or two significant decimal digits.” A million draws buys one or two digits of the interval's length. That is the price of the square root.

The standard also gives a general rule rather than a single default. For a coverage probability p it requires “a value of M that is large compared with 1/(1-p), e.g. M at least 10^4 times greater than 1/(1-p)”, and it allows M to be selected adaptively as the trials progress, because no pre-assigned number is guaranteed to suffice. The UK's National Physical Laboratory built software to the standard in 2011. It reproduces the four Supplement 1 example problems at M = 10^6 and applies the procedure “either non-adaptively or adaptively”.

Good Monte Carlo work therefore depends on variance reduction, diagnostics, and honest accounting of dependence and rare events.

The rare event does not have to be invented. FAA Advisory Circular 25.1309-1B, issued 30 August 2024, defines the top severity category: “Extremely Improbable Failure Condition. A failure condition having an average probability per flight hour on the order of 1 x 10-9 or less.” That is the target a certification argument has to reach.

Now try to reach it by sampling. Butler and Finelli, at NASA Langley Research Center, did the arithmetic in 1993 in IEEE Transactions on Software Engineering. Take a system required to fail with probability 10^-9 over a 10-hour mission. Their Table 1 gives an expected test duration of 10^10 hours — 1,141,550 years — with a single test article, and 10^9 hours, or 114,155 years, with ten replicates. Their text calls life testing at these targets “prohibitively impractical”. Their abstract states the conclusion without hedging: “This paper affirms that the quantification of life-critical software reliability is infeasible using statistical methods whether applied to standard software or fault-tolerant software.” Naive sampling here is not merely inefficient. It is off the table by six orders of magnitude in wall-clock centuries.

Variance can often be reduced without changing the estimand. Three of the four strategies below leave the sampling distribution alone and work on the arithmetic around it. The fourth, importance sampling, changes the distribution itself. That is why its failure is the one with published numbers attached: a fitted tail parameter with a cut-off at 0.7, and an ensemble requirement that grows exponentially with dimension. The two sections that follow supply both.

An expectation from random computation earns trust through this protocol. Two of its six steps rest on published thresholds rather than on judgement, and one of those two is skipped by most of the literature that would need it.

Step 4, inspect weights, has a pass/fail statistic. Pareto smoothed importance sampling fits a generalised Pareto distribution to the tail of the importance ratios and reads off the shape parameter, k̂. Vehtari and four co-authors set the method out in the Journal of Machine Learning Research in 2024. The algorithm's final step reports a warning “If the estimated shape parameter k̂ exceeds min(1 - 1/log10(S), 0.7)”, and above two thousand samples the formula collapses to a single number: “For S > 2000 this threshold is 0.7”. Beyond it the estimates are unstable, biased, and their Monte Carlo error is itself unreliable. Marc Box, in his own paper on importance sampling with partitioned data, summarises the scale: “They find that k̂ < 0.5 is an indicator of good performance, but that Pareto smoothed importance weights will still provide reliable results for k̂ < 0.7 and that importance sampling is unreliable beyond this.” Step 4 is a number to report. It is not an instruction to look at a histogram and feel reassured.

Step 5, estimate Monte Carlo error, is the one that gets left out. Morris and colleagues, writing in Statistics in Medicine in 2019, read 100 randomly selected simulation studies from volume 34 of that same journal and counted: “In our review of Volume 34, seven articles presented Monte Carlo SEs for estimated performance: three in the text, two in a table, one in a graph, and one in a float caption.” Seven of a hundred. The budget was not defended either: “Of the 87 studies reporting n_sim, four provided any justification of the choice.”

An independent audit in a different field reached the same place. Six authors screened 321 articles from three psychology journals published in 2021-2022 and found a simulation study in 100 of them. Their report, in Psychological Methods in November 2024, says that “The vast majority of studies (77%) did not report the uncertainty of performance measures (Panel B in Figure 2), despite our liberal approach of including visualizations such as box plots as indicative of Monte Carlo uncertainty.” Only 8% justified the number of repetitions.

The protocol below is a correction of documented practice. It is not a description of what everyone already does.

When direct samples are unavailable, a Markov chain can be designed so that its stationary distribution is the target. After convergence, chain states can estimate expectations. The samples are dependent, so the raw number of iterations overstates the information they carry.

The practice of checking this has a date. Gelman and Rubin introduced the potential scale reduction factor in Statistical Science in November 1992, together with the design that makes it meaningful. Their abstract puts it plainly: “Our recommended strategy is to use several independent sequences, with starting points sampled from an overdispersed distribution.” A single chain cannot report on the territory it never left.

The threshold that came with it has since been tightened by an order of magnitude, and on simulation evidence rather than taste. Vehtari and four co-authors, in Bayesian Analysis in 2021, recommend “only using the sample if R̂ < 1.01. This threshold is much tighter than the one recommended by Gelman and Rubin (1992), reflecting lessons learnt over more than 25 years of use”. They put a floor under the effective sample size in the same breath: “To ensure reliable estimates of variances and autocorrelations needed for R̂ and ESS, we recommend requiring that the rank-normalized ESS is greater than 400, a number we chose based on practical experience and simulations (see Appendix A) as typically sufficient to get a stable estimate of the Monte Carlo standard error.”

Both numbers are necessary and neither is sufficient. Rank-normalised R̂ below 1.01 and rank-normalised ESS above 400 are what current practice asks before a chain's output is used at all. Convergence diagnostics can reveal problems, but they cannot prove that every relevant mode was explored.