Graphical Profile Likelihoods

• Graphical profile likelihoods are a statistical method that visualizes the likelihood function by varying parameters of interest while optimizing out nuisance parameters.
• They enable direct derivation of confidence intervals through plots of likelihood-ratio statistics calibrated by chi-squared thresholds, bridging likelihood geometry with uncertainty quantification.
• Their applications span disciplines such as cosmology, geostatistics, and particle physics, where they diagnose asymmetries, boundary issues, and non-Gaussian behaviors in parameter estimates.

Searching arXiv for recent and foundational papers on graphical profile likelihoods and related computational methods. arXiv search query: "graphical profile likelihood cosmology profile likelihood" Graphical profile likelihoods are likelihood-based inferential displays in which one or more parameters of interest are varied while all remaining nuisance parameters are maximized out, and the resulting likelihood-ratio statistic is plotted as a curve, contour, or surface. In the notation used in contemporary cosmology, if the likelihood is $L(\theta,\psi)=P(\mathrm{data}\mid \theta,\psi)$, then the profile likelihood is $L_p(\theta)=\max_{\psi}L(\theta,\psi)$ and the graphical construction is based on $q(\theta)=-2\ln\lambda(\theta)$, where $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$; confidence sets are then read off from the region where $q(\theta)\le c_\alpha$ (Herold et al., 2024). In Gaussian settings, $q$ coincides with the familiar $\Delta\chi^2$, which makes profile plots a natural bridge between likelihood geometry, asymptotic testing, and visual uncertainty quantification (Herold et al., 2024).

1. Formal definition and inferential object

The central object is the profile likelihood obtained by partitioning the parameter vector into a parameter of interest and nuisance parameters. With scalar or vector parameter of interest $\theta$ and nuisance vector $\psi$, one defines

$L_p(\theta)=\max_{\psi}L(\theta,\psi),$

and the likelihood-ratio statistic

$L_p(\theta)=\max_{\psi}L(\theta,\psi)$0

For Gaussian problems, $L_p(\theta)=\max_{\psi}L(\theta,\psi)$1 equals $L_p(\theta)=\max_{\psi}L(\theta,\psi)$2, and in practice one plots $L_p(\theta)=\max_{\psi}L(\theta,\psi)$3 and reads off intervals from fixed “iso-$L_p(\theta)=\max_{\psi}L(\theta,\psi)$4” levels (Herold et al., 2024).

An equivalent formulation used in several application domains is based on the log-likelihood. In diffusion models, the normalized log-likelihood is written as

$L_p(\theta)=\max_{\psi}L(\theta,\psi)$5

and the normalized profile log-likelihood for a parameter block $L_p(\theta)=\max_{\psi}L(\theta,\psi)$6 is

$L_p(\theta)=\max_{\psi}L(\theta,\psi)$7

with the likelihood-ratio statistic

$L_p(\theta)=\max_{\psi}L(\theta,\psi)$8

used for graphical cutoff construction (Simpson et al., 2020).

In geostatistics, the same structure appears after concentrating the Gaussian likelihood over regression parameters or over both regression and variance parameters. For a scalar parameter of interest $L_p(\theta)=\max_{\psi}L(\theta,\psi)$9, one may define the relative profile likelihood

$q(\theta)=-2\ln\lambda(\theta)$0

and the deviance

$q(\theta)=-2\ln\lambda(\theta)$1

then read confidence intervals from the set where $q(\theta)=-2\ln\lambda(\theta)$2 remains below the relevant $q(\theta)=-2\ln\lambda(\theta)$3 threshold (Xu et al., 2023).

These formulations are mathematically equivalent expressions of the same inferential idea: graphical profile likelihoods visualize the best attainable fit as the parameter of interest is displaced away from its maximum-likelihood value. This gives a direct geometric view of curvature, asymmetry, flatness, and truncation that is often obscured by summaries based only on local curvature or by marginalization over nuisance parameters.

2. Graphical construction in one and several dimensions

The standard one-dimensional construction proceeds by choosing a grid of values $q(\theta)=-2\ln\lambda(\theta)$4, maximizing over nuisance parameters at each fixed $q(\theta)=-2\ln\lambda(\theta)$5, computing

$q(\theta)=-2\ln\lambda(\theta)$6

and plotting the resulting curve (Herold et al., 2024). The confidence interval is then the set $q(\theta)=-2\ln\lambda(\theta)$7, where $q(\theta)=-2\ln\lambda(\theta)$8 is determined by the desired confidence level and the assumed asymptotic distribution of $q(\theta)=-2\ln\lambda(\theta)$9 (Herold et al., 2024).

A closely related graphical convention plots normalized profile log-likelihoods rather than $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$0 directly. In the stochastic diffusion setting, univariate profiles are displayed as $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$1 against the parameter value, with a horizontal cutoff at $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$2 for approximate $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$3 univariate regions because $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$4 and $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$5 is equivalent to $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$6 (Simpson et al., 2020). Bivariate profiles are displayed as shaded regions where the profiled surface remains above $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$7, corresponding to $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$8 (Simpson et al., 2020).

In cosmology, graphical profile likelihoods are frequently displayed as one-dimensional $$\lambda(\theta)=L(\theta,\hat\psi_\theta)/L(\hat\theta,\hat\psi)$$9 curves or as two-dimensional confidence contours. The CONNECT framework was designed specifically to make “typical triangle plots normally associated with Bayesian marginalisation” feasible for profile likelihoods, by computing large numbers of one-dimensional and two-dimensional optimized profile points (Nygaard et al., 2023). Procoli follows the same visual convention, placing one-dimensional profile curves on diagonal panels and two-dimensional $q(\theta)\le c_\alpha$0 contours off-diagonal, and additionally decomposes the total profile into per-experiment contributions (Karwal et al., 2024).

For explicitly joint regions, the radial profile log-likelihood method provides a direct construction of the boundary. In two dimensions, one writes

$q(\theta)\le c_\alpha$1

and solves

$q(\theta)\le c_\alpha$2

for the smallest positive $q(\theta)\le c_\alpha$3 along each ray; the resulting points trace the confidence boundary (Jaeger, 2015). In three dimensions, the same idea extends using spherical coordinates and the unit direction

$q(\theta)\le c_\alpha$4

to construct a mesh for the joint confidence surface (Jaeger, 2015).

A different higher-dimensional construction is based on constrained optimization and ordinary differential equations. For a scalar objective $q(\theta)\le c_\alpha$5, the profile bound can be defined as maximizing or minimizing $q(\theta)\le c_\alpha$6 subject to the constraint that the log-likelihood remains above $q(\theta)\le c_\alpha$7, and the associated Karush–Kuhn–Tucker system yields ODEs that move along the likelihood contour itself (Deville, 2024). The same ODE-based technique applies to profile likelihood contours for couples of parameters (Deville, 2024).

3. Calibration, Wilks’ theorem, and coverage

The graphical method becomes a confidence procedure only after calibration of the cutoff $q(\theta)\le c_\alpha$8. Under Wilks’ theorem in the asymptotic regime, $q(\theta)\le c_\alpha$9 with $q$0 equal to the dimension of the parameter of interest (Herold et al., 2024). For $q$1, the commonly used thresholds are $q$2, $q$3, and $q$4; many cosmology applications also use $q$5 and $q$6 as $q$7 and $q$8 heuristics based on $q$9 (Herold et al., 2024). For $\Delta\chi^2$0, one uses $\Delta\chi^2$1 quantiles, with values such as $\Delta\chi^2$2 and $\Delta\chi^2$3 (Herold et al., 2024).

The same calibration appears across fields, with minor differences in reporting convention. Planck profile plots use horizontal lines at $\Delta\chi^2$4, $\Delta\chi^2$5, and $\Delta\chi^2$6 for one-dimensional intervals, and at $\Delta\chi^2$7, $\Delta\chi^2$8, and $\Delta\chi^2$9 for two-dimensional contours (Collaboration et al., 2013). CONNECT similarly adopts $\theta$0, $\theta$1, and $\theta$2 in one dimension, and $\theta$3, $\theta$4, and $\theta$5 in two dimensions, while noting that interval construction itself is not the focus of that work (Nygaard et al., 2023).

Setting	68% cutoff	95% cutoff
1 parameter under Wilks	$\theta$6	$\theta$7
2 parameters under Wilks	$\theta$8	$\theta$9
Gaussian near a boundary, $\psi$0	$\psi$1	$\psi$2

Correct coverage requires more than the formal $\psi$3 rule. The conditions emphasized in cosmology are large-sample asymptotics, regularity and identifiability, parameters in the interior of the allowed region, a unique MLE, differentiability, near-Gaussian behavior around the MLE, and absence of strong degeneracies that induce non-regular behavior when nuisance parameters are profiled out (Herold et al., 2024). When these conditions hold, the graphical intervals can have correct frequentist coverage (Herold et al., 2024).

Concrete validation results illustrate the point. For $\psi$4CDM with Planck 2018 Plik_lite data, mock tests showed $\psi$5 distributions consistent with $\psi$6 for each $\psi$7CDM parameter both with nuisance parameters fixed and with nuisance profiling, Wald’s relation held, and the alternative-hypothesis tests yielded non-central $\psi$8 as expected; the conclusion was that graphical profile-likelihood intervals had correct coverage for $\psi$9CDM under Plik_lite (Herold et al., 2024). In the Planck 2013 analysis, the base $L_p(\theta)=\max_{\psi}L(\theta,\psi),$0CDM profiles were described as “mostly parabolic,” and frequentist intervals agreed closely with Bayesian posterior constraints (Collaboration et al., 2013).

Coverage can also be better than Wald intervals in other model classes. In Gaussian geostatistical models with Matérn covariance, simulation studies and real-data applications showed that profile-based confidence intervals of covariance parameters and regression parameters had superior coverage to traditional standard Wald type confidence intervals (Xu et al., 2023). This is especially relevant when the global likelihood surface is flat or asymmetric in directions that a local Hessian cannot summarize adequately.

4. Boundaries, non-regularity, and failure modes

A major limitation of graphical profile likelihoods is that the Wilks calibration can fail in non-regular problems. The most common causes are physical boundaries, non-identifiability, discrete limits, non-Gaussian or multimodal likelihoods, very weak constraints, flat directions, and strong parameter degeneracies that make the likelihood asymmetric or not well approximated by a parabola (Herold et al., 2024). In geostatistics, analogous difficulties arise when the nugget is on the boundary $L_p(\theta)=\max_{\psi}L(\theta,\psi),$1 or when anisotropy approaches isotropy, producing flat or mixture $L_p(\theta)=\max_{\psi}L(\theta,\psi),$2 asymptotics (Xu et al., 2023).

Near a physical boundary, the asymptotic distribution of the likelihood-ratio statistic changes. For a single bounded parameter, the Chernoff boundary result gives the asymptotic distribution at the boundary as the mixture $L_p(\theta)=\max_{\psi}L(\theta,\psi),$3, implying more probability mass at $L_p(\theta)=\max_{\psi}L(\theta,\psi),$4 and lower effective cut thresholds than naive $L_p(\theta)=\max_{\psi}L(\theta,\psi),$5 (Herold et al., 2024). The boundary-aware profile-likelihood ratio in one dimension is written as

$L_p(\theta)=\max_{\psi}L(\theta,\psi),$6

and the modified Wald relation gives a linear branch when the unconstrained MLE lies in the unphysical region and a parabolic branch otherwise (Herold et al., 2024).

The Feldman–Cousins graphical method addresses this problem by using the likelihood-ratio statistic as the ordering rule in a Neyman construction, yielding unified intervals that avoid empty intervals and switch automatically between one-sided and two-sided behavior (Herold et al., 2024). Planck’s neutrino-mass analysis used the Feldman–Cousins prescription precisely because the unconstrained profile minimum for $L_p(\theta)=\max_{\psi}L(\theta,\psi),$7 with CMB-only data lay well within the unphysical negative-mass region; adding CMB lensing regularized the profile and the final robust frequentist upper limit was reported as $L_p(\theta)=\max_{\psi}L(\theta,\psi),$8 at $L_p(\theta)=\max_{\psi}L(\theta,\psi),$9 confidence for CMB+lensing+BAO (Collaboration et al., 2013).

The 2024 cosmology review refines this picture. In $L_p(\theta)=\max_{\psi}L(\theta,\psi)$00CDM with free sum of neutrino masses $L_p(\theta)=\max_{\psi}L(\theta,\psi)$01, the behavior was described as Gaussian near the physical boundary $L_p(\theta)=\max_{\psi}L(\theta,\psi)$02, many mocks had $L_p(\theta)=\max_{\psi}L(\theta,\psi)$03 when profiling over nuisance parameters, and the recommendation was to use boundary-corrected graphical intervals or Feldman–Cousins belts (Herold et al., 2024). By contrast, for $L_p(\theta)=\max_{\psi}L(\theta,\psi)$04CDM with free $L_p(\theta)=\max_{\psi}L(\theta,\psi)$05, the profile $L_p(\theta)=\max_{\psi}L(\theta,\psi)$06 was asymmetric and not well approximated by a parabola, mocks showed bimodality and an absence of $L_p(\theta)=\max_{\psi}L(\theta,\psi)$07 behavior, and the empirical $L_p(\theta)=\max_{\psi}L(\theta,\psi)$08 cut was $L_p(\theta)=\max_{\psi}L(\theta,\psi)$09 rather than $L_p(\theta)=\max_{\psi}L(\theta,\psi)$10; the graphical intervals were therefore flagged as approximate unless a full Neyman construction was carried out (Herold et al., 2024).

These cases correct a persistent misconception: a smooth-looking profile curve is not, by itself, evidence that standard $L_p(\theta)=\max_{\psi}L(\theta,\psi)$11 thresholds have correct coverage. The literature repeatedly emphasizes that near-boundary or strongly non-Gaussian profiles require explicit diagnostic tests, alternative cutoffs, or a full construction with guaranteed coverage (Herold et al., 2024).

5. Numerical realization and software ecosystems

The computational burden of graphical profile likelihoods comes from the need to solve many constrained or conditional optimization problems. In cosmology, early high-precision work interfaced the CLASS Boltzmann solver to the Planck likelihood and used Minuit2, specifically MIGRAD followed by HESSIAN, to build one-dimensional profiles by minimizing $L_p(\theta)=\max_{\psi}L(\theta,\psi)$12 over all remaining parameters at each grid point (Collaboration et al., 2013). Precision settings in CLASS were tuned so that numerical noise in $L_p(\theta)=\max_{\psi}L(\theta,\psi)$13 remained well below the $L_p(\theta)=\max_{\psi}L(\theta,\psi)$14 scales used for confidence intervals, and outliers were rejected with Minuit’s expected distance to minimum diagnostic (Collaboration et al., 2013).

The 2024 cosmology review released the code pinc (“profiles in cosmology”), a minimal profile-likelihood framework that computes profile-likelihood curves and intervals, supports boundary corrections, produces graphical outputs, and uses simulated-annealing minimization via MontePython’s step-size $L_p(\theta)=\max_{\psi}L(\theta,\psi)$15 and temperature $L_p(\theta)=\max_{\psi}L(\theta,\psi)$16 schedules while interfacing to CLASS (Herold et al., 2024). The same review recommends tuning grid density near the MLE, using conservative stopping criteria, re-running failed points, seeding multiple starts to avoid local minima, and using emulators such as CosmoPower, CosmicNet, and CONNECT when many mocks or full Neyman belts are required (Herold et al., 2024).

Two recent cosmology packages operationalize these recommendations in different ways. Procoli integrates with MontePython and CLASS, uses simulated annealing as a sequence of tempered Metropolis–Hastings optimizations, and constructs profiles sequentially by stepping the parameter of interest while re-optimizing nuisance parameters from the previous solution; it also records the component $L_p(\theta)=\max_{\psi}L(\theta,\psi)$17 contributions of individual experiments, allowing profile decomposition by dataset (Karwal et al., 2024). CONNECT replaces CLASS evaluations with a neural-network emulator and combines auto-differentiable likelihoods with a modified gradient-based basin-hopping scheme using parallel local BFGS solves; with default precision settings it achieves $L_p(\theta)=\max_{\psi}L(\theta,\psi)$18 relative to CLASS, a single likelihood-plus-gradient evaluation takes $L_p(\theta)=\max_{\psi}L(\theta,\psi)$19 once the TensorFlow graph is built, and the profile computation gains an additional speed-up of $L_p(\theta)=\max_{\psi}L(\theta,\psi)$20–$L_p(\theta)=\max_{\psi}L(\theta,\psi)$21 orders of magnitude relative to simulated annealing while maintaining excellent agreement (Nygaard et al., 2023).

Outside cosmology, numerical strategies differ with model structure. In the layered diffusion random-walk model, univariate and bivariate profiles were computed on fixed grids using MATLAB’s fmincon with bound constraints, warm starts from previous optimal nuisance values, and either an exact discrete-time Markov chain likelihood or a fast moment-matched Gamma approximation (Simpson et al., 2020). In Gaussian geostatistics, the computational bottleneck is repeated dense covariance factorizations; the proposed methodology uses internal reparameterizations, representative points obtained from quadratic approximations, and GPU-batched kernels for covariance construction, Cholesky factorization, triangular solves, and crossproducts (Xu et al., 2023).

High-dimensional multimodal models expose yet another computational issue: samplers configured for Bayesian posterior reconstruction may not resolve the sharp peaks relevant for profile likelihoods. In multi-dimensional SUSY scans, standard MultiNest settings with modest live points and loose tolerance were found to reconstruct posteriors accurately but to approximate the profile likelihood poorly; the tuned configuration required substantially larger n_live and much smaller tol, at much higher computational cost, to recover the global maximum and narrow spike regions reliably (Feroz et al., 2011).

6. Applications across disciplines

In cosmology, graphical profile likelihoods have become a standard complement to Bayesian credible intervals because they can diagnose prior sensitivity, volume effects, and boundary issues directly from the likelihood surface (Herold et al., 2024). For base $L_p(\theta)=\max_{\psi}L(\theta,\psi)$22CDM, both the 2013 Planck analysis and the 2024 review found close agreement between frequentist and Bayesian constraints, with nearly parabolic profiles and good coverage diagnostics (Collaboration et al., 2013). For neutrino mass, however, the profile often behaves differently because the unconstrained minimum may lie in the unphysical negative-mass region; boundary-aware treatment then becomes essential (Herold et al., 2024). In early dark energy, Procoli showed that the profile on $L_p(\theta)=\max_{\psi}L(\theta,\psi)$23 peaks near $L_p(\theta)=\max_{\psi}L(\theta,\psi)$24, while Bayesian posteriors peak near $L_p(\theta)=\max_{\psi}L(\theta,\psi)$25 because of prior-volume effects as $L_p(\theta)=\max_{\psi}L(\theta,\psi)$26 (Karwal et al., 2024).

In stochastic transport through heterogeneous media, graphical profile likelihoods are used primarily to assess practical identifiability rather than asymptotic coverage. Exact and approximate profiles for layer-specific hopping rates showed that identifiability depends strongly on release design and sample size, that parameters near the absorbing boundary are better identified, and that deeper layers often exhibit one-sided intervals hitting the upper bound (Simpson et al., 2020). When parameters were not practically identifiable, the same profile machinery supported a reduced-model strategy in which several layers were combined into a single effective layer, producing more peaked profiles and exit-time distributions visually indistinguishable from those of the full model under the observed design (Simpson et al., 2020).

In Gaussian geostatistics, graphical profiles are especially valuable for covariance parameters that are weakly identified under Matérn models. The paper’s key finding is that the Matérn shape parameter is often “quite flat but still identifiable,” meaning that the deviance has a unique minimum and can usually rule out very small values even when the confidence interval is wide (Xu et al., 2023). Real-data applications showed that profile-based intervals are often markedly wider than Wald intervals for covariance parameters, while two-dimensional contours in internal anisotropy coordinates are more regular than contours in the original range-angle parameterization (Xu et al., 2023).

In supersymmetric parameter inference, the application is not merely graphical convenience but discovery of high-likelihood structures of negligible posterior mass. The CMSSM study showed that the focus-point region could define the global best fit while lying largely outside Bayesian $L_p(\theta)=\max_{\psi}L(\theta,\psi)$27 credible regions, and that accurate profile maps required a sampler configuration fundamentally different from a posterior-oriented one (Feroz et al., 2011). This use case underscores that graphical profile likelihoods can reveal “spiky” regions that marginal posterior summaries systematically de-emphasize.

In extreme-value statistics, the emphasis shifts from static intervals to confidence bands and contours defined by the geometry of the log-likelihood surface. Constrained optimization and ODEs were used to derive profile-likelihood confidence limits for scalar functions such as GEV return levels and to trace profile contours for parameter pairs, with the ODE evolving tangentially along the likelihood contour itself (Deville, 2024). This suggests a broader role for graphical profile likelihoods as geometric objects, not merely as one-dimensional cut plots.

7. Relation to Bayesian inference and interpretive controversies

Graphical profile likelihoods and Bayesian posterior summaries answer different inferential questions. Frequentist intervals based on profiling depend only on the likelihood and do not incorporate prior-volume effects, whereas Bayesian credible intervals depend on priors, parameterization, and marginalization over nuisance dimensions (Herold et al., 2024). The literature therefore treats the two approaches as complementary rather than interchangeable (Nygaard et al., 2023).

Where the likelihood is close to Gaussian and interior regularity holds, the two methods often agree closely. This was the conclusion for Planck $L_p(\theta)=\max_{\psi}L(\theta,\psi)$28CDM analyses in both 2013 and 2024, where profile curves were nearly parabolic and frequentist intervals matched Bayesian results well (Collaboration et al., 2013). In such cases, agreement itself is informative because it indicates that neither prior specification nor marginalization volume is substantially distorting the main constraints.

Disagreement, however, is diagnostically important rather than pathological. In the massive-neutrino case, frequentist upper limits can be tighter than Bayesian limits for Planck-only and Plik_lite data because the extrapolated parabolic minimum lies in the unphysical negative-mass region; once BAO is added, the minimum shifts close to zero and the two approaches agree more closely (Herold et al., 2024). In early dark energy, the Bayesian posterior is driven toward vanishing EDE by prior-volume effects, whereas the profile tracks the likelihood preference directly (Karwal et al., 2024). In SUSY and other multimodal settings, profile maps can emphasize narrow high-likelihood spikes with tiny posterior mass (Feroz et al., 2011).

A second recurrent controversy concerns whether approximate profile curves extracted from generic samplers are adequate. The SUSY study showed that standard nested-sampling configurations suited to Bayesian evidence and posterior mapping were not sufficient for accurate profile-likelihood reconstruction, because they under-sampled narrow peaks and did not reliably refine the global maximum (Feroz et al., 2011). The practical implication is that a “profile” obtained by binning posterior samples is not automatically a reliable graphical profile likelihood. Procoli makes this point explicitly, warning against binned MCMC “profiles” because prior-volume effects and sampling noise can bias the curve and miss the true minimum (Karwal et al., 2024).

The broader methodological consensus is cautious. Graphical profile likelihoods are most informative when accompanied by diagnostics: curvature checks near the MLE, mock or Monte Carlo comparisons of the empirical $L_p(\theta)=\max_{\psi}L(\theta,\psi)$29 distribution to $L_p(\theta)=\max_{\psi}L(\theta,\psi)$30, Wald-relation checks, nuisance-insensitivity tests, and explicit handling of boundaries (Herold et al., 2024). When those diagnostics fail, the graph remains useful as a descriptive object, but its thresholded interval should be regarded as approximate unless supported by Feldman–Cousins or a full Neyman construction (Herold et al., 2024).