Reproducible Cosmological Inference

Reproducing a cosmological posterior means recovering the same probability model and demonstrating numerically controlled sampling—not merely obtaining a similar contour. The record must bind data, likelihood, theory, priors, parameter basis, environment, configuration, random streams, chains, convergence diagnostics, and derived-quantity code. This page gives that contract and reports a transparent correlated-Gaussian benchmark with an analytic answer, deterministic quadrature, and seeded Metropolis chains. It does not claim to reproduce a collaboration-scale likelihood.

Required background. Probability and Statistical Inference supplies Gaussian conditioning, Markov chains, autocorrelation, and convergence diagnostics.

Helpful background. Einstein Equations supplies the field-equation model behind production predictions; Relativistic Matter and Compact Objects supplies source and nuisance physics; Relativistic Cosmology supplies observable calculations and parameter maps; and Computational Relativity supplies refinement, residual, and independent-benchmark discipline.

An executable record has four layers

Scientific identity records the model equations, parameter basis, fixed parameters, priors including Jacobians, neutrino and recombination prescriptions, nonlinear model, nuisance definitions, likelihood normalization, and derived maps. Data identity records release names, retrieval dates, ordered selections, masks, response matrices, covariance, calibration, units, and hashes. Software identity records source commits, package lockfiles, compilers or interpreters, platform, configuration, caches, and executable hashes. Run identity records seeds, initial states, proposal or sampler settings, stopping rules, chain filtering, weights, and all attempted runs.

The posterior artifact then includes raw chains, column metadata, sample weights, log prior and component log likelihoods, convergence history, and scripts that regenerate summaries. A plot is an output, not the record. A chain that ended without an error is not thereby converged.

Analytic benchmark: product of correlated Gaussians

The local benchmark uses a two-parameter likelihood and prior,

$$\boldsymbol\mu_L=(1.2,-0.7), \qquad \mathsf C_L= \begin{pmatrix} 1.0&0.6\\ 0.6&1.5 \end{pmatrix},$$ $$\boldsymbol\mu_P=(0,0), \qquad \mathsf C_P= \begin{pmatrix} 4.0&-0.4\\ -0.4&2.5 \end{pmatrix}.$$

Completing the square gives the independent reference solution

$$\mathsf C_*=\left(\mathsf C_L^{-1}+\mathsf C_P^{-1}\right)^{-1},$$ $$\boldsymbol\mu_*=\mathsf C_* \left(\mathsf C_L^{-1}\boldsymbol\mu_L +\mathsf C_P^{-1}\boldsymbol\mu_P\right).$$

In IEEE-754 binary64 arithmetic, the script obtains

$$\boldsymbol\mu_*=(1.0649298597,-0.5674348697),$$ $$\mathsf C_*= \begin{pmatrix} 0.7214428858&0.2729458918\\ 0.2729458918&0.8822645291 \end{pmatrix}.$$

This tests covariance algebra, correlation signs, prior inclusion, and derived moments without requiring an external scientific package.

Two independent numerical routes

The numerical algorithm has two independent implementations: midpoint-grid quadrature and random-walk Metropolis sampling, both compared against analytic completion of the square.

The deterministic route evaluates the unnormalized posterior on midpoint grids over $[-8,8]^2$ with $41^2$, $81^2$, and $161^2$ cells, normalizes the weights, and computes first and second moments. All three grids agree with the analytic moments at roughly binary64 rounding precision; the finest maximum mean and covariance errors are $5.77\times10^{-15}$ and $6.42\times10^{-14}$.

The stochastic route uses four random-walk Metropolis chains. The proposal covariance is the analytic posterior covariance multiplied by scale $1.5^2$. Initial points are displaced by $(\pm3,\pm3)$ in whitened coordinates. Each chain takes 60,000 steps and discards 10,000 as burn-in. The hexadecimal seeds are 0x13579bdf, 0x2468ace0, 0x1badb002, and 0xc001d00d. This uses the analytic result to design an efficient proposal, but the acceptance decision and sampled moments are computed independently.

The sampling rule is short enough to understand without treating a sampler as a black box. From the current point $\boldsymbol\theta$, draw $\boldsymbol\theta'=\boldsymbol\theta+\boldsymbol\epsilon$ with $\boldsymbol\epsilon\sim\mathcal N(0,s^2\mathsf C_*)$. Let $\ell(\boldsymbol\theta)$ be the sum of the Gaussian log likelihood and log prior, and draw $u$ uniformly on $(0,1)$. Accept the proposal when

$$\ln u\lt\min\{0,\ell(\boldsymbol\theta')-\ell(\boldsymbol\theta)\}.$$

Otherwise keep the old point as the next chain sample. Keeping rejected states is essential: discarding them would sample the jump process rather than the desired posterior. The proposal is symmetric, so its forward and reverse densities cancel in this rule; a state-dependent or asymmetric proposal requires their ratio. Thus proposal tuning changes efficiency, while the likelihood, prior, and acceptance rule determine the target.

The benchmark was rerun on 2026-08-28 with Node.js 24.18.0, darwin/arm64, and JavaScript Number binary64. It is dependency-free: standard library only, with no external scientific packages. The linked MIT-licensed executable has SHA-256 94ca19a2dc80cb7cede79bfc0e1f548e76f21d7e03e3da4ac2894732c4ab0743.

Convergence and error results

The acceptance fractions were $0.3997$, $0.3988$, $0.3982$, and $0.4019$. The table recomputes pooled moments and split-$\widehat R$ from increasing retained prefixes of each post-burn-in chain.

Kept samples per chain	Maximum mean error	Maximum covariance error	Split $\widehat R_x$	Split $\widehat R_y$
5,000	1.301208e-2	5.059313e-2	1.00167	1.00135
15,000	9.463134e-3	1.566573e-2	1.00023	1.00025
50,000	2.260865e-3	7.609589e-3	1.00011	1.00003

At 50,000 retained samples per chain, the initial-positive-pair effective sample estimates, truncated at lag 250, are $25{,}735.9$ for $x$ and $25{,}595.0$ for $y$. The posterior-mean Monte Carlo standard errors are $5.29\times10^{-3}$ and $5.90\times10^{-3}$. The observed mean errors are smaller than those one-standard-error scales. Split-$\widehat R$ near one, adequate effective sample size, stable moments, plausible acceptance, and agreement with two independent references jointly support the result; none alone is sufficient.

The same script reruns all four seeds while changing burn-in and proposal scale. Each comparison below uses 40,000 retained samples per chain; all runs take 60,000 total steps, so only the declared analysis window or proposal is changed.

Setting	Burn-in	Proposal scale	Mean acceptance	Maximum mean error	Maximum covariance error	Maximum split $\widehat R$
Earlier burn	5,000	1.5	0.3997	2.727478e-3	1.330674e-2	1.00014
Baseline	10,000	1.5	0.3997	3.547962e-3	6.151950e-3	1.00013
Later burn	15,000	1.5	0.3997	1.992142e-3	8.335121e-3	1.00024
Smaller proposal	10,000	0.8	0.6280	2.577168e-3	8.725141e-3	1.00015
Larger proposal	10,000	2.4	0.2324	3.211830e-3	1.165612e-2	1.00038

Acceptance changes in the expected direction when the proposal changes, but the posterior moments and split diagnostics remain stable. This comparison does not prove a general sampler is well mixed: the target is a unimodal Gaussian and the proposal orientation uses its analytic covariance. It does show that the reported agreement is not an artifact of one retained length, one burn-in cut, one proposal scale, or one random stream.

Exact rerun and expected output

Use the standalone download-and-run recipe:

curl -fL https://blackhole.org/code/cosmology/cosmology-benchmarks.mjs -o cosmology-benchmarks.mjs
node cosmology-benchmarks.mjs --case inference

The run should identify Node and platform, print the analytic moments, report the three grid results, list four seeds and acceptance fractions, print the retained-length and five-row burn-in/proposal tables, report effective sample sizes and Monte Carlo errors, and end with PASS: 1 benchmark case. A different runtime or executable hash is a new execution record, even if it passes.

Production-scale inference requires a stronger pin

Cobaya provides a modular graph of parameters, theory calculations, likelihoods, priors, and samplers, with caching and output metadata. Torrado and Lewis 2021, §§II–V and Appendix A describes the architecture, fast–slow parameter blocking, MCMC convergence monitoring, and reproducibility metadata. The official Cobaya 3.6.2 documentation describes input, components, installation, samplers, post-processing, and cosmological likelihood interfaces.

A production example for this chapter must pin Cobaya v3.6.2 at commit 899f30a; archive the exact documentation snapshot and its SHA-256; record Python, compiler, MPI, numerical libraries, theory code and commit, likelihood packages and data, YAML, environment lockfile, source or executable hashes, package path, random seeds, proposal covariance, stopping criteria, and raw chains; and verify the likelihood value and posterior against an independent implementation or released reference.

The official 208-page PDF Cobaya Documentation, Release 3.6.2, generated 2026-03-27, was retrieved on 2026-08-19 from Read the Docs and hashed as 5dc1a4206db0571871a93b3ec4af8b0f0197ce820bd9f7b4d0f5f0dbf095499d. Torrado and Lewis 2026, Cobaya Documentation 3.6.2, printed pp. 11–13, 25–43, and 79–80 specifies the Gaussian YAML, run and outputs; effective input, output metadata, and priors; and MCMC convergence with its necessary-not-sufficient boundary. This documents the software interface and convergence criteria, but the page does not report a production Cobaya run: the external theory and likelihood data, Python environment, configuration, chains, and independent result checks were not executed here.

Failure modes and diagnostic responses

Failure mode	Symptom	Required response
Unpinned data or hidden defaults	Same YAML changes after reinstall	Hash every input and expand effective configuration
Covariance normalization error	Biased moments or unstable evidence	Check analytic Gaussian and retained determinant terms
Prior truncation	Posterior mass or derived tails reach a boundary	Broaden justified support and rerun rather than reweight outward
Poor mixing	Stable acceptance but large autocorrelation or split disagreement	Reparameterize, improve proposal, extend independent chains
Seed dependence	Summaries shift beyond Monte Carlo error	Add seeds, inspect multimodality, and compare independent samplers
Stale cache	Theory output disagrees with clean execution	Hash cache keys, clear only documented derived caches, rerun cleanly
Basis mismatch	Chain labels and derived code use different units or conventions	Store parameter schema and test transformations with Jacobians
Omitted derived code	Base posterior reproduces but reported observables do not	Version and test every derived mapping alongside the chain

The benchmark deliberately avoids evidence estimation: tail-sensitive normalization needs additional algorithms and convergence tests. Trotta 2008, §§4.2–4.3 explains why evidence is a separate numerical integral rather than a by-product of an ordinary posterior chain.

Exercise: interpret convergence jointly

At 50,000 retained samples per chain, the maximum mean error is $2.26\times10^{-3}$, while the larger componentwise Monte Carlo standard error is $5.90\times10^{-3}$ and both split-$\widehat R$ values are below $1.001$. Does this prove that any cosmological chain with the same diagnostics is correct?

Solution

No. Here the observed error being below the Monte Carlo standard-error scale is consistent with sampling variation, and split-$\widehat R$ plus the effective sample estimates support stationarity and mixing for the monitored moments. Correctness is stronger: it follows here only because the same target has an analytic solution and an independent grid calculation. A production chain can converge precisely to the wrong posterior if its likelihood, covariance, prior, units, selection, or theory code is wrong; it also may miss an unmonitored tail or mode.

Scope boundary

This page establishes a reproducibility record and a bounded numerical benchmark. It does not publish undated cosmological best-fit values, reproduce Planck, ACT, DESI, or supernova likelihoods, certify Cobaya merely from a version string, or treat structural validation as independent scientific review.

References

• Cobaya Project. 2026. Cobaya Documentation. Release 3.6.2, generated 27 March 2026. Jesús Torrado, Antony Lewis, et al. 208 pp. Official PDF.
• Torrado, Jesús, and Antony Lewis. 2021. “Cobaya: Code for Bayesian Analysis of Hierarchical Physical Models.” Journal of Cosmology and Astroparticle Physics 05: 057. arXiv:2005.05290v2. DOI. Open PDF. Official v3.6.2 documentation.
• Trotta, Roberto. 2008. “Bayes in the Sky: Bayesian Inference and Model Selection in Cosmology.” Contemporary Physics 49: 71–104. DOI. Open PDF.