Prepared for submission to JCAP

CLiENT: A new tool for emulating

cosmological likelihoods using deep
arXiv:2512.17509v1 [[Link]] 19 Dec 2025

neural networks

Luca Janken,a Steen Hannestad, a

Thomas Tram,a Andreas
Nygaardb
a Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark
b Department of Astrophysics, University of Zurich, CH-8057 Zurich, Switzerland
E-mail: [Link]@[Link], steen@[Link], [Link]@[Link],
andreas@[Link], [Link]@[Link]

Abstract. Cosmological emulation of observables such as the Cosmic Microwave Back-

ground (CMB) spectra and matter power spectra have become increasingly common in recent
years because of the potential for saving computation time in connection with cosmological
parameter inference or model comparison. In this paper we present client (Cosmological
Likelihood Emulator using Neural networks with TensorFlow), a new method which circum-
vents the computation of observables in favour of directly emulating the likelihood function
for a data set given a model parameter vector. We find that the method is competitive with
observable emulators in terms of the required number of function evaluations, but has the dis-
tinct advantage of producing a surrogate likelihood which is completely auto-differentiable.
Using less than 2 × 104 function evaluations client typically achieves credible intervals
within better than 0.1σ of those obtained using the true likelihood and single-point emulator
precision better than ∆χ2 ∼ 0.5 across relevant regions in parameter space.
Contents

1 Introduction 1

2 Likelihood functions used 2

3 Emulator framework 3
3.1 Loss function 3
3.2 Activation function 4
3.3 Sampling methods for training data 5
3.4 Generating new proposals for training data 6

4 Metrics for testing the emulator precision 7

4.1 KL divergence and credible intervals 7
4.2 Profile likelihoods 8

5 Training convergence 9

6 Results 9
6.1 Tests on synthetic likelihoods 10
6.2 Tests on the Planck likelihood 11
6.3 Runtime considerations 12

7 Conclusion and outlook 13

A Quantiles of the χ2 -distribution. 17

1 Introduction

The traditional pipeline for performing cosmological Bayesian parameter inference or model
comparison involves the computation of one or more observables such as the Cosmic Mi-
crowave Background (CMB) spectrum or the large scale matter power spectrum for a very
large number of different model parameters - often requiring millions of model computations
(see e.g. [1, 2]). For CMB analysis a single point requires running an Einstein–Boltzmann
solver such as camb [3] or class [4] to compute the CMB power spectrum, followed by the
evaluation of the likelihood on this particular point. The CPU time requirement is typically
tens of core seconds for a single point, making parameter inference computationally expensive
and Bayesian model comparison prohibitively expensive. Analysing large scale structure in
the non-linear regime carries a much higher cost because it requires running an N -body sim-
ulation to compute the power spectrum, de facto making a brute force approach to Bayesian
model comparison impossible.
This problem has made emulation of observables increasingly popular in cosmology and a
number of different frameworks are now available for this purpose (see e.g. [5–14]). However,
since emulation of cosmological observables still requires a subsequent direct computation
of the likelihood, parameter inference or model comparison can still be very expensive in
terms of computational resources in the event that the likelihood is complex. Furthermore,
while the emulator possesses the desirable feature of being auto-differentiable, this property

–1–
is typically lost in the likelihood computation. This problem makes the direct emulation
of the likelihood function from the model parameter vector a desirable goal (note also that
direct emulation of data from a model makes likelihood-free inference possible. However, this
requires direct handling of data rather than a compressed version necessary).
Here, we present client, a general framework for emulating likelihood functions. The
current version allows for the emulation of any likelihood currently supported by the industry
standard Markov Chain Monte Carlo (MCMC) samplers MontePython [1] or Cobaya [2]
. The algorithm works by iteratively training the network on a cloud of points gathered so
that their density asymptotes to follow a tempered version of the actual likelihood function 1 .
This allows it to gather training data which maximises knowledge about the target likelihood
function while minimising the total computational cost. In Section 2 we outline the likelihood
functions used in the current work, and in Section 3 we provide details of the emulator
framework. Section 4 contains a discussion of the metrics we use for benchmarking the
emulator performance, and Section 5 describes the convergence criterion for the training
procedure. Section 6 contains numerical results for both synthetic and actual likelihoods,
and finally, in Section 7 we provide our conclusions.

2 Likelihood functions used

Our ultimate aim is to emulate the likelihood functions typically used in cosmological data
analysis, and the most complex likelihood function pipeline currently used is that of the CMB
anisotropies measured by the Planck mission [16]. Depending on the specific implementation,
this likelihood function depends on all cosmological model parameters (i.e. 6 in the standard
ΛCDM model) plus of order 20 experimental nuisance parameters. For the analysis presented
here we will use TT,TE,EE+lowE+lensing CMB data from Planck 2018 [16] in combination
with BAO measurements from BOSS DR12 [17], 6dF [18] and MGS [19], corresponding to
base2018TTTEEE_lensing_bao in MontePython. In the remainder of this paper we will
refer to this likelihood combination simply as the “Planck likelihood”.
Using this configuration, each single likelihood evaluation then consists of a call to class
taking approximately 10 core-seconds for the simplest models plus a call to the likelihood
which takes of order 200 ms on a single core. This typically means several hours or more of
wall time consumption for each test performed.
In order to speed up testing of the code we instead perform all our initial tests using
a synthetic and fully analytic likelihood function. The simplest choice for this function is a
multi-variate Gaussian of dimension N = Ncosmology + Nnuisance using the same parameter
covariance matrix, C, as the actual Planck likelihood.
1 1
 
L(θ) = q exp − (θ − µ)T C−1 (θ − µ) . (2.1)
(2π) det C
N 2

In order to also test a non-Gaussian yet simple analytic likelihood function we have
extended this setup with 2 additional parameters which are uncorrelated with all other pa-
rameters and for which the contribution to the likelihood is of the form e−θN +1 θN +2 An
example of such a likelihood could be the presence of an additional sterile neutrino of mass
ms and effective number Neff,s so that the total contribution to the energy density at late
1
Note that the method presented here bears some similarities to the tempered gathering of training data
described in [15] for the purpose of emulating gravitational wave signals.

–2–
times is ρ ∝ ms Neff,s . However, note that the two non-Gaussian parameter directions do
not correspond to the actual sterile neutrino model we use as one of the real data test cases
so that inferred parameter values for this case do not correspond to the sterile neutrino pa-
rameters estimated from Planck data. For the analytic likelihood we will refer to these two
non-Gaussian parameters simply as x28 and x29 . The similarity of this analytic function to
that of the actual Planck likelihood makes it very efficient for tuning hyperparameters of
both the network and the training procedure, and we will see later that the network and
training structure perform very well on actual Planck data. This means that in total we use
four different likelihoods to test the client framework:

• A 27-dimensional analytic Gaussian with covariance set to be identical to that of the

actual Planck data as described above.

• A 29-dimensional analytic likelihood with two new parameters added which are uncor-
related with all other parameters, but for which the contribution to the likelihood is
non-Gaussian.

• The 27-dimensional Planck likelihood for ΛCDM, corresponding to Ωb h2 , Ωcdm h2 , 100θs ,

ln(1010 As ), ns , τreio as cosmological parameters with an additional 21 nuisance param-
eters.

• The 29-dimensional Planck likelihood for ΛCDM + (Neff,s , ms )

3 Emulator framework

3.1 Loss function

Emulation of the likelihood function differs somewhat from general function emulation in the
sense that use of the emulated function typically comes in the form of Bayesian parameter
inference, maximum likelihood analysis, or Bayesian evidence calculations. All three methods
are sensitive to the exact value of the likelihood function in regions of large likelihood, but
insensitive in regions of moderate or low likelihood. This makes the use of a traditional mean
squared error loss non-optimal. Rather, it seems more appropriate to use a loss function
based on a relative error
!2 +
χ2surrogate − χ2exact
*
loss ∝ , (3.1)
χ2exact
where we have defined an effective χ2 in terms of the likelihood as

χ2 ≡ −2 (log(L) − log(Lbest )) , (3.2)

and log(Lbest ) is the largest value of L discovered so far. However, the loss-function in
equation (3.1) would diverge as χ2exact → χ2best = 0, so we regularise it in the following way,
!2 +
χ2surrogate − χ2exact
*
loss ∝ , (3.3)
χ2exact + ϵ

where ϵ > 0. For points where χ2exact ≪ ϵ, the loss reduces to the standard mean squared
error, whereas it becomes a relative error whenever χ2exact ≫ ϵ. Thus, ϵ decides how far away

–3–
 1000
σ = 23.0
3 σ = 3.0 500
Perfect prediction
2 0
true
2 − χ2

−500
χsurr

1
−1000
σ = 3.0
0 −1500
σ = 23.0
−2000
Perfect prediction
−80 −75 −70 −65 −60 −55 −50 −45 0 10000 20000 30000 40000
2
χtrue 2
χtrue

Figure 1: Surrogate error as a function of exact χ2 value for networks trained on the
27-dimensional analytic Gaussian likelihood. The left panel shows the region closest to the
best-fit while the right panel shows the error deep in the tails. A lower value of the parameter
κ in the loss function clearly reduce the scatter of the emulated likelihood around the best-fit
at the expense of increasing the scatter further away.

from the best-fit a point should be before it is penalised. This ϵ parameter will depend on
the dimensions of the problem, so instead we define a range based on the number of standard
deviations away from the best-fit. Using the well-known Wilson–Hilferty transformation [20]
we find the formula
!3
2 2
r
ϵ = ∆χ ≃ n 1 −
2
+κ , (3.4)
9n 9n
where κ is the number of standard deviations and n is the number of dimensions.
(Additional details can be found in appendix A). A priori, one would think that 1 ≲ κ ≲ 10
would be reasonable.
Figure 1 shows the surrogate error evaluated on a cloud of points with varying χ2 for
two different values of κ. As expected, a small value of κ pushes the network to emulate
points close to the likelihood maximum very well, but at the price of larger uncertainties in
the tails. It should also be noted that the network trained using a large value of κ exhibits
a systematic bias at small χ2 .
In the current work the main focus will be on emulating the likelihood well enough to
extract accurate posteriors, and we will therefore use κ = 3 for our actual runs. If the main
goal was profile likelihoods an even smaller value of κ might be favourable, whereas if the
goal was to calculate the Bayesian evidence it might be better to use a larger value in order
to be more accurate in the tails.

3.2 Activation function

The emulation problem discussed here consists of mapping a set of real model (and nuisance)
parameters to a real-valued function which has a minimum at some finite value, and which
diverges quadratically for values of model parameters far removed from the best-fit point.
Such a problem is typically not well handled by activation functions which are bounded, i.e.
tanh or similar, meaning that it is typically better to use e.g. a ReLU or similar. However,
as was also discussed in [6], ReLU is also not an ideal choice because it lacks expressivity
at negative values. Instead we have settled on the same choice of activation function as in
the connect framework, namely the Alsing et al. activation function [21] which replaces
ReLU with a linear function at negative values and is characterised by two hyperparameters:
The asymptotic slope, γ, of the linear function at negative values and the broadness, β, of

–4–
the transition region between the asymptotic functions at negative and positive values. The
vectorised function can be written as
  −1 
−β⊙x
f (x) = γ + 1 + e ⊙ (1 − γ) ⊙ x, (3.5)

such that β and γ can be trained for each node in each layer. This allows the network to find
optimal configurations of the activation function tailored for the output of each individual
node.
When using this activation function instead of ReLU or a sigmoid type function such
as tanh, it becomes abundantly clear that it leads to much better results given the same
amount of training (see [6] for a more detailed discussion).

3.3 Sampling methods for training data

It is well known that the precision of emulators is seriously affected by differences in avail-
able training data density across the model parameter domain. Assuming for simplicity that
the model parameter domain is given by a rectangular shape in Rn , the simplest selection
of training data points which has constant expected density is that of a latin hypercube.
Indeed this choice is often made when generating the initial training data set for an emula-
tor. However, the latin hypercube suffers greatly from the curse of dimensionality - in high
dimensions the number of training data points corresponding to regions of high likelihood
will be exceedingly small, and most of the emulating power of the network will be spent on
emulating regions of little actual interest.
Various methods for remedying this problem have been applied in the literature and
typically involve selecting training points with a density proportional to some function of the
posterior. While this method does lead to a close to optimal training point selection, it suffers
from the problem that the true posterior is unknown and requires running a prohibitively
large number of function evaluations to compute. However, assuming that a reasonable
approximation to the true posterior can be computed using only the surrogate likelihood
from a previously trained network, this can be used to sample new points in good locations
in the model parameter space. As an example, connect employs such a method to iteratively
train a network using new points selected from a MCMC run using the network trained on
the previous iteration of points.
The method we use here is conceptually similar to the one used in connect, but
employs a more sophisticated scheme for deciding whether a new point should be added.
At each iteration the aim is to acquire new training data such that the density, q(θ), of
training points is proportional to a tempered version of the true likelihood,

q(θ) ∝ L1/TT , (3.6)

with the “training temperature” TT a free and tunable hyper parameter. At this stage the
pre-existing training data consists of N points for which the true likelihood is known

S = {(θj , L(θj )), j ∈ {1, ..., N }}. (3.7)

The integral of q(θ) across the input parameter domain should fulfill the condition
Z
q(θ)dθ = N, (3.8)

–5–
and if we write q(θ) = cL1/TT we find that c = N/ L1/TT dθ. In order to relate this
R

continuous integral to properties of the existing training point cloud a number of different
approaches are possible. Here we will simply use a (very crude) approximation of estimating
the local density of points at the location of point p via the k th nearest neighbour
k
ρp ∼ ∼ rk−d , (3.9)
Vk
where d is the dimensionality of the input parameter space, Vk is the d-dimensional hypervol-
ume enclosing the k nearest neighbours, and rk is the distance to the k th nearest neighbour
and thus the radius of a hypersphere enclosing all k points. In practice, we find this crude
estimate to yield good results and to be independent of k for 5 ≤ k ≤ 50. We can now write
the integral as
Z X L(θp )1/TT
L1/TT dθ ∼ , (3.10)
p∈S
ρp

and substitute this into the equation for c

 −1
X L(θp )1/TT
c∼N  . (3.11)
p∈S
ρp

A given new point p∗ should only be added if the local estimated density is lower than the
target density, i.e., ρp∗ < q(θp∗ ) = cL(θp∗ )1/TT . Dividing both sides of this inequality with
the tempered likelihood function evaluated at the point yields the acceptance criterion for
the new point:
ρp∗
< c. (3.12)
L(θp∗ )1/TT

We note that the evidence integral can be written as

N N
L(θj )
Z
E= L(θj )Vj = (3.13)
X X
L(θ)dθ ∼ ,
j=1 j=1
ρj

so that for TT = 1, c ∼ N/E and the selection procedure will make each point contribute
roughly equally to the evidence integral. Conversely, for TT → ∞ we have L1/TT → 1 and
hence c ∼ N/ N j=1 ρj so that points are selected simply to be homogeneously spaced within
1
P

the input parameter volume.

3.4 Generating new proposals for training data

Having established a selection criterion for whether a given new point should be added to
the existing pool of training data we also need a method for generating a new set of points
from which to choose. In order to do this we employ a method essentially identical to what
is used in connect: Based on the existing data a network is trained and then used to
run a Markov Chain Monte Carlo process. The chain then directly samples the posterior of
L1/TMCMC , where TMCMC is the temperature of the chain and is a controllable hyperparameter.
This method for generating candidate points has the advantage of being cheap because it
only requires an evaluation of the surrogate likelihood and can therefore be run to high

–6–
convergence. Ensemble methods are well suited to this purpose because they can efficiently
utilize the surrogate likelihood’s parallelizability, which is why we employ the affine-invariant
ensemble sampler emcee [22]. The process works as follows:

1. Prior Restriction: Before sampling begins, the prior bounds are restricted using the
nstd parameter, which defines a hypercube of size ±nstd × σi centered on the fiducial
value for the ith parameter, where σi is a reference standard deviation. The actual
bounds used for sampling are determined by taking the intersection of this hypercube
with the original prior bounds—i.e., the more restrictive of the two is always chosen for
each parameter. This ensures that physically motivated prior constraints are always
respected while allowing tighter restrictions when desired.

2. Initialization: Starting positions for the ensemble of walkers are drawn uniformly
from these effective prior bounds.

3. MCMC Sampling: The sampler runs adaptively until convergence criteria are met:
a minimum number of effective samples (NESS × max(τ ), where NESS is the effective
sample size per walker and τ is the integrated autocorrelation time), stability in the au-
tocorrelation time (max ∆ττ < δτ,tol ), or until a maximum number of steps is reached.

4. Resampling: Once the MCMC run is complete, Ncand candidate points are sampled
1 1
−
from the chain using weights w = L TT TMCMC . The target value c is then computed
separately for the current training set S and the candidate set Scand , and we use the
larger of the two values to ensure adequate sampling in the region of highest posterior
density.

5. Candidate Evaluation: For each candidate point p∗ , we first evaluate the acceptance
criterion in eq. (3.12) using the surrogate likelihood Lsurr . Only candidates that pass
this initial test are then evaluated using the computationally expensive true likelihood
Ltrue , and the acceptance criterion is re-evaluated. This two-stage process significantly
reduces the number of expensive likelihood evaluations required.

6. Dynamic Update: After each point is accepted, the target concentration c is dynam-
ically updated to reflect the addition of the new training point. This ensures that the
acceptance criterion adapts as the training set grows, maintaining the desired sampling
density throughout the iterative process.

The typical configuration for this process uses nstd = 10, 216 walkers (i.e. 8 × Ndim ), NESS =
50, δτ,tol = 0.05, and Ncand = 10000.

4 Metrics for testing the emulator precision

Having developed the procedure for generating data and training the emulator network we
now proceed to test the precision of the trained emulator when used for statistical analysis
using a number of different metrics

4.1 KL divergence and credible intervals

A typical usage of likelihood emulators in the context of cosmological data analysis comes in
the form of Bayesian posterior distributions, typically presented as a set of 1D posteriors for

–7–
model and nuisance parameters obtained by integrating over all other parameters. Therefore,
an excellent benchmark for how well the emulator performs in typical tasks is to compare
the set of 1D posteriors obtained using the surrogate likelihood from a similar run using the
true likelihood.
The most often used benchmark for comparing two probability distributions, P and Q,
is the Kullback-Leibler divergence [23]
P (θ)
Z
DKL (P ||Q) = P (θ) log dθ, (4.1)
Q(θ)
where integration is performed over all parameters. The numerical complexity involved in
this computation is similar to the computation of the Bayesian evidence and for typical runs
can involve the computation of several million individual likelihood values. While this is
not overly expensive for the surrogate likelihood it is prohibitively expensive for the true
likelihood. Furthermore, since typical parameter estimation usage involves computing 1D
posterior distributions we instead use the set of 1D marginalised KL divergences, defined as
pi (θi )
Z
i
DKL (pi ||qi ) = pi (θi ) log dθi , (4.2)
qi (θi )
where pi (θi ) and qi (θi ) are the 1D marginalised posterior densities estimated from MCMC
samples using kernel density estimation (KDE). This heuristic benchmark is comparatively
cheap to compute and furthermore better suited to typical emulator usage.
The set of 1D marginalised KL divergences are combined into the rms KL divergence
!1/2
1 X i
DKL,rms = DKL (pi ||qi )2 , (4.3)
N i
which we will use as a simple benchmark for the precision of the emulator.
Since the computation of the marginalised KL divergence requires the computation of
all 1D marginalised posteriors we automatically also recover 1D credible intervals for all
parameters. These can also be used as a simple benchmark for following the precision of
the emulator during the training iterations. In order to arrive at a simple scalar expression,
we define the following “credible metric” as a benchmark for the precision with which the
emulator can recover credible intervals
|θi+ − θi,0
+ −
| + |θi− − θi,0 |
∆CM
i ≡ + − . (4.4)
θ0,i − θ0,i
±
Here, θi± and θi,0 are the upper and lower credible interval values for parameter i for the
emulator and target likelihoods respectively. In practise we will typically use maxi (∆CM ),
i.e. the parameter with the poorest quality of emulation as the actual benchmark for credible
intervals. We will refer to ∆CM for the 68% and 95% intervals as ∆CM68 and ∆CM95 respec-
|θi± −θi,0
±
|
tively. Finally, note that in the case of one-sided credible intervals we use ∆CM
i ≡ ±
|θ0,i |
instead. ± then refers to the one credible interval which is defined.

4.2 Profile likelihoods

Another quantity often used in analysis of cosmological data is the 1D profile likelihood,
which for parameter θj is defined as
Lprof (θj ) = max L(θj ), (4.5)
θi ̸=θj

–8–
i.e. for a given value of parameter θj the profile likelihood is obtained by maximising over
all other parameters. Unlike the marginalised posterior, the emulated profile likelihood is
typically very sensitive to the emulated likelihood function at very high values close to the
global best-fit (see e.g. [24]). Here, we will simply use the precision of the emulator at the
maximum likelihood point in parameter space as a useful and simple metric.

5 Training convergence

Once the iterative training algorithm has been set up the next question is how to stop the
iterations at an appropriate time so that sufficient accuracy has been reached, but without
using excessive computational resources. The connect framework uses the Gelman-Rubin
statistic, R − 1, calculated between chains from each iteration to estimate the difference
between iterations and terminates the procedure once R−1 crosses below a certain threshold.
We employ a similar approach, taking advantage of the fact that at each iteration
the emcee sampler generates an ensemble of walkers that can be combined into a single
flattened chain with well-defined statistical properties [22]. For each iteration i, we discard
the first 50% of steps to avoid burn-in effects, then flatten the remaining chain and compute
the mean vector µi and covariance matrix Ci from a random subsample of min(105 , Nchain )
points, where Nchain is the size of the remaining flattened chain.
The multivariate Gelman-Rubin statistic is then computed by comparing the statistics
from consecutive iterations. Following the approach of [25], we compute the between-chain
covariance
B = Cov(µi−1 , µi ), (5.1)
and the within-chain covariance
Ci−1 + Ci
W= . (5.2)
2
The convergence diagnostic R − 1 is then given by the largest eigenvalue of W−1/2 BW−1/2
(appropriately normalized). Training is terminated when R − 1 falls below a specified thresh-
old, typically set to 0.01.

6 Results

Having finalised the set-up of the training pipeline we subsequently tested client on the
four separate likelihoods described in Section 2. In each case the full training pipeline is run
using 10 iterations (rather than the R − 1 stopping criterion) in order to better visualise how
the emulator performance improved with increasing iterations. Results from emcee runs for
the emulator are compared to either the known analytic function for the two first likelihoods
and to standard MCMC runs performed with class and MontePython.
All runs presented in this work were performed with essentially identical hyperparameter
settings. The neural network architecture consisted of 5 hidden layers with 512 neurons each,
using the Alsing activation function with trainable β and γ parameters as described in Section
3.2. Training was performed for up to 5000 epochs with a batch size of 128, using the MSRE
loss function with κ = 3 (see Section 3.1), a learning rate of 10−4 , a random 10% validation
split, and early stopping with patience of 250 epochs. The initial training set consisted of
5000 points drawn from a latin hypercube restricted to ±10σ around the fiducial values. For
the iterative training procedure we used TT = 7 for the training temperature, k = 20 for
the k-nearest neighbour density estimation, and 5000 candidate points per iteration. The

–9–
 RMS KL (all)
10 0 R−1
Metric Value

10 −1

10 −2

∆ CM68
10 0 ∆ CM95
max(∆ CM )

10 −1

10 2 Surrogate at Reference
Linear region
best

10 Perfect agreement
2 − χ2

1
χsurr

0
-1
0 1 2 3 4 5 6 7 8 9
Iteration

Figure 2: Convergence metrics for the synthetic Gaussian likelihood as a function of it-
erations. Top panel: RMS of the KL divergence as defined in equation (4.3) and the Gel-
man–Rubin statistic of samples between iterations. Middle panel: Maximum of the 68% and
95% Credible Metrics defined in equation (4.4). Bottom panel: ∆χ2 as evaluated by the
network at the maximum likelihood point. Note that the plot is linear from −1 < ∆χ2 < 1
and logarithmic otherwise.

MCMC sampling was performed using emcee with 8 × Ndim walkers (i.e., 216 walkers for
the Gaussian and ΛCDM cases, and 232 walkers for the non-Gaussian and ΛCDM+(Neff,s ,
ms ) cases), a maximum of 105 steps per walker, and a burn-in of 5000 steps per walker.
Sampling used TMCMC = 7 and convergence criteria of NESS = 50 effective samples per
walker and δτ,tol = 0.05 for the relative change in autocorrelation time. Additional technical
parameters include a rebuild frequency of 50 candidates for the k-d tree, updating the k-NN
radii used for density estimation, chunk size of 5000 steps for autocorrelation calculations,
and autocorrelation thinning factor of 10.

6.1 Tests on synthetic likelihoods

Figure 2 shows the emulator performance as a function of iterations in the training pipeline.
From the figure it is clear that already after 2 iterations the network is able to recover posterior

– 10 –
credible intervals at a precision of around 0.05σ and achieves a single-point precision of better
than ∆χ2 ∼ 1. Running the training for 10 iterations pushes the precision somewhat further
up, but it should be noted that in later iterations comparatively fewer points are added
to the training data. We also note that R − 1 drops below 0.01 already after 2 iterations
indicating that training data already samples the relevant parameter space quite well at this
stage. Finally, from the figure it can be seen that 95% credible intervals are in general more
robustly recovered than 68% intervals because the likelihood is steeper further away from the
best-fit point. This makes it easier to localise cuts in parameter space.
Next, figure 3 shows a standard triangle plot for a subset of parameters in the 27-
dimensional Gaussian (named so that they match the corresponding cosmological parame-
ters). As can be seen from both figures, client is able to emulate the target likelihood to
a precision exceeding what is required for standard MCMC runs for parameter inference.
Typically, a total of around 2 × 104 are used in total to train the network in iteration 10,
with more than 90% of those points having been added already by iteration 2.
For the 29-dimensional synthetic likelihood with two non-Gaussian directions we find
essentially the same performance. However, as can be seen from figure 4 the training takes
somewhat longer to fully converge the non-Gaussian parameter directions . This happens
because the required pool of training points takes longer to populate along the increasingly
narrow funnel directions for the non-Gaussian directions. The ∆CM benchmarks are in this
case also dominated by the two non-Gaussian directions. The final result of the training is,
however, on par with the fully Gaussian likelihood. This can also be seen from the triangle
plot in figure 5 in which we see excellent agreement also for the two non-Gaussian directions.

6.2 Tests on the Planck likelihood

Next, we proceeded to use client on actual Planck data. The hyperparameter settings for
these runs were identical to those in the previous analytic likelihood tests. From figures 6
and 7 we see that client performs almost identically well on the real Planck data with the
ΛCDM model as it did with the 27-dimensional analytic likelihood. Convergence is achieved
quite quickly and in terms of our metrics, the emulator again recovers posterior intervals at
better than 0.1σ after just 2 iterations. Single-point accuracy of the emulator is at around
∆χ2 ∼ 2 after two iterations, but improves to around 0.5 in later iterations - again roughly
comparable to the analytic case.
Our final case, Planck data with the ΛCDM + (Neff,s , ms ) model is shown in figures 8
and 9. client again achieves excellent emulator precision in this case, although as was also
the case for the 29-dimensional analytic likelihood, convergence is slightly slower than in the
more Gaussian ΛCDM case. The credible intervals do exhibit a somewhat larger deviation
than in the previous cases, with the 68% intervals only recovered at around 0.2σ precision.
However, we note here that this deviation is most likely not entirely due to lack of precision
of the emulator, but rather a lack of convergence of the class-based run with which we
compare. The single-point emulator precision is also somewhat worse in this case than in
ΛCDM, but even here does not exceed ∆χ2 ∼ 1. Finally, we note that we have run with the
same hyperparameter settings for all four cases. In the ΛCDM + (Neff,s , ms ) this leads to
slightly fewer accepted training points than for our ΛCDM case due to the increased difficulty
of exploring parameter space. With more aggressive setting a larger training sample can be
collected and better emulator precision achieved.

– 11 –
 Surrogate (emcee)
Reference (cobaya)
Training Data

0.125

0.120
Ω cdm h 2

0.115

1.043
1.042
100 θs

1.041
1.040
1.039
3.10

3.05
ln(10 10 As )

3.00

2.95

0.98

0.96
ns

0.94

0.08

0.06
τreio

0.04

0.022 0.023 0.024 0.115 0.120 1.040 1.042 2.95 3.00 3.05 3.10 0.93 0.96 0.99 0.04 0.06 0.08
Ωb h 2 Ω cdm h 2 100 θs ln(10 10 As ) ns τreio

Figure 3: Triangle plot for the synthetic 27-dimensional Gaussian likelihood after the full
10 iterations. 19334 evaluations of the true likelihood with 14076 accepted points beyond the
initial latin hypercube of 5000 points. Only the first 6 parameters are shown in the figure.
The points show the final point cloud on which the iteration 10 network is trained.

6.3 Runtime considerations

The iterative pipeline described in section 3 proceeds in three phases:

1. Evaluate likelihood on new points.

2. Train the network.

3. Sample from network and draw new points.

– 12 –
 RMS KL (all)
10 0 RMS KL (non-Gaussian)
Metric Value

R−1
10 −1
10 −2
10 −3

∆ CM68
10 0 ∆ CM95
max(∆ CM )

10 −1

10 −2

10 2 Surrogate at Reference
Linear region
best

10 Perfect agreement
2 − χ2

1
χsurr

0
-1
0 1 2 3 4 5 6 7 8 9
Iteration

Figure 4: Convergence metrics for the synthetic likelihood that has two non-Gaussian pa-
rameters as a function of iterations. Top panel: RMS of the KL divergence as defined in
equation (4.3) and the Gelman–Rubin statistic of samples between iterations. Middle panel:
Maximum of the 68% and 95% Credible Metrics defined in equation (4.4). Bottom panel:
∆χ2 as evaluated by the network at the maximum likelihood point. Note that the plot is
linear from −1 < ∆χ2 < 1 and logarithmic otherwise.

We performed our runs on a node with two 32-core AMD EPYC 9334 CPU cores and
two Nvidia L40S-48GB GPUs. We used 62 CPU cores and 1 of the GPUs for our runs. In
phase 1 we used only CPU cores while the training in phase 2 was done on the GPU. In phase
3 we used the GPU for evaluating the network and the CPUs for the rest. With this setup,
each phase of the tree phases took a comparable amount of wall clock time. Specifically, for
the ΛCDM-run with the full Planck likelihood, the runtime was 100 minutes. The runtimes
for the individual phases are shown in table 1.

7 Conclusion and outlook

We have presented a novel framework for emulating cosmological likelihood functions directly
rather than through the normal two-stage approach of an observable emulation followed by a

– 13 –
 Surrogate (emcee)
Reference (cobaya)
Training Data

0.125

0.120
Ω cdm h 2

0.115

1.043
1.042
100 θs

1.041
1.040
1.039
3.10

3.05
ln(10 10 As )

3.00

2.95

4
3
x28

2
1

4
3
x29

2
1

0.022 0.023 0.024 0.115 0.120 1.040 1.042 2.95 3.00 3.05 3.10 1 2 3 4 1 2 3 4
Ωb h 2 Ω cdm h 2 100 θs ln(10 10 As ) x28 x29

Figure 5: Triangle plot for the synthetic 29-dimensional non-Gaussian likelihood after the
full 10 iterations. 17094 evaluations of the true likelihood with 11806 accepted points beyond
the initial latin hypercube of 5000 points. Only parameters 1, 2, 3, 4, 28, and 29 are shown
in the figure. The points show the final point cloud on which the iteration 10 network is
trained.

likelihood evaluation. The method was demonstrated to achieve excellent emulator precision
on the 29D case of Planck CMB data for the sterile neutrino extension of the ΛCDM model,
chosen specifically because the likelihood contains highly non-Gaussian directions in the
cosmological parameter space. For the specific setup presented here a total of less than
2 × 104 function evaluations (i.e. a class call followed by a likelihood evaluation) call is
required to achieve this precision, comparable to what is required for observable emulation
frameworks such as connect. With these settings client achieved a typical precision of

– 14 –
 10 1 RMS KL (all)
R−1
Metric Value

10 0
10 −1
10 −2

∆ CM68
∆ CM95
10 0
max(∆ CM )

10 −1

10 2 Surrogate at Reference
Linear region
best

10 Perfect agreement
2 − χ2

1
χsurr

0
-1
0 1 2 3 4 5 6 7 8 9
Iteration

Figure 6: Convergence metrics for the Planck data with the ΛCDM model as a function
of iterations. Top panel: RMS of the KL divergence as defined in equation (4.3) and the
Gelman–Rubin statistic of samples between iterations. Middle panel: Maximum of the 68%
and 95% Credible Metrics defined in equation (4.4). Bottom panel: ∆χ2 as evaluated by the
network at the maximum likelihood point. Note that the plot is linear from −1 < ∆χ2 < 1
and logarithmic otherwise.

Phase Total runtime

1: Evaluating likelihood on new points 37m (22,728 points)
2: Training the network 45m (1724 epochs on average)
3: Sampling from the network 18m

Table 1: Example of typical runtimes for a ΛCDM-run.

∆χ2 ∼ 0.3 − 0.5 in regions close to the best-fit point, and posterior credible intervals within
approximately 0.1σ of the true target credible intervals.
While client was tested only on Planck data as a worked example, it works on any
cosmological likelihood currently defined within the context of MontePython or Cobaya

– 15 –
 Surrogate (emcee)
Reference (montepython)
Training Data

0.13
ω cdm

0.12

0.11

1.043
1.042
100θ s

1.041
1.040
1.039

3.15
3.10
ln10 10 As

3.05
3.00
2.95

1.00

0.98
ns

0.96

0.94

0.10
τ reio

0.06

0.02
2.2 2.3 0.11 0.12 0.13 1.040 1.042 3.0 3.1 0.95 1.00 0.02 0.06 0.10
10 −2 ω b ω cdm 100θ s ln10 10 As ns τ reio

Figure 7: Triangle plot for the ΛCDM likelihood after the full 10 iterations. 22728 evalua-
tions of the true likelihood with 12654 accepted points beyond the initial latin hypercube of
5000 points. Only the 6 cosmological parameters are shown in the figure. The points show
the final point cloud on which the iteration 10 network is trained.

without modification. It is also worth pointing out that the client framework itself is
not restricted to use in cosmology, but will work on any likelihood function defined on a
continuous subset of Rn .

Reproducibility

We have used the publicly available client framework available at [Link]

AarhusCosmology/client_public to perform all runs presented in this work.

– 16 –
 10 1 RMS KL (all)
RMS KL (non-Gaussian)
Metric Value

10 0 R−1
10 −1
10 −2

10 1
∆ CM68
∆ CM95
max(∆ CM )

10 0

10 −1

10 2 Surrogate at Reference
Linear region
best

10 Perfect agreement
2 − χ2

1
χsurr

0
-1
0 1 2 3 4 5 6 7 8 9
Iteration

Figure 8: Convergence metrics for the Planck likelihood with the ΛCDM + (Neff,s , ms ).
Top panel: RMS of the KL divergence as defined in equation (4.3) and the Gelman–Rubin
statistic of samples between iterations. Middle panel: Maximum of the 68% and 95% Credible
Metrics defined in equation (4.4). Bottom panel: ∆χ2 as evaluated by the network at the
maximum likelihood point. Note that the plot is linear from −1 < ∆χ2 < 1 and logarithmic
otherwise .

Acknowledgements

We acknowledge computing resources from the Centre for Scientific Computing Aarhus
(CSCAA). A.N. was supported by the Carlsberg Foundation, grant CF24-1944

A Quantiles of the χ2 -distribution.

To relate our notion of "κ sigmas" to a concrete value of χ2 we start by computing the
probability integral inside κσ from the mean in a normal distribution. It is given by the
cumulative distribution function (CDF) Φ(x) of the normal distribution as

– 17 –
 Surrogate (emcee)
Reference (montepython)
Training Data

0.13
ω cdm

0.12

0.11

1.043
1.042
100θ s

1.041
1.040
1.039

3.15
3.10
ln10 10 As

3.05
3.00
2.95

0.6
Neff, s

0.4

0.2

2.5
2.0
1.5
ms

1.0
0.5

2.2 2.3 0.11 0.12 0.13 1.040 1.042 3.0 3.1 0.2 0.4 0.6 1 2
10 −2 ω b ω cdm 100θ s ln10 10 As Neff, s ms

Figure 9: Triangle plot for the ΛCDM + (Neff,s , ms ) likelihood after the full 10 iterations.
22024 evaluations of the true likelihood with 11467 accepted points beyond the initial latin
hypercube of 5000 points. Only 6 of the cosmological parameters are shown in the figure.
The points show the final point cloud on which the iteration 10 network is trained.

p = Φ(µ + κσ) − Φ(µ − κσ) , (A.1)

1 µ + κσ − µ 1 µ − κσ − µ
     
= 1 + erf √ − 1 + erf √ , (A.2)
2 σ 2 2 σ 2
κ
 
= erf √ . (A.3)
2

– 18 –
 We can now use the CDF of a χ2 -distribution with n degrees of freedom,
1 n x n x
   
CDFχ2 = n γ , ≡P , (A.4)
Γ 2 2 2 2


2

where γ(a, z) is the lower incomplete gamma-function and P (a, z) is the regularised lower
incomplete gamma-function. We now wish to find x such that
n x κ
   
P , = p = erf √ , (A.5)
2 2 2

and the solution to that equation is usually written as P −1 (p). Specifically, we have

n κ
  
χ2 (κ) = 2P −1 , erf √ . (A.6)
2 2
Evaluating this numerically is not very stable for large values of κ, so we will use a simple
1
approximation. If we define X ≡ χ2 , then (X/n) 3 is approximately normally distributed with
mean µ = 1 − 9n
2
and variance σ 2 = 9n
2
. This transformation is known as the Wilson–Hilferty
transformation [20]. This means that we have the approximate value
 1
x 3
≃ µ + κσ , (A.7)
n
if we are κσ away from the best-fit. Solving for x gives

∆χ2 = x ≃ n(µ + κσ)3 , (A.8)

!3
2 2
r
≃n 1− +κ . (A.9)
9n 9n

Note, however, that because the χ2 -distribution is one-sided, equation (A.9) is strictly
speaking finding the κ that solves Φ(µ + σ) − Φ(−∞) = p2 . Thus, we should in principle be
using a modified κ∗ to correct for this, where

∗
√ −1

κ
  
κ ≡ 2erf 2erf √ −1 . (A.10)
2
In the asymptotic limit one finds κ∗ ∼ κ2 − 2 log 2 so we have κ∗ → κ for κ ≫ 1.
p

In figure 10 we have shown the exact computation of ∆χ2 compared to the Wilson–Hilferty
approximation in equation (A.9), with and without the correction from κ∗ . When κ is very
large the approximation becomes less good since we are sampling deep in the tail of the
distribution where the normality of the cubically transformed variable no longer holds. Even
though the simple equation (A.9) breaks down in both ends, it is sufficient for our purpose
because we only care about the interpretation of κ in the range 1 ≲ κ ≲ 10 anyway.
Lastly, evaluating equation (A.6) numerically is non-trivial deep in the tails. Using e.g.
SciPy we may evaluate it as

2*[Link](n/2, [Link](kappa/[Link](2)))

– 19 –
 Exact (n=6) Exact (n=29)
104 ∆χ2 (κ), Wilson–Hilferty ∆χ2 (κ), Wilson–Hilferty
∆χ2 (κ∗ ), Wilson–Hilferty ∆χ2 (κ∗ ), Wilson–Hilferty
3
10
∆χ2

102

101

10−1 100 101 102

κ

Figure 10: ∆χ2 at x = κσ for n = 6 and n = 29. The exact calculation compared to the
Wilson-Hilferty approximation in equation (A.9), and the same equation where we have used

but for κ ≳ 38 the complementary error function underflows in double precision. For
κ > 38 we instead asymptotically expand the defining equation (A.5) using [Link]
[Link]/8.11.E2, which leads to the equation

x n−2
 
x = (n − 2) log +2+ A, (A.11)
2 x
2
!
n πκ2
 
A ≡ −2 log Γ + log + κ2 + 2 , (A.12)
2 2 κ
for x ≡ ∆χ2 . Equation (A.11) can be solved using fixed-point iteration and converges
rapidly for large κ.

References
[1] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, “Conservative Constraints on Early
Cosmology: an illustration of the Monte Python cosmological parameter inference code,”
JCAP 02 (2013) 001, arXiv:1210.7183 [[Link]].
[2] J. Torrado and A. Lewis, “Cobaya: Code for Bayesian Analysis of hierarchical physical
models,” JCAP 05 (2021) 057, arXiv:2005.05290 [[Link]].
[3] A. Lewis, A. Challinor, and A. Lasenby, “Efficient computation of CMB anisotropies in closed
FRW models,” Astrophys. J. 538 (2000) 473–476, arXiv:astro-ph/9911177.
[4] D. Blas, J. Lesgourgues, and T. Tram, “The Cosmic Linear Anisotropy Solving System
(CLASS) II: Approximation schemes,” JCAP 07 (2011) 034, arXiv:1104.2933
[[Link]].
[5] A. Spurio Mancini, D. Piras, J. Alsing, B. Joachimi, and M. P. Hobson, “CosmoPower:
emulating cosmological power spectra for accelerated Bayesian inference from next-generation
surveys,” Mon. Not. Roy. Astron. Soc. 511 (2022) no. 2, 1771–1788, arXiv:2106.03846
[[Link]].
[6] A. Nygaard, E. B. Holm, S. Hannestad, and T. Tram, “CONNECT: a neural network based
framework for emulating cosmological observables and cosmological parameter inference,”
JCAP 05 (2023) 025, arXiv:2205.15726 [[Link]].

– 20 –
 [7] S. Günther, L. Balkenhol, C. Fidler, A. R. Khalife, J. Lesgourgues, M. R. Mosbech, and R. K.
Sharma, “OLÉ — Online Learning Emulation in cosmology,” JCAP 09 (2025) 059,
arXiv:2503.13183 [[Link]].
[8] T. Auld, M. Bridges, M. P. Hobson, and S. F. Gull, “Fast cosmological parameter estimation
using neural networks,” Mon. Not. Roy. Astron. Soc. 376 (2007) L11–L15,
arXiv:astro-ph/0608174.
[9] W. A. Fendt and B. D. Wandelt, “Pico: Parameters for the Impatient Cosmologist,” Astrophys.
J. 654 (2006) 2–11, arXiv:astro-ph/0606709.
[10] A. Lazanu, “Power spectrum emulators from neural networks and tree-based methods,” Phys.
Dark Univ. 51 (2026) 102186, arXiv:2506.07514 [[Link]].
[11] K. Heitmann, E. Lawrence, J. Kwan, S. Habib, and D. Higdon, “The Coyote Universe
Extended: Precision Emulation of the Matter Power Spectrum,” Astrophys. J. 780 (2014) 111,
arXiv:1304.7849 [[Link]].
[12] Euclid Collaboration, M. Knabenhans et al., “Euclid preparation: II. The EuclidEmulator – A
tool to compute the cosmology dependence of the nonlinear matter power spectrum,” Mon.
Not. Roy. Astron. Soc. 484 (2019) 5509–5529, arXiv:1809.04695 [[Link]].
[13] G. Aricò, R. E. Angulo, and M. Zennaro, “Accelerating Large-Scale-Structure data analyses by
emulating Boltzmann solvers and Lagrangian Perturbation Theory,” arXiv:2104.14568
[[Link]].
[14] K. R. Moran, K. Heitmann, E. Lawrence, S. Habib, D. Bingham, A. Upadhye, J. Kwan,
D. Higdon, and R. Payne, “The Mira–Titan Universe – IV. High-precision power spectrum
emulation,” Mon. Not. Roy. Astron. Soc. 520 (2023) no. 3, 3443–3458, arXiv:2207.12345
[[Link]].
[15] L. Negri and A. Samajdar, “Neural likelihood estimators for flexible gravitational wave data
analysis,” arXiv:2509.17606 [[Link]].
[16] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. V. CMB power spectra and
likelihoods,” Astron. Astrophys. 641 (2020) A5, arXiv:1907.12875 [[Link]].
[17] BOSS Collaboration, S. Alam et al., “The clustering of galaxies in the completed SDSS-III
Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample,”
Mon. Not. Roy. Astron. Soc. 470 (2017) no. 3, 2617–2652, arXiv:1607.03155 [[Link]].
[18] F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker,
W. Saunders, and F. Watson, “The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the
Local Hubble Constant,” Mon. Not. Roy. Astron. Soc. 416 (2011) 3017–3032,
arXiv:1106.3366 [[Link]].
[19] A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, “The
clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at
z = 0.15,” Mon. Not. Roy. Astron. Soc. 449 (2015) no. 1, 835–847, arXiv:1409.3242
[[Link]].
[20] E. B. Wilson and M. M. Hilferty, “The distribution of chi-square,” Proceedings of the National
Academy of Sciences of the United States of America 17 (1931) no. 12, 684–688.
[21] J. Alsing, H. Peiris, J. Leja, C. Hahn, R. Tojeiro, D. Mortlock, B. Leistedt, B. D. Johnson, and
C. Conroy, “SPECULATOR: Emulating Stellar Population Synthesis for Fast and Accurate
Galaxy Spectra and Photometry,” The Astrophysical Journal Supplement Series 249 (2020) 5,
arXiv:1911.11778 [[Link]].
[22] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The MCMC Hammer,”
Publ. Astron. Soc. Pac. 125 (2013) 306–312, arXiv:1202.3665 [[Link]].

– 21 –
[23] S. Kullback and R. A. Leibler, “On Information and Sufficiency,” The Annals of Mathematical
Statistics 22 (1951) no. 1, 79–86.
[24] A. Nygaard, E. B. Holm, S. Hannestad, and T. Tram, “Fast and effortless computation of
profile likelihoods using CONNECT,” JCAP 11 (2023) 064, arXiv:2308.06379
[[Link]].
[25] S. Brooks, “Markov chain Monte Carlo method and its application,” Journal of the Royal
Statistical Society: Series D (The Statistician) 47 no. 1, 69–100,
[Link]
[Link]

– 22 –