Cosmological Likelihood Construction

A cosmological likelihood is a probability model for a versioned data product given a physical prediction and its nuisance parameters. It is not merely a $\chi^2$ expression copied from a paper. A trustworthy construction binds the data vector to its response, calibration, selection, covariance, foreground model, and preprocessing history; it also retains every normalization that changes across the comparison being made. This page develops that construction and uses the public Pantheon+ release as a concrete example.

Required background. Probability and Statistical Inference supplies conditional densities, nuisance marginalization, covariance, and coverage.

Helpful background. Einstein Equations supplies the gravitational prediction being tested; Relativistic Matter and Compact Objects supplies matter and source models; Relativistic Cosmology supplies light-cone observables and transfer functions; and Computational Relativity supplies numerical validation and convergence tests.

Start from the data-generating process

Let $\mathbf d$ be the released data vector. A useful schematic model is

$$\mathbf d =\mathsf R(\boldsymbol\alpha) \,\mathbf s(\boldsymbol\theta) +\mathsf T\boldsymbol\eta+\mathbf n,$$

where $\mathbf s(\boldsymbol\theta)$ is the cosmological signal, $\mathsf R$ contains the mask, window, beam, transfer, binning, or fiducial-coordinate response, $\boldsymbol\alpha$ are calibration parameters, $\mathsf T\boldsymbol\eta$ represents foreground or nuisance templates, and $\mathbf n$ is stochastic noise. The likelihood is then

$$\mathcal L(\boldsymbol\theta,\boldsymbol\alpha,\boldsymbol\eta) =p(\mathbf d\mid\boldsymbol\theta,\boldsymbol\alpha,\boldsymbol\eta, \text{selection},\text{release}).$$

Every word after the conditioning bar is part of the model. A map mask, catalog veto, redshift cut, multipole range, or object-quality threshold changes which random variable was observed. A calibration correction changes the response, and a fiducial-cosmology conversion can change both the mean prediction and its covariance. Hidden preprocessing cannot be recovered by adding a nuisance parameter later.

The Gaussian form and the terms one must not drop

For residual $\mathbf r=\mathbf d-\boldsymbol\mu(\boldsymbol\theta,\boldsymbol\eta)$ and positive-definite covariance $\mathsf C(\boldsymbol\theta,\boldsymbol\eta)$, the normalized Gaussian likelihood is

$$-2\ln\mathcal L =\mathbf r^{\mathsf T}\mathsf C^{-1}\mathbf r +\ln\det\mathsf C+n\ln(2\pi).$$

If $\mathsf C$ and the data dimension $n$ are fixed, the last two terms are constant with respect to the parameters. They may be dropped for posterior shape, but an absolute evidence requires the normalized likelihood; a common fixed factor cancels a Bayes factor only when both models genuinely use the same data and normalization. If $\mathsf C$ depends on cosmology or nuisance parameters, $\ln\det\mathsf C$ is no longer constant, and omitting it rewards models simply for inflating their errors. Likelihoods built from different data dimensions describe different observed random variables and should not be compared by silently retaining or dropping incompatible constants.

Numerically, factor $\mathsf C=\mathsf L\mathsf L^{\mathsf T}$ and solve triangular systems. Explicitly forming $\mathsf C^{-1}$ is slower and less stable. A failed Cholesky factorization is a diagnostic: it can reveal an inconsistent mask on one covariance axis, roundoff in an estimated covariance, duplicated modes, or an invalid systematic construction.

Nuisance parameters are either sampled jointly or marginalized under declared priors:

$$\mathcal L_{\rm marg}(\boldsymbol\theta) =\int d\boldsymbol\eta\, \mathcal L(\boldsymbol\theta,\boldsymbol\eta) \pi(\boldsymbol\eta).$$

A precomputed correction and a sampled nuisance describing the same uncertainty must not both be applied. Conversely, fixing a correction while omitting its uncertainty makes the posterior too narrow.

For linear additive nuisances with a Gaussian prior, $\boldsymbol\eta\sim\mathcal N(\boldsymbol\eta_0,\mathsf S)$, independent of $\mathbf n\sim\mathcal N(0,\mathsf C)$, this integral has a transparent exact answer:

$$\boldsymbol\mu_{\rm marg} =\mathsf R\mathbf s+\mathsf T\boldsymbol\eta_0, \qquad \mathsf C_{\rm marg}=\mathsf C+\mathsf T\mathsf S\mathsf T^{\mathsf T}.$$

The extra covariance is the propagated uncertainty of the same nuisance amplitude shared across its template. Adding this term and also sampling that amplitude would count its uncertainty twice. The identity assumes a linear response and the stated independence; nonlinear calibration generally requires the integral itself.

Response, covariance, and compression belong together

A released compression $\mathbf y=\mathsf W\mathbf d$ has mean $\mathsf W\boldsymbol\mu$ and covariance $\mathsf W\mathsf C\mathsf W^{\mathsf T}$. Using the compressed mean with an uncompressed covariance, or changing the scale cuts without rebuilding both, is inconsistent. Window matrices also mix physical scales: a nominal cut on the plotted coordinate is not necessarily a sharp cut on contributing modes.

Estimated covariance introduces another stochastic layer. Finite simulations can bias its inverse and propagate uncertainty into parameter widths. The release must identify the simulation ensemble, number of realizations, estimator, conditioning or shrinkage, and any correction applied to the precision matrix. Non-Gaussian tails, connected higher-point covariance, parameter-dependent sample variance, and uncertain responses should be tested rather than hidden inside a fixed Gaussian approximation.

Planck’s 2018 high-multipole likelihood is a particularly clear example of an inseparable package: frequency cross-spectra, masks, beams, polarization efficiencies, foreground templates, multipole cuts, nuisance priors, and covariance enter together. Planck Collaboration 2020, §§2–4 and appendices documents that construction and its validation; its cosmological parameters cannot be reproduced from a plotted spectrum alone.

Worked example using the Pantheon+ release

Pantheon+ provides 1,701 light-curve rows for 1,550 distinct supernovae, so repeated objects carry off-diagonal statistical covariance. Its cosmology release builds

$$\mathsf C_{\rm tot}=\mathsf C_{\rm stat}+\mathsf C_{\rm syst}, \qquad -2\ln\mathcal L=\Delta\mathbf D^{\mathsf T} \mathsf C_{\rm tot}^{-1}\Delta\mathbf D,$$

with systematic covariance obtained from calibrated response variations. Brout et al. 2022, §2.2, Eqs. (5)–(9), and Appendix C gives the construction and warns that plotting-only diagonal errors are not the likelihood covariance.

For the supernova-only Hubble diagram, use the release coordinate m_b_corr, the redshifts zHD and zHEL, and the selection cut zHD > 0.01. Reshape the full $1701\times1701$ covariance first, then apply the same ordered Boolean mask to both axes; this leaves 1,590 rows. The data-release paper records the quality cuts and repeated-light-curve structure in Scolnic et al. 2022, §2.2, Table 2, and Appendices A–B.

Supernovae without an absolute calibrator constrain relative distances. If $\mathbf r(\boldsymbol\theta)$ is the vector before subtracting a common intercept $M$, analytic profiling gives

$$\widehat M =\frac{\mathbf 1^{\mathsf T}\mathsf C^{-1}\mathbf r} {\mathbf 1^{\mathsf T}\mathsf C^{-1}\mathbf 1},$$ $$\chi^2_{\rm prof} =\mathbf r^{\mathsf T}\mathsf C^{-1}\mathbf r -\frac{(\mathbf 1^{\mathsf T}\mathsf C^{-1}\mathbf r)^2} {\mathbf 1^{\mathsf T}\mathsf C^{-1}\mathbf 1}.$$

Profiling chooses the best intercept for each cosmology; marginalization averages over its prior. To see the difference, write $A=\mathbf1^{\mathsf T}\mathsf C^{-1}\mathbf1$ and complete the square: $\chi^2(M)=\chi^2_{\rm prof}+A(M-\widehat M)^2$. For a uniform prior of finite width $\Delta M$, broad compared with $A^{-1/2}$ and with $\widehat M$ far from its endpoints,

$$\mathcal L_{\rm marg}(\boldsymbol\theta) \simeq\mathcal L_{\rm prof}(\boldsymbol\theta) \frac{\sqrt{2\pi/A}}{\Delta M}.$$

With fixed covariance this factor does not change posterior shape, but it does enter an absolute evidence. An improper flat prior has no normalized evidence. If $A$ changes with cosmology or a prior endpoint is approached, even the posterior shape differs from profiling.

Substituting a Cepheid-calibrated distance column would answer a different question and import an absolute calibration. The exact release repository is commit c447f0fea703fcd0fff57de5000947b5ca81286b; it has no license file, so this page gives a retrieval recipe rather than redistributing the large data objects.

Selection and blinding change the conditional density

If an object is observed only when event $A$ occurs, the correct density is

$$p(d\mid\boldsymbol\theta,A) =\frac{p(A\mid d,\boldsymbol\theta) p(d\mid\boldsymbol\theta)}{P(A\mid\boldsymbol\theta)}.$$

The denominator is the selection efficiency. Ignoring it generally biases population or distance parameters because detectable objects are not a random sample of the latent population. Completeness simulations, trigger thresholds, target selection, quality cuts, sky masks, and redshift failures therefore belong in the release record.

Blinding protects choices from result-driven tuning, but the record must say what was hidden, which diagnostics were examined, and whether any cut changed after unblinding. A post-unblinding change is not automatically wrong; it must be identified and stress-tested on simulations or independent data.

Validation and clean-room reproduction

Likelihood validation has three distinct layers. Unit tests check normalization, dimensions, ordering, and known limiting cases. Injection tests generate data through the complete forward model and verify bias and interval coverage across the parameter regime. Data checks reproduce a released likelihood value, best-fit residual, or posterior summary and run posterior predictive or goodness-of-fit statistics that were not tuned on the same discrepancy.

For a clean rerun, archive the data vector, covariance, response matrices, masks, nuisance templates, simulations, likelihood source, environment lockfile, configuration, and hashes. Record units, fiducial transformations, calibration and foreground conventions, all selection cuts, random seeds, optimizer or sampler settings, and convergence criteria. ACT DR6.02 demonstrates the scale of this contract: its public products and Louis et al. 2025, §§3–7 couple maps, passbands, beams, masks, transfer corrections, simulations, covariance, foregrounds, and approximately two thousand null tests.

Results and limitations through 2026-08-19

As of 2026-08-19, Pantheon+, Planck 2018, ACT DR6.02, and DESI DR2 offer versioned public likelihood or compressed-likelihood examples with different calibration domains and selections. None licenses the practice of mixing a vector from one release with covariance, masks, nuisance priors, or corrections from another. This page has specified the Pantheon+ algebra and retrieval boundary but has not claimed a clean-environment reproduction of its full posterior.

The most consequential failure modes are parameter-dependent covariance without its determinant, finite-simulation precision bias, non-Gaussian tails, inconsistent masks, hidden preprocessing, doubled corrections, hard nuisance boundaries, post-selection inference, look-elsewhere choices, and a normalization silently discarded before model comparison.

The durable scope is likelihood construction and validation. It does not provide undated best-fit values, replace instrument-specific release documentation, or certify that a collaboration pipeline has been reproduced merely because its public likelihood can be evaluated.

Exercise: the cost of ignoring shared noise

Two measurements of a common mean $\mu$ have unit variance and correlation $\rho=0.5$. Both report the value $1$. Find the generalized-least-squares estimate and its variance. Compare with the result obtained by multiplying two independent unit-variance likelihoods.

Solution

For

$$\mathsf C= \begin{pmatrix} 1&\rho\\ \rho&1 \end{pmatrix},$$

symmetry gives $\widehat\mu=1$. Its variance is $(\mathbf1^{\mathsf T}\mathsf C^{-1}\mathbf1)^{-1}=(1+\rho)/2=0.75$. Treating the data as independent sets $\rho=0$ and gives variance $1/2$. The central value happens not to move because the two observations agree, but the independent product overstates the information and reduces the standard deviation from $\sqrt{0.75}=0.866$ to $\sqrt{0.5}=0.707$. A shared calibration nuisance or explicit cross-covariance restores the correct uncertainty.

References

• Brout, Dillon, et al. 2022. “The Pantheon+ Analysis: Cosmological Constraints.” Astrophysical Journal 938: 110. arXiv:2202.04077v2. DOI. Open PDF.
• Louis, Thibaut, et al. (ACT Collaboration). 2025. “The Atacama Cosmology Telescope: DR6 Power Spectra, Likelihoods and ΛCDM Parameters.” arXiv:2503.14452v2. Open PDF. Official DR6.02 release.
• Planck Collaboration. 2020. “Planck 2018 Results. V. CMB Power Spectra and Likelihoods.” Astronomy & Astrophysics 641: A5. DOI. Open PDF.
• Scolnic, Dan, et al. 2022. “The Pantheon+ Analysis: The Full Data Set and Light-curve Release.” Astrophysical Journal 938: 113. arXiv:2112.03863v2. DOI. Open PDF.