Posterior Sampling and Convergence

A posterior sampler is a numerical integration method, not a certificate that the posterior has been explored. Gravitational-wave posteriors can be curved, multimodal, bounded, periodic, and strongly degenerate; evidence calculations add a global-normalization problem. This page compares the principal sampler families, defines diagnostics appropriate to weighted and correlated samples, and provides an executed importance-sampling benchmark whose exact mode probability and mean expose support loss and weight degeneracy.

Required background. Probability and Statistical Inference supplies Monte Carlo integration, conditional densities, Markov chains, importance weights, and sampling variance.

Helpful background. Gravitational Radiation and Waveforms supplies the parameter degeneracies and waveform likelihoods being sampled.

Sampling is integration

For normalized posterior

$$p(\theta\mid d) =\frac{\mathcal L(d\mid\theta)\pi(\theta)}{Z},$$

an expectation is

$$\mathbb E[g] =\int g(\theta)p(\theta\mid d)\,d\theta.$$

A sampler approximates this integral by correlated draws, weighted particles, or a sequence of constrained prior volumes. Its output is scientifically useful only when the approximation error is small for the quantities actually reported. Agreement of one median does not establish accurate tails, secondary modes, evidence, or joint credible regions.

Common families make different compromises:

• Markov-chain Monte Carlo constructs a correlated chain with the posterior as its stationary distribution. Efficiency depends on proposal geometry, adaptation, and mixing across modes.
• Nested sampling replaces low-likelihood prior mass by higher-likelihood points, estimating both posterior weights and evidence. Accuracy depends on constrained-prior sampling and resolving all important modes.
• Sequential Monte Carlo transports weighted particles through intermediate distributions, with resampling and rejuvenation. Temperature schedules and particle diversity are central.
• Importance sampling draws from a proposal $q$ and weights by the target-to-proposal ratio. It is parallel and diagnostic-rich, but fails when $q$ misses support or yields extreme weights.

The gravitational-wave use and validation of MCMC and nested sampling are described in Veitch et al. 2015, §§III–V, PDF; an introductory treatment of these samplers appears in Thrane and Talbot 2019, §4, PDF.

Effective samples are quantity dependent

For self-normalized importance samples $\theta_i\sim q$ with

$$w_i=\frac{\mathcal L(d\mid\theta_i)\pi(\theta_i)}{q(\theta_i)}, \qquad \widetilde w_i=\frac{w_i}{\sum_jw_j},$$

the estimate is

$$\widehat{\mathbb E[g]} =\sum_i\widetilde w_i g(\theta_i).$$

A common weight diagnostic is

$$N_\mathrm{eff}^{(w)} =\frac{1}{\sum_i\widetilde w_i^2}.$$

It ranges from one to $N$, but it does not detect proposal regions never sampled. A large nominal $N_\mathrm{eff}^{(w)}$ within one mode cannot reveal a missing second mode.

For independent proposal draws and finite importance-weight variance, a quantity-specific Monte Carlo variance estimate is

$$\widehat{\operatorname{Var}}\!\left(\widehat{\mathbb E[g]}\right) \simeq\frac{N}{N-1}\sum_i\widetilde w_i^2 \left[g(\theta_i)-\widehat{\mathbb E[g]}\right]^2.$$

This estimates sampling error of the weighted mean, not posterior uncertainty in $g$. The posterior spread remains $\operatorname{Var}_p(g)$ no matter how many samples are drawn.

For a stationary Markov chain and scalar summary $g$, an asymptotic effective size is

$$N_\mathrm{eff}^{(g)} \simeq \frac{N}{1+2\sum_{k\ge1}\rho_g(k)},$$

where $\rho_g(k)$ is the lag-$k$ autocorrelation of $g(\theta)$. Different parameters and tail indicators have different autocorrelation times. Rank-normalized split-chain diagnostics, between-chain comparisons, trace plots, and repeated initializations complement rather than replace effective size.

A bimodal importance-sampling benchmark

Take the normalized target

$$p(x) =0.7\,\mathcal N(x;-2,0.5^2) +0.3\,\mathcal N(x;2,0.5^2)$$

and proposal $q(x)=\mathcal N(x;0,3^2)$. The exact mean is

$$\mathbb E_p[x]=0.7(-2)+0.3(2)=-0.8.$$

The positive-half-line probability is

$$P_p(x\gt0) =0.7[1-\Phi(4)]+0.3\Phi(4) =0.3000127\ldots ,$$

essentially the right-component weight but with the tiny Gaussian tails included. With Python 3.9.6, NumPy==2.0.2, and seed 424242 reset for each row, the executed estimates are:

Proposal draws	$P(x\gt0)$ result	Mean result	Weight effective size
2000	0.302281	-0.815770	633.7
10000	0.300604	-0.810565	3296.6
50000	0.296815	-0.819760	16228.9

The estimates fluctuate around the exact values rather than moving monotonically toward them; only about one third of the nominal draws are weight-effective. At 50,000 draws, the quantity-specific standard errors are approximately $0.00304$ for $P(x\gt0)$ and $0.0127$ for the mean, placing the errors at about $1.05$ and $1.56$ standard errors respectively. This controlled success relies on the broad proposal covering both modes.

Executable rerun

import numpy as np

def normal_pdf(x, mean, sigma):
    z = (x - mean) / sigma
    return np.exp(-0.5 * z**2) / (sigma * np.sqrt(2 * np.pi))

for size in (2_000, 10_000, 50_000):
    rng = np.random.default_rng(424242)
    x = rng.normal(0.0, 3.0, size)
    target = (0.7 * normal_pdf(x, -2.0, 0.5)
              + 0.3 * normal_pdf(x, 2.0, 0.5))
    proposal = normal_pdf(x, 0.0, 3.0)
    weight = target / proposal
    weight /= np.sum(weight)
    positive_mass = np.sum(weight * (x > 0))
    mean = np.sum(weight * x)
    effective = 1.0 / np.sum(weight**2)
    print(size, f"{positive_mass:.6f}", f"{mean:.6f}",
          f"{effective:.1f}")

# Quantity-specific errors for the final 50,000-draw run.
for label, g in (("positive mass", x > 0), ("mean", x)):
    estimate = np.sum(weight * g)
    variance = size / (size - 1) * np.sum(
        weight**2 * (g - estimate)**2
    )
    print(label, "Monte Carlo SE", f"{np.sqrt(variance):.6f}")

Save this as importance_bimodal_check.py, then run:

python3 -m venv .venv
. .venv/bin/activate
python3 -m pip install numpy==2.0.2
python3 importance_bimodal_check.py

The benchmark should also be rerun with a narrow proposal centered at $-2$. It can return apparently stable left-mode summaries while assigning effectively no samples to the right mode. That deliberate failure demonstrates why convergence within visited support is insufficient.

A sampler validation ladder

Validation should proceed from cases with known answers to the target problem:

1. recover moments and evidence for Gaussian and low-dimensional analytic likelihoods;
2. recover injected parameters over repeated Gaussian-noise realizations;
3. compare independent seeds, dispersed starts, live-point counts, temperature schedules, or proposals;
4. inspect marginal and joint effective sizes, weights, traces, mode occupancy, and tail indicators;
5. compare a subset with a method based on different numerical assumptions; and
6. perturb waveform, PSD, calibration, prior bounds, and data segment to separate numerical from model sensitivity.

Stopping because the evidence increment is small does not guarantee accurate posterior tails. Conversely, a chain-based convergence statistic does not validate nested-sampling evidence. Diagnostics must match the estimator.

Event-analysis reproduction and failures

For a real event, record the GWOSC or controlled data release and hashes, event version, detectors and GPS interval, calibration and PSD products, quality selections and gates, sample rate and frequency band, waveform and likelihood commits, prior coordinate files, sampler implementation and version, parallelization, seeds, initialization, proposal or live-point settings, adaptation or temperature schedule, stopping rule, wall-clock precision, checkpoint policy, and post-processing weights. Retain logs and independent-run summaries.

At the evidence cutoff of 2026-08-18, the O4b analyses use versioned likelihoods, priors, and sampling configurations described in LVK Collaboration 2026, §§4–5, PDF. This page does not rerun those products. Its software claim is limited to the synthetic NumPy benchmark above.

Failure modes include a missing mode, poor proposal support, weight overflow, unrecognized chain nonstationarity, premature evidence stopping, too few independent likelihood evaluations, parameter-boundary sticking, incorrect Jacobians, stochastic likelihood noise, and post-processing that discards weights. Report repeated-run dispersion and any failed configuration rather than selecting only the most agreeable run.

Check your understanding

A chain has $N=10000$ retained draws. For the scalar $g(\theta)$ of interest, the estimated integrated autocorrelation factor is $1+2\sum_{k\ge1}\rho_g(k)=25$. Estimate the effective sample size. If $\operatorname{Var}_p(g)=1$, estimate the Monte Carlo standard error of the sample mean. Why might this still fail to detect a missing mode?

Solution

The effective size is

$$N_\mathrm{eff}^{(g)}\simeq\frac{10000}{25}=400.$$

The standard error is approximately

$$\operatorname{SE}(\bar g) \simeq\sqrt{\frac{\operatorname{Var}_p(g)} {N_\mathrm{eff}^{(g)}}} =\frac{1}{20}=0.05.$$

These formulas describe correlations within the sampled stationary distribution. If every chain remains in the same mode, its within-mode autocorrelations can be short and its apparent effective size large even though the global posterior is not represented. Dispersed starts, mode-sensitive summaries, and an independent exploration method are needed.

References

• LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration. “GWTC-5.0: Methods for Identifying and Characterizing Gravitational-Wave Transients.” arXiv:2605.27224 [gr-qc] (2026). Open PDF.
• Thrane, Eric, and Colm Talbot. “An Introduction to Bayesian Inference in Gravitational-Wave Astronomy: Parameter Estimation, Model Selection, and Hierarchical Models.” Publications of the Astronomical Society of Australia 36 (2019): e010. DOI. Open PDF, corrected arXiv v9.
• Veitch, John, et al. “Parameter Estimation for Compact Binaries with Ground-Based Gravitational-Wave Observations Using the LALInference Software Library.” Physical Review D 91 (2015): 042003. DOI. Open PDF.