Recipes#

Short user stories for the ground-up rxmc described in groundup_design.md. Each recipe says what a user wants, how it is spelled, and what behaviour they should expect. Together they are the minimal set of use cases the rewrite must support; every one of them is backed by a notebook or a test in the current repository, or by a design decision recorded in the design document. Every recipe is also a test: test/recipes/test_recipe_NN_<slug>.py asserts its Expected behaviour bullets, and the richest recipes are tutorial notebooks (design document §7 and §9). Snippets assume:

import numpy as np, rxmc as rx
from rxmc import terms as T, transforms as tf
from scipy import stats
from sklearn.gaussian_process.kernels import Matern, RBF

1. Fit a model to data with reported statistical errors#

I have x, y, and a statistical error per point, and a model with a few parameters. I want the posterior.

m, b = rx.Parameter("m", prior=stats.norm(0, 5)), rx.Parameter("b", prior=stats.norm(0, 5))
line = rx.Model(lambda x, m, b: m * x + b, [m, b])
d = rx.Dataset(x, y, y_err)
problem = rx.Problem([rx.Constraint([rx.Comparison(d, line)])])

Expected behaviour:

The covariance is diag(y_err**2) and nothing else. problem.chi2(theta) equals sum(((y - ym) / y_err)**2) exactly.
problem.ndim == 2, problem.names == ["m", "b"], in declaration order.
problem.log_posterior, problem.sample_prior(n), problem.prior_transform are all that emcee or dynesty need (recipe 16).

2. Infer an unknown noise level#

My data have no usable error bars, or I do not trust them. I want to infer the noise magnitude alongside the model.

log_eps = rx.Parameter("log_eps", prior=stats.norm(-2, 2))
c = rx.Constraint([rx.Comparison(d, line)], terms=[T.noise(log_eps)], statistical=False)
# or fractional noise:  T.noise_fraction(log_eps)
# or model error on the average of data and prediction:  T.model_error(log_gamma)
# or noise growing along x:  T.noise(log_eps, basis=T.exp_growth(np.pi), basis_params=(slope,))

Expected behaviour:

With statistical=False the inferred noise replaces the reported errors; with the default statistical=True it is added to them.
The posterior of log_eps reflects the residual scatter. In the sampling_algos scenario its truth is recovered.
noise_fraction and model_error scale with the prediction, so the covariance changes with the model parameters. That is allowed and costs nothing extra.

3. Use the reported systematic errors#

The measurement reports a fractional normalisation error and an absolute offset error. I want them in the likelihood as correlated modes.

d = rx.from_measurement(m, reaction=reaction, quantity="dXS/dA")   # norm_err, offset_err filled
comp = rx.Comparison(d, omp)
c = rx.Constraint([comp], terms=comp.reported_terms())

Expected behaviour:

reported_terms() returns the offset mode first, then the prediction-scaled normalisation mode. Zero magnitudes yield no term.
Without the call, neither systematic enters the covariance. Nothing is folded in silently.
Re-adding both terms reproduces the pre-0.1 auto-folded covariance to the last digit (test_regression, log-likelihood 1.195784087817536).
Under space=tf.log the modes are propagated by the delta method: the offset at the data, the normalisation at the prediction.

4. Infer a normalisation or offset the experiment did not report#

I suspect an unreported normalisation (or background offset) and want its magnitude as a nuisance parameter.

log_eta = rx.Parameter("log_eta", prior=stats.norm(-3, 1))
c = rx.Constraint([comp], terms=[T.normalization(parameter=log_eta)])
# absolute offset instead:  T.offset(parameter=log_omega)
# a mode with any shape:     T.systematic(log_s, basis=T.x_basis(np.pi))

Expected behaviour:

One rank-one mode exp(log_eta)**2 * outer(ym, ym) is added.
The model parameters decorrelate from the overall scale of the data; the data’s normalisation pull moves into log_eta.

6. One latent scale per dataset#

Each dataset has its own unknown normalisation. I want a Kennedy–O’Hagan scale factor on the model prediction, one per dataset, sampled with the model parameters.

rhos = [rx.Parameter(f"log_rho_{i}", prior=stats.norm(0, np.log1p(s))) for i, s in enumerate(sys_errs)]
comps = [rx.Comparison(d_i, omp | tf.scale(rho_i)) for d_i, rho_i in zip(datasets, rhos)]
# one global scale instead:  omp | tf.scale(rho) on every comparison

Expected behaviour:

The scale multiplies the mean, so it is not a covariance term. Its parameters follow the model’s in problem.names.
A masked or held-out view of the comparison keeps the same rho_i because the comparison, not the data, carries the model.
Fitting the ρᵢ recovers the true normalisations by MAP (normalization_inference); folding the reported σ_sys in as fixed normalisation modes (recipe 3) gives statistically the same model posterior.

7. Absorb model deficiency with a Gaussian process#

My model is missing physics. I want a smooth correlated discrepancy, in angle or in momentum transfer, learned from the residuals.

log_A = rx.Parameter("log_A", prior=stats.norm(0, 2))
gp = T.kernel(Matern(1.0, nu=2.5), on=comp, coords=lambda x: x / np.pi,
              amplitude=T.constant_amplitude, amplitude_params=(log_A,))
# in momentum transfer with amplitude A q^(r/2):
q = lambda x: rx.reactions.momentum_transfer(x, d.meta["k"])
gp_q = T.kernel(RBF(1.0), on=comp, coords=q, amplitude=lambda c, lA, r: np.exp(lA) * c.x ** (r / 2),
                amplitude_params=(log_A, r))
c = rx.Constraint([comp], terms=[gp])
band = rx.predictive.total_predictive_band(problem, gp, omp.bind(x_fine, d.meta), x_fine, samples)

Expected behaviour:

One parameter per free kernel hyperparameter element, named discrepancy_<hyperparameter>, sampled in sklearn’s log-theta space.
The model parameters relax from their biased values toward the truth (gp_discrepancy); the learned discrepancy tracks the true defect.
total_predictive_band finds the kernel’s columns from the term itself; no column arithmetic.
Every error-model form of the α+Ca study reproduces a hand-built dense matrix (TestStudyForms).

8. Sample a mean discrepancy explicitly#

Instead of marginalising the discrepancy, I want to sample an additive correction with a parametric shape.

phi = [rx.Parameter(f"c{i}", prior=stats.norm(0, 1)) for i in range(3)]
delta = rx.Model(lambda x, c0, c1, c2: c0 + c1 * x + c2 * x**2, phi)
comp = rx.Comparison(d, omp + delta)

Expected behaviour:

The prediction is omp(x) + delta(x) in physical space; problem.names lists omp.params then phi.
(omp + delta) | tf.scale(rho) scales the sum; (omp | tf.scale(rho)) + delta scales only the model.
A multiplicative, x-dependent correction is omp * Model(g_fn, phi); an additive discrepancy in log space is exactly that. scale(rho) is its constant special case.
The GP of recipe 7 may be used alongside it.

9. Heavy tails#

A few points are gross outliers. I do not want them to drag the fit; I want a likelihood that widens instead of breaking.

nu = rx.Parameter("nu", bounds=(1.0, 100.0))          # uniform on the bounds
c = rx.Constraint([comp], likelihood=rx.StudentT(nu=nu))

Expected behaviour:

Same covariance, different functional. The Gaussian version puts the truth at several σ; the Student-t version covers it (robust_likelihoods), and the posterior of nu is small.
The multivariate t applies one radial tail to the whole stacked residual of the constraint. A comparison that needs Gaussian tails goes in its own constraint.
rx.Chi2() drops the log-determinant for a pure chi-squared objective.
The Student-t value is pinned to its closed form (test_likelihood).

10. Compare in log space#

Cross sections span orders of magnitude. I want the Gaussian to live in log space, with the model still written in physical units.

comp = rx.Comparison(d, omp, space=tf.log)
c = rx.Constraint([comp], terms=[T.noise(log_eps)], statistical=False)   # constant noise in log space

Expected behaviour:

comp.y == log(d.y) and comp.y_err == d.y_err / d.y (delta method). The prediction is transformed the same way at every evaluation, so it can never be double-transformed.
A non-positive prediction gives log_likelihood == -inf and chi2 == +inf. A non-positive datum on an active point is a compile error naming the dataset.
problem.log_jacobian() is -sum(log y) over active points; add it to the log-evidence of the log-space fit before comparing with a linear-space fit of the same data.
Constant noise in log space and fractional noise in linear space are different models with different evidences. That is the point.

11. Hold out data and score it#

I want to fit below an angular cut and score the prediction above it, or compare error models by their held-out predictive density.

fit  = c.masked_where(lambda x: x < cut)
held = fit.complement()
problem = rx.Problem([fit], priors=priors)
heldout = rx.Problem([held], priors=problem.priors)
lp = rx.diagnostics.heldout_log_predictive(heldout, samples)
score = rx.diagnostics.log_posterior_predictive(lp, logw=res.logwt)   # weighted for nested sampling

Expected behaviour:

fit and held share every Comparison, Term, and Parameter; the two problems have identical names, so a chain from one scores the other.
Active sets are disjoint, their union is every point, and ll(fit) + ll(held) == ll(full) for a block-local covariance.
Terms are authored once over all points; masking selects rows, it never rebuilds anything.

12. Temper the likelihood#

I have many points and worry the posterior is overconfident, or I want a power posterior.

c = rx.Constraint([comp], weight=2 / n)

Expected behaviour:

The log-likelihood of that constraint is multiplied by weight; the prior is never touched. One knob, honoured by every driver.
Empirical coverage of the posterior predictive (recipe 17) moves toward the nominal level as the weight decreases (overconfidence).

13. Declare priors#

I want each nuisance parameter to carry its own prior, the optical potential to have a correlated multivariate normal prior, and everything to work with nested sampling.

log_eps = rx.Parameter("log_eps", prior=stats.halfnorm(scale=1))            # marginal
nu      = rx.Parameter("nu", bounds=(1, 100))                               # uniform on bounds
V       = rx.Parameter("V", prior=stats.norm(50, 5), bounds=(30, 70))       # truncated marginal
problem = rx.Problem([c], priors=[(omp.params, stats.multivariate_normal(mu, cov))])   # joint block

Expected behaviour:

Every slot is covered exactly once, or Problem raises naming the parameter (uncovered) or the pair of priors (double covered).
problem.prior_transform(u) exists for all of the above: rescaled ppf for marginals, whitening for the multivariate normal. It is finite at u = 0 and u = 1.
problem.sample_prior(n) gives an (n, ndim) array suitable as emcee’s starting positions.

14. From an EXFOR measurement to a dataset#

I have an exfor_tools distribution in mb/sr, or as a ratio to Rutherford, with its reported errors. I want a dataset in the library’s units with nothing lost.

d = rx.from_measurement(m, reaction=reaction, quantity="dXS/dA")      # or "dXS/dRuth", "Ay"
d_ias = rx.from_measurement(m, reaction=reaction, ExIAS=Ex)           # (p,n) IAS channel

Expected behaviour:

y, y_err and offset_err are divided by the unit conversion factor; norm_err (fractional) passes through untouched. Angles are stored in radians.
Requesting dXS/dRuth from a dXS/dA measurement (or the reverse) uses the Rutherford cross section on the data angles, a closed form of the kinematics; the factor is then a per-angle array.
d.meta carries reaction, Elab, quantity, k (and ExIAS), which is everything a reaction model needs to bind.
Incompatible units, or a quantity the measurement cannot be converted to, raise at conversion time.

15. Evaluate a reaction model on any grid#

I want the model on the data angles for the likelihood and on a fine grid for plotting, with the solver set up once per grid.

omp = rx.reactions.ElasticXS("dXS/dA", central, spin_orbit, args_from_params, params,
                             lmax=20, wavelengths_beyond_range=2.0, zeros_per_node=5)
comp = rx.Comparison(d, omp)                                   # bound to d.x, d.meta
fine = omp.bind(np.deg2rad(np.linspace(0.5, 179.5, 200)), d.meta)
ys = [fine(*s[problem.columns(omp.params)]) for s in samples[::50]]
band = rx.predictive.predictive_band(ys, levels=(5, 50, 95))

Expected behaviour:

The model owns the jitr workspace; the dataset never does. Solver settings on the model reach jitr (TestSolverSettingsForwarding).
The predictor is a pure function of the potential. A compound-elastic contribution is subtracted from d.y before the dataset is built (recipe 20); there is no correction hook on the model.
Cross sections are in b/sr; dXS/dRuth and Ay are dimensionless.
problem.columns(omp.params) is the only chain bookkeeping a user does.

16. Drive the calibration with an external sampler#

I want to use emcee, dynesty, or the black-box-bayes CLI, and read the chain back without positional arithmetic.

# emcee
p0 = problem.sample_prior(32)
sampler = emcee.EnsembleSampler(32, problem.ndim, problem.log_posterior)
samples = sampler.run_mcmc(p0, 5000) and sampler.get_chain(discard=1000, flat=True)

# dynesty
ns = dynesty.NestedSampler(problem.log_likelihood, problem.prior_transform, problem.ndim)
ns.run_nested(); samples = ns.results.samples_equal()
logz = rx.diagnostics.logz_summary(ns.results.logz[-1], ns.results.logzerr[-1])

# black-box-bayes
dill.dump(problem, open("problem.pkl", "wb"))        # posterior.py delegates six names to it

Expected behaviour:

Chains from all three are (n, ndim) arrays in problem.names order. Every analysis function takes (problem, samples).
log_posterior evaluates the prior first and never calls the forward model when the prior is -inf.
dill.loads(dill.dumps(problem)).log_posterior(theta) equals problem.log_posterior(theta), workspaces included.
problem.NDIM, problem.parameter_names, problem.starting_location, problem.log_posterior_batch are the bbb spellings of the same things.

17. Check the posterior predictive#

I want to know whether my error model is calibrated: do 68 % intervals contain 68 % of the points, and how wide are they?

draws = rx.diagnostics.predictive_draws(problem, samples, constraint=0, n_rep=4)
cov = rx.diagnostics.coverage_curve(draws, problem.constraints[0].y[problem.constraints[0].active])
err = rx.diagnostics.coverage_error(draws, y_active)
width = rx.diagnostics.sharpness(draws, transform=np.exp)      # widths in physical space for a log fit

Expected behaviour:

Draws are ym(theta) + L z on the active points in comparison space; model_only=True returns ym(theta) and assembles no covariance.
Coverage is near nominal for a correct error model and clearly below it for an overconfident one.
logz_summary reports the max of the replicate half-range and the sampler’s own error; compare_logz returns "tie" unless the difference exceeds sigma * hypot(err_a, err_b).

18. Compare error models by evidence#

The α+Ca study: several error models for data without reported errors. I want the evidence for each, comparable across comparison spaces.

models = {
    "L0": rx.Constraint([comp_log], terms=[T.noise(log_eps)], statistical=False),
    "E0": rx.Constraint([comp_lin], terms=[T.noise_fraction(log_eps)], statistical=False),
    "L2y": rx.Constraint([comp_log], terms=[T.noise(log_eps), T.normalization(log_sys)], statistical=False),
    "Lgp": rx.Constraint([comp_log], terms=[T.noise(log_eps), gp], statistical=False),
    "L0t": rx.Constraint([comp_log], terms=[T.noise(log_eps)], statistical=False, likelihood=rx.StudentT()),
}
logz = {}
for name, c in models.items():
    p = rx.Problem([c.masked_where(lambda x: x < cut)], priors=priors)
    res = run_dynesty(p)
    logz[name] = rx.diagnostics.logz_summary(res.logz[-1] + p.log_jacobian(), res.logzerr[-1])
verdict = rx.diagnostics.compare_logz(logz["Lgp"], logz["L0"])

Expected behaviour:

The same log_eps object is reused across models without conflict: each Problem is compiled independently.
Adding log_jacobian makes the log-space and linear-space evidences comparable.
Every term list in the table reproduces the corresponding hand-built covariance in TestStudyForms.

19. Bring my own covariance or term#

I have a full covariance matrix from a correlated measurement, or a noise model no helper expresses.

fixed = rx.Term(C, on=comp)                                       # any symmetric PD matrix
stat  = rx.Term(sig, kind="diag", on=comp)                        # a fixed std-dev vector
custom = rx.Term(lambda c, e, l: np.exp(e) * np.exp(l * c.x / np.pi), (log_e, slope), kind="diag")

Expected behaviour:

A plain array is a fixed contribution, factored once. Shape and symmetry are checked at construction against the term’s on.
A callable sees a TermContext with x (through coords), y, ym, and len(c); it returns a vector for diag/mode or a matrix.
Fitting correlated data with the correct Term(C) instead of its diagonal is the difference between an honest and an overconfident posterior (normalization_inference gallery).

20. Preprocess instead of asking for a feature#

I want to mean-subtract, standardise, or project both data and model onto principal components of a prior predictive ensemble.

mu = ens.mean(0); A = np.linalg.svd(ens - mu, full_matrices=False)[2][:k]
d_pc = rx.Dataset(np.arange(k), A @ (d.y - mu), np.zeros(k), label=f"{d.label} PC")
native = omp.bind(d.x, d.meta)
proj = rx.Model(lambda x, *theta: A @ (native(*theta) - mu), omp.params)
stat = rx.Term(A @ np.diag(d.y_err**2) @ A.T, on=d_pc)
c = rx.Constraint([rx.Comparison(d_pc, proj)], terms=[stat, T.noise(log_eps, on=d_pc)], statistical=False)

Expected behaviour:

Dataset.x is opaque to the library; a model may ignore its x and close over a predictor on another grid; y_err may be zero when statistical=False. None of these is “fixed” later.
Subtracting a known contribution from the data, such as a compound- elastic cross section, is the same pattern: replace(d, y=d.y - cn). A reported fractional normalisation error refers to the measured value, so build that mode from the unsubtracted prediction if it matters; the offset error is unaffected.
Pointwise maps (mean subtraction, standardisation) can equally be a space= transform built from closed-over arrays, which keeps the delta method and log_jacobian correct.
Terms that scale with the native prediction see the projected ym; they must call the native predictor themselves if they need it.
A projected fit’s evidence is not comparable to the unprojected one.

21. Get a useful error, not a `LinAlgError`#

My measurement reports no statistical error and I forgot to add any term.

d = rx.from_measurement(m_without_errors, reaction=reaction, quantity="dXS/dA")
rx.Problem([rx.Constraint([rx.Comparison(d, omp)])])
# ValueError: covariance of constraint 0 is singular on comparison 'E1234-002': the dataset
# reports zero statistical error and no term covers its points.  Add
# comp.reported_terms(), a noise term, or a fixed Term; or set statistical=False and
# compose the covariance explicitly.

Expected behaviour:

Raised by Problem, before any sampler runs, naming the comparison by its label. A parametric covariance is not checked eagerly (it may be fine at some parameter values).
The same compile step reports duplicate parameter names, an uncovered prior slot, an on= outside its constraint, and a non-finite comparison-space value on an active point.

23. A discrepancy correlated across energies and angles#

I believe the model’s defect varies smoothly in both energy and angle. I want one GP over (E, θ) that correlates the datasets at different energies.

kE, kθ = RBF(10.0), Matern(0.3, nu=2.5)
lE, lθ, log_A = (rx.Parameter(n, prior=stats.norm(0, 1)) for n in ("log_lE", "log_ltheta", "log_A"))

def fn(c, lE, lθ, lA):
    E, θ = c.meta("Elab")[:, None], c.x[:, None]
    K = kE.clone_with_theta([lE])(E) * kθ.clone_with_theta([lθ])(θ)     # separable product
    return np.exp(2 * lA) * K + 1e-10 * np.eye(len(c))

md = rx.Term(fn, (lE, lθ, log_A), kind="matrix", on=comps)
c = rx.Constraint(comps, terms=[md])

Expected behaviour:

The term spans comparisons, so the constraint’s covariance has non-zero off-diagonal blocks between energies and takes the dense path. Cost grows as the cube of the total number of points; per-comparison GPs with a linked amplitude (recipe 22) are the cheap alternative when cross-energy correlation is not needed.
Any kernel built from c.meta and c.x is allowed; the product of two one-dimensional kernels is the factorised form. There is no Kronecker speed-up because datasets share no angle grid.
The dense matrix equals the hand-built product kernel on the stacked (E, θ) inputs.

24. Per-dataset parameters drawn from a sampled hyperprior#

Each dataset has its own normalisation, and I want to learn how spread out those normalisations are, rather than fix the spread.

rhos = [rx.Parameter(f"log_rho_{i}") for i in range(len(comps))]    # no marginal: the joint covers them
log_tau = rx.Parameter("log_tau")

class RhoHierarchy:
    def logpdf(self, v):                       # v ordered as (rhos..., log_tau)
        *r, lt = v
        return (stats.norm(0, np.exp(lt)).logpdf(r).sum()
                + stats.halfnorm(scale=0.3).logpdf(np.exp(lt)) + lt)    # hyperprior on tau, with Jacobian
    def prior_transform(self, u):              # hyperparameter first, children given it
        lt = np.log(stats.halfnorm(scale=0.3).ppf(u[-1]))
        return np.append(stats.norm(0, np.exp(lt)).ppf(u[:-1]), lt)

comps = [rx.Comparison(d_i, omp | tf.scale(rho_i)) for d_i, rho_i in zip(datasets, rhos)]
problem = rx.Problem([rx.Constraint(comps)], priors=[(rhos + [log_tau], RhoHierarchy())])

Expected behaviour:

log_tau is an ordinary column of the chain. Its posterior is the learned spread of the normalisations; the rho_i shrink toward zero when tau is small.
Forgetting the joint block is a compile error: the rho_i carry no marginal and have infinite bounds, so no prior covers them.
Nested sampling works because the comparison supplies its own unit-cube map, drawing the hyperparameter before the children.
The same mechanism gives per-dataset GP amplitudes with a shared spread, or any other partially pooled parameter.

25. SafeBayes: learn the tempering exponent#

I suspect my model is misspecified and want the tempering exponent η chosen by the data rather than by hand, following Grünwald’s SafeBayes (arXiv:1412.3730): minimise the prequential log-loss of the η-generalised posterior over prefixes of the data.

from dataclasses import replace
order = rng.permutation(n)                                   # one ordering; average a few
prefix = lambda i: [np.isin(np.arange(n), order[:i])]
etas = [1.0, 0.5, 0.25, 0.125]
loss = dict.fromkeys(etas, 0.0)

for i in range(n0, n):
    before = rx.Problem([c.masked(prefix(i))], priors)       # untempered, for scoring
    after  = rx.Problem([c.masked(prefix(i + 1))], priors)
    for eta in etas:
        post = rx.Problem([replace(c.masked(prefix(i)), weight=eta)], priors)
        samples = run(post)                                   # emcee or dynesty
        ll_i = np.array([after.log_likelihood(s) - before.log_likelihood(s) for s in samples])
        loss[eta] -= ll_i.mean()                              # R-log-loss

eta_hat = min(loss, key=loss.get)
final = rx.Problem([replace(c, weight=eta_hat)], priors)

Expected behaviour:

replace(c, weight=eta) and c.masked(...) share every comparison, term and parameter with c, so no workspace is rebuilt across the n × |η| problems and every problem has the same names.
after.log_likelihood(θ) − before.log_likelihood(θ) is log p(y_i | y_<i, θ) exactly, also under a correlated covariance. For a diagonal covariance it reduces to the single-point marginal, which is the paper’s setting.
With dynesty one nested run per prefix serves every η: the η-posterior weights over the same samples are logvol + η·logl. The cost is then n runs, not n × |η|.
weight is a float by type. Putting a prior on η and sampling it is not SafeBayes and is not coherent, since p(D|θ)^η is not a normalised likelihood in η; the design refuses it.
The same loop at fixed η with the held-out score of recipe 11 is the cheaper “choose η by held-out predictive”; SafeBayes is its principled form.

The recipes below come from a survey of the calibration literature, general and nuclear. Each is a pure user-level pattern over the skeleton as it stands; none needs a library addition. They are ordered by how often the pattern appears in nuclear-physics calibration practice. A closing section lists what the API does not express.

26. Unaccounted-for model error per data type (KDUQ)#

I am calibrating a global optical potential to many datasets of several observable types. I want one fractional “unaccounted-for” uncertainty per type, sampled with the potential, added in quadrature to the reported errors and scaled with the average of datum and prediction.

delta = {t: rx.Parameter(f"delta_{t}", prior=stats.halfnorm(scale=s0[t])) for t in ("dxs", "ay", "sig_tot")}
comps = [rx.Comparison(d, omp_for(d)) for d in datasets]
terms = [T.model_error(delta[d.meta["type"]], averaging=True, on=comp) for d, comp in zip(datasets, comps)]
c = rx.Constraint(comps, terms=terms)          # statistical=True: reported errors are a floor

# KDUQ additionally scales the whole log-likelihood by k/N ("democratic"), or
# each data type's share by k/(n_types N_t) ("federal").  That is tempering:
c_dem = rx.Constraint(comps, terms=terms, weight=n_params / n_data)
c_fed = [rx.Constraint([comp for comp in comps if comp.data.meta["type"] == t],
                       terms=[tt for tt in terms if tt.on.data.meta["type"] == t],
                       weight=n_params / (len(types) * n_pts[t]))
         for t in types]                                     # one constraint per type, one delta_t each

Expected behaviour:

One delta_t column per type in problem.names, shared by every comparison of that type through the same Parameter object.
Σ_ii = y_err_i² + (δ_T (y_i + ym_i) / 2)². Because the term scales with the prediction it changes with the model parameters; that is allowed.
The democratic and federal variants are Constraint(weight=). The published text describes them as a rescaling of the covariance, which would also change logdet; the authors have confirmed that what was actually done is a scaling of the whole log-likelihood by k/N, which is tempering. The federal form is one constraint per data type, holding every comparison of that type, with its own weight; each delta_t lives in its type’s constraint, and only the potential parameters are common to all constraints, as in any multi-constraint problem.
Reported uncertainties are respected, then augmented; a dataset’s known normalisation systematic can still be added as a mode (recipe 3).

Reference: Pruitt, Escher, Rahman, Uncertainty-quantified phenomenological optical potentials for single-nucleon scattering, Phys. Rev. C 107, 014602 (2023), arXiv:2211.07741.

27. Peelle’s Pertinent Puzzle: normalise the prediction, never the data#

My datasets carry a common fractional normalisation error and I want the fit not to be biased low.

# right: the mode is built from the prediction (this is what normalization() does)
c_ok = rx.Constraint([comp], terms=[T.normalization(magnitude=d.norm_err, on=comp)])

# wrong on purpose, to see the bias: the mode built from the data
c_bad = rx.Constraint([comp], terms=[rx.Term(lambda c: d.norm_err * c.y, kind="mode", on=comp, constant=True)])

# t0 variant: a fixed mode from a reference prediction, refit once
t0 = comp.predictor(*theta_hat)
c_t0 = rx.Constraint([comp], terms=[rx.Term(d.norm_err * comp.space(t0), kind="mode", on=comp)])

Expected behaviour:

With the data-built mode, a fit of a constant to n points with fractional normalisation error s is biased low by the factor 1 / (1 + n s²), growing without bound in n. This is D’Agostini’s bias and the origin of Peelle’s Pertinent Puzzle. With the prediction-built mode the estimate is unbiased; an additive offset mode has no such bias either way.
normalization() reads c.ym, so the default spelling is the safe one. A free log_eta (recipe 4) also multiplies the prediction.
The t0 mode makes the covariance constant, so it is factored once; one refit after a zeroth-order fit suffices in practice.

References: D’Agostini, On the use of the covariance matrix to fit correlated data, Nucl. Instrum. Meth. A 346, 306 (1994); Frühwirth, Neudecker, Leeb, Peelle’s Pertinent Puzzle and its solution, EPJ Web Conf. 27, 00008 (2012); Ball et al. (NNPDF), Fitting parton distribution data with multiplicative normalization uncertainties, JHEP 05 (2010) 075.

28. Stacking by leave-one-dataset-out#

Evidence weights assume the true model is among my candidates. I would rather weight models by how well they predict each dataset when it is left out.

def loo_scores(model):
    out = []
    for i in range(len(comps)):
        fit_c = c.masked([np.full(comp.data.y.shape, j != i) for j, comp in enumerate(comps)])
        fit, held = rx.Problem([fit_c], priors), rx.Problem([fit_c.complement()], priors)
        s = run(fit)
        out.append(rx.diagnostics.log_posterior_predictive(rx.diagnostics.heldout_log_predictive(held, s)))
    return np.array(out)                                   # one log score per held-out dataset

S = np.stack([loo_scores(m) for m in models])              # (n_models, n_datasets)
w = maximise(lambda w: np.sum(logsumexp(np.log(w)[:, None] + S, axis=0)), simplex)

Expected behaviour:

Held-out log densities are joint over the held-out comparison and exact under a correlated covariance; PSIS-LOO per point is not available without per-point likelihood factors (closing section).
Stacking weights need not sum to the evidence weights; in the M-open setting they are the ones to prefer.
The cost is one fit per dataset per model; for a diagonal covariance and many datasets the importance-sampling shortcut of Vehtari et al. is the cheaper path, implemented from problem.predict and y_err by the user.

References: Yao, Vehtari, Simpson, Gelman, Using stacking to average Bayesian predictive distributions, Bayesian Analysis 13, 917 (2018); Vehtari, Gelman, Gabry, Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC, Stat. Comput. 27, 1413 (2017).

29. Cut (modular) posterior by multiple imputation#

One module of my model, say a systematic-error parameter or a GP hyperparameter, should be learned from its own data only and not be contaminated by the primary data, which I trust less.

stage1 = rx.Problem([rx.Constraint([rx.Comparison(d_aux, aux_model)])], priors_aux)     # phi from Z only
phis = run(stage1)[::thin][:T]

pooled = []
for phi in phis:
    model_t = rx.Model(functools.partial(f, phi=phi), theta_params)                  # phi fixed
    stage2 = rx.Problem([rx.Constraint([rx.Comparison(d, model_t)], terms=terms_at(phi))], priors_theta)
    pooled.append(run(stage2))
samples = np.concatenate(pooled)

# partial feedback instead of a cut: per-module powers
c_aux = rx.Constraint([rx.Comparison(d_aux, aux_model)], weight=gamma1)
c_pri = rx.Constraint([rx.Comparison(d, full_model)], weight=gamma2)

Expected behaviour:

The marginal for phi equals its stage 1 posterior; the primary data never update it. Uncertainty in phi still propagates into theta.
The cut distribution is not a single posterior, so it cannot be handed to emcee or dynesty as one log_posterior; the two-stage loop is the definition, not an approximation. A single alternating chain converges to updater-dependent limits and is not a substitute.
The power-weighted alternative is one Problem with two weights and is a legitimate posterior; it gives partial feedback control.
Bayarri’s “modular MLE” is the same shape with T = 1: fix GP hyperparameters at their stage 1 estimates and pass them as constants.

References: Plummer, Cuts in Bayesian graphical models, Stat. Comput. 25, 37 (2015); Jacob, Murray, Holmes, Robert, Better together? Statistical learning in models made of modules, arXiv:1708.08719; Bayarri et al., A framework for validation of computer models, Technometrics 49, 138 (2007).

30. Leave-one-experiment-out prediction#

I want to know whether the calibrated model, with its discrepancy, predicts an experiment it was not fit to, and with what tolerance.

for i, comp in enumerate(comps):
    fit  = c.masked([np.full(other.data.y.shape, j != i) for j, other in enumerate(comps)])
    p    = rx.Problem([fit], priors)
    s    = run(p)
    held = rx.Problem([fit.complement()], priors)
    draws = rx.diagnostics.predictive_draws(held, s, n_rep=4)
    tol = np.percentile(np.abs(draws - draws.mean(0)), 90, axis=0)      # tolerance bound per point
    cov = rx.diagnostics.coverage_curve(draws, held.constraints[0].y[held.constraints[0].active])

Expected behaviour:

The held-out comparison’s covariance terms are the same objects as in the fit; a GP discrepancy conditioned on the other experiments carries into the prediction through predictive_draws.
Coverage on the held-out experiment is the honest check; in-sample coverage is not.
Tolerance bounds are empirical percentiles of the draws, componentwise.

References: Higdon, Gattiker, Williams, Rightley, Computer model calibration using high-dimensional output, J. Am. Stat. Assoc. 103, 570 (2008); Bayarri et al., Technometrics 49, 138 (2007).

31. Simulation-based calibration of the sampler#

Before trusting a chain, I want to check that the sampler recovers parameters drawn from the prior when the data are simulated from the model.

ranks = []
for theta0 in problem.sample_prior(n_sims, rng):
    y_sim = rx.diagnostics.predictive_draws(problem, theta0[None], n_rep=1)[0]
    d_sim = dataclasses.replace(d, y=comp.space.inverse(y_sim))            # back to physical units
    p_sim = rx.Problem([rx.Constraint([rx.Comparison(d_sim, model)], terms=terms)], priors)
    s = run(p_sim)[::thin][:L]
    ranks.append((s < theta0).sum(0))                                 # one rank per column
ranks = np.array(ranks)                                               # uniform on 0..L if correct

Expected behaviour:

Every rank histogram is uniform when the sampler is correct. An inverted U means the computed posterior is too wide, a U too narrow, a skew a bias.
Chains must be thinned to roughly independent draws first, or spurious boundary spikes appear.
SBC validates the computation under the assumed model; it says nothing about whether the model fits real data; that is the posterior predictive coverage check of recipe 17.

Reference: Talts, Betancourt, Simpson, Vehtari, Gelman, Validating Bayesian inference algorithms with simulation-based calibration, arXiv:1804.06788.

32. Emulator as the model, emulator variance as a term#

My model is too expensive to run in the chain. I have a GP or PCA emulator trained on a design of runs, and I want its predictive variance in the likelihood.

emu = rx.Model(lambda x, *theta: emulator.mean(theta), params)             # x ignored: the grid is the design's
emu_var = rx.Term(lambda c, *theta: np.sqrt(emulator.var(theta)), params, kind="diag", on=comp)
c = rx.Constraint([rx.Comparison(d, emu)], terms=[emu_var, T.noise(log_eps)])

Expected behaviour:

The term declares the same Parameter objects as the model, so it receives the current theta and can evaluate the emulator variance there. No special mechanism.
Bayarri et al. recommend fixing emulator hyperparameters at their estimates from the design runs (recipe 29 with T = 1) because emulator uncertainty is usually dominated by calibration and bias uncertainty.
The only full-posterior nuclear EDF calibration to 2015 replaced the code by a GP response surface exactly this way.

References: Higdon et al., J. Am. Stat. Assoc. 103, 570 (2008); McDonnell, Schunck, Higdon, Sarich, Wild, Nazarewicz, Phys. Rev. Lett. 114, 122501 (2015); Schunck et al., Uncertainty quantification and propagation in nuclear density functional theory, Eur. Phys. J. A 51, 169 (2015); Bayarri et al., Technometrics 49, 138 (2007).

33. MAP and Laplace approximation#

I want a quick Gaussian approximation to the posterior, and to know when it is good enough.

from scipy.optimize import minimize
res = minimize(lambda t: -problem.log_posterior(t), problem.sample_prior(1)[0], method="L-BFGS-B",
               bounds=problem.bounds)
H = numerical_hessian(lambda t: -problem.log_posterior(t), res.x)
cov = np.linalg.inv(H)                                   # Laplace covariance

Expected behaviour:

For a near-Gaussian posterior (a five-parameter optical potential with a flat prior) the Laplace covariance reproduces the emcee uncertainties. For a six-parameter potential the posterior is banana-shaped and it does not; compare against a chain before reporting.
This is the dominant uncertainty-propagation practice in nuclear DFT (inverse Hessian at the optimum, propagated linearly).
problem.bounds feeds the optimiser; the prior handles the rest.

References: Pruitt, Lovell, Hebborn, Nunes, The role of the likelihood for elastic scattering uncertainty quantification, arXiv:2403.00753; Schunck et al., Eur. Phys. J. A 51, 169 (2015).

34. Global error scale factor and unrecognised sources of uncertainty#

Repeated measurements scatter more than their stated errors. I want a global scale on the reported errors, or a fully correlated unknown component per experimental technique.

log_s = rx.Parameter("log_s", prior=stats.norm(0, 0.5))
scaled = rx.Term(lambda c, ls: np.exp(ls) * comp.y_err, (log_s,), kind="diag", on=comp)   # closes over the comparison's errors
c = rx.Constraint([comp], terms=[scaled], statistical=False)               # s multiplies every stated error

# USU: one unknown, fully correlated component per technique, shared across its datasets
log_delta = {tech: rx.Parameter(f"log_usu_{tech}") for tech in techniques}
usu = [T.offset(parameter=log_delta[tech], on=[comp for comp in comps if comp.data.meta["technique"] == tech])
       for tech in techniques]

Expected behaviour:

Under a Student-t likelihood only a modest scale is needed where least squares needs a large one, because the tail absorbs the outliers.
A global Birge-type rescaling inflates every point equally and is judged poor evaluation practice; the targeted USU component is preferred, added only when other explanations are exhausted.
A USU component inside the fit shifts the evaluated means, not only the widths, whenever more than one quantity is evaluated. It must be in the covariance, not added afterwards, which is why it is a sampled term here.
With statistical=False the scaled term is the diagonal; it closes over comp.y_err, the comparison-space errors, so it is right under a space transform too.

References: Hanson, Lessons about likelihood functions from nuclear physics, AIP Conf. Proc. 954 (2007), arXiv:0712.0021; Capote et al., Unrecognized sources of uncertainties (USU) in experimental nuclear data, Nucl. Data Sheets 163, 191 (2020), arXiv:1911.01825.

35. Energy-dependent parameters and per-comparison model instances#

A potential depth depends on energy through a few coefficients I want to share across datasets at different energies.

V0, V1 = rx.Parameter("V0", prior=stats.norm(50, 5)), rx.Parameter("V1", prior=stats.norm(-0.3, 0.1))

def depth_model(E):
    return rx.Model(lambda x, v0, v1, *rest: potential(x, V=v0 + v1 * E, *rest), [V0, V1, *rest_params])

comps = [rx.Comparison(d, depth_model(d.meta["Elab"])) for d in datasets]

Expected behaviour:

One model instance per comparison, all sharing the coefficient objects, so problem.names has one V0 and one V1.
Energy enters by closure at construction; nothing in the library reads it. The same pattern gives energy-dependent systematic errors in a term through c.meta("Elab") (recipe 22).
A smooth energy dependence with more freedom is a discrepancy on a basis (recipe 36) or a GP over energy (recipe 23).

References: Schnabel, Capote, Koning, Brown, Nuclear data evaluation with Bayesian networks, arXiv:2110.10322; Pruitt, Escher, Rahman, Phys. Rev. C 107, 014602 (2023).

36. Discrepancy on a physically constrained basis#

I know the shape the model defect can take, say a few Legendre modes in angle, and want the discrepancy restricted to that basis.

from scipy.special import eval_legendre
modes = [T.systematic(rx.Parameter(f"log_s{k}", prior=stats.norm(-3, 1)),
                      basis=lambda c, k=k: eval_legendre(k, np.cos(c.x)), on=comp) for k in range(1, 4)]
c_marg = rx.Constraint([comp], terms=modes)                        # marginalised: rank-3 covariance

coeffs = [rx.Parameter(f"c{k}", prior=stats.norm(0, 0.1)) for k in range(1, 4)]
delta = rx.Model(lambda x, *cs: sum(ck * eval_legendre(k, np.cos(x)) for k, ck in enumerate(cs, 1)), coeffs)
c_samp = rx.Constraint([rx.Comparison(d, omp + delta)])              # sampled: mean correction

Expected behaviour:

Each mode is one rank-one term; the covariance is low rank and stays on the structured path.
The sampled form gives the coefficients’ posterior directly; the marginalised form gives their scales. Both confound with the model parameters if the basis overlaps the model’s own response, so put a real prior on the amplitudes.
Process-convolution and kernel bases are the same pattern with a different basis.

Reference: Higdon, Gattiker, Williams, Rightley, J. Am. Stat. Assoc. 103, 570 (2008).

37. Correlated systematics between observables of one measurement#

One experiment reports both a cross section and an analysing power, and they share a normalisation or an angle calibration.

comp_xs = rx.Comparison(d_xs, rx.reactions.ElasticXS("dXS/dA", *pot, params=p))
comp_ay = rx.Comparison(d_ay, rx.reactions.ElasticXS("Ay", *pot, params=p))
log_eta = rx.Parameter("log_eta", prior=stats.norm(-3, 1))
c = rx.Constraint([comp_xs, comp_ay], terms=[T.normalization(log_eta, on=comp_xs),        # the ratio is unaffected
                                       T.systematic(log_dtheta, basis=dydtheta, on=[comp_xs, comp_ay])])

Expected behaviour:

Two datasets, two comparisons, one constraint; the shared systematic is a mode spanning both comparisons (case A of recipe 5).
A normalisation error affects the cross section and not a ratio observable; an angle-calibration error affects both through their angular derivatives, which the basis supplies from c.ym and c.x.
The multi-quantity extension of the Peelle treatment applies: build the mode from predictions, not data.

Reference: Neudecker, Frühwirth, Kawano, Leeb, Adequate treatment of correlated experimental data in nuclear data evaluations, Nucl. Data Sheets 118, 364 (2014).

38. The classic normal hierarchical model (eight schools)#

Several groups each report an estimate y_j with a known standard error σ_j. I believe the group effects θ_j are drawn from a common distribution N(μ, τ²) and want to learn μ, τ, and the shrunken θ_j.

d = rx.Dataset(x=np.arange(J), y=estimates, y_err=std_errors)          # x is the group index
mu      = rx.Parameter("mu", prior=stats.norm(0, 25))
log_tau = rx.Parameter("log_tau", prior=stats.norm(1, 1))

# 1. marginalised: integrate the theta_j out; y_j ~ N(mu, sigma_j^2 + tau^2)
c_marg = rx.Constraint([rx.Comparison(d, rx.Model(lambda x, mu: np.full(len(x), mu), [mu]))],
                       terms=[T.noise(log_tau)])
# theta_j afterwards, from the conditional normal at each draw (mu, tau):
#   mean = (y_j / s_j^2 + mu / tau^2) / (1 / s_j^2 + 1 / tau^2), var = 1 / (1 / s_j^2 + 1 / tau^2)

# 2. non-centred: theta_j = mu + tau * eta_j, eta_j ~ N(0, 1); marginal priors only
etas = [rx.Parameter(f"eta_{j}", prior=stats.norm(0, 1)) for j in range(J)]
school = rx.Model(lambda x, mu, lt, *eta: mu + np.exp(lt) * np.asarray(eta), [mu, log_tau, *etas])
c_nc = rx.Constraint([rx.Comparison(d, school)])

# 3. centred: theta_j as parameters, hyperprior as a joint block including mu, tau (recipe 24)
thetas = [rx.Parameter(f"theta_{j}") for j in range(J)]
c_c = rx.Constraint([rx.Comparison(d, rx.Model(lambda x, *th: np.asarray(th), thetas))])
problem_c = rx.Problem([c_c], priors=[(thetas + [mu, log_tau], SchoolHierarchy())])

Expected behaviour:

All three spellings give the same posterior for mu and tau. The posterior of tau piles up near zero, as in the book, and the theta_j shrink toward mu as tau falls.
The marginalised form has two columns and is the design’s own philosophy: a Gaussian latent belongs in the covariance. The non-centred form has J + 2 columns, all with marginal priors, so prior_transform is available with no joint block. The centred form has the funnel geometry the book warns about and samples worst with emcee or dynesty.
A parameter attached only to a fully masked comparison is sampled from its prior. Holding out one school therefore leaves its eta_j in the chain, driven by tau alone, so complement() and predictive_draws give the predictive for a new group with no extra code.
Meta-analysis on log-odds with per-study standard errors is the same structure. The beta-binomial rat-tumour example is not expressible: a binomial likelihood is not elliptical (see the closing section). A normal approximation with y_err = sqrt(n_j p_j (1 - p_j)) as a prediction-scaled diag term is expressible but is not the book’s model.

The notebook hierarchical_calibration (design document §7) is the parameter-space version of this hierarchy: per-dataset deviation vectors on the physics parameters with a hyperprior on their spread, used to repair a misspecified energy dependence.

Reference: Gelman, Carlin, Stern, Dunson, Vehtari, Rubin, Bayesian Data Analysis, 3rd ed., CRC Press (2013), Chapter 5.

What this API does not express#

Each item names the assumption that breaks, the nearest workaround, and the size of the addition that would lift it.

Non-elliptical likelihoods. Poisson counts, censored points and upper limits, and two-component good/bad mixtures (1 − β) N + β t (Hanson 2007) are not functionals of (d2, logdet, n). Workaround: none that is faithful. Addition: a Likelihood that receives the residual and the factor, plus a second evaluation path in CompiledConstraint; about 100 lines.
Chain-dependent masks. KDUQ’s iterative rejection of points more than 3σ from the current model, updated during the walk, needs a mask that depends on chain state. Masks are compiled. Workaround: an outer loop of problems with the mask refit between runs. Addition refused by design: it is mutable state inside a spec.
Per-point latent variables. Errors-in-variables in x (Berkson), explicit latent function values on a mesh (Schnabel et al. 2021), or a sampled per-point scale in a scale mixture. Expressible in principle as one Parameter per point, but the dimension defeats the drivers. Addition: latent nodes with analytic marginalisation; large.
Per-point log-likelihood factors for PSIS-LOO and WAIC under a correlated covariance. Workaround: the prefix-difference of recipe 25, one problem per point. Addition: CompiledConstraint.pointwise_log_likelihood from the factor; about 30 lines.
Input-dependent mixture likelihoods. Bayesian model mixing with a per-point two-component Gaussian likelihood (Semposki et al. 2022). Mean mixing is a Model (w(x; θ) f1 + (1 − w) f2); the likelihood form is the first bullet.
Correlation across constraints. By construction. Merge the constraints.
A sampled tempering exponent. weight is a float by type; a sampled η is not a coherent posterior. SafeBayes chooses it outside the sampler (recipe 25).
Nonlinear, non-pointwise comparison maps. space is pointwise so the delta method and the log-Jacobian stay exact. A linear projection is preprocessing (recipe 20); a nonlinear feature of the whole vector would need a full Jacobian. Not planned.
Kronecker-structured covariances. A separable GP over energy and angle on a common grid could factor as a Kronecker product; real data share no grid, so the dense path is used (recipe 23). Addition: a fourth term kind; not planned until a case needs it.
Joint integration over emulator hyperparameters alongside calibration. Bayarri et al. recommend fixing them; not a gap worth filling.

Recipes#

1. Fit a model to data with reported statistical errors#

2. Infer an unknown noise level#

3. Use the reported systematic errors#

4. Infer a normalisation or offset the experiment did not report#

5. Share an error model between datasets, or couple them#

6. One latent scale per dataset#

7. Absorb model deficiency with a Gaussian process#

8. Sample a mean discrepancy explicitly#

9. Heavy tails#

10. Compare in log space#

11. Hold out data and score it#

12. Temper the likelihood#

13. Declare priors#

14. From an EXFOR measurement to a dataset#

15. Evaluate a reaction model on any grid#

16. Drive the calibration with an external sampler#

17. Check the posterior predictive#

18. Compare error models by evidence#

19. Bring my own covariance or term#

20. Preprocess instead of asking for a feature#

21. Get a useful error, not a LinAlgError#

22. Share hyperparameters across datasets, with values that depend on the dataset#

23. A discrepancy correlated across energies and angles#

24. Per-dataset parameters drawn from a sampled hyperprior#

25. SafeBayes: learn the tempering exponent#

26. Unaccounted-for model error per data type (KDUQ)#

27. Peelle’s Pertinent Puzzle: normalise the prediction, never the data#

28. Stacking by leave-one-dataset-out#

29. Cut (modular) posterior by multiple imputation#

30. Leave-one-experiment-out prediction#

31. Simulation-based calibration of the sampler#

32. Emulator as the model, emulator variance as a term#

33. MAP and Laplace approximation#

34. Global error scale factor and unrecognised sources of uncertainty#

35. Energy-dependent parameters and per-comparison model instances#

36. Discrepancy on a physically constrained basis#

37. Correlated systematics between observables of one measurement#

38. The classic normal hierarchical model (eight schools)#

What this API does not express#

21. Get a useful error, not a `LinAlgError`#