CUQIpy Bayesian unfolding: uncertainty-quantified MCMC#

The unfold_cuqi method integrates CUQIpy — Computational Uncertainty Quantification for Inverse Problems (Technical University of Denmark, DTU) — as a full Bayesian back end for neutron spectrum unfolding. Instead of a point estimate, the method draws samples from the joint posterior distribution of the (log-)spectrum and reports the posterior mean, per-bin standard deviations, highest-posterior-density (HPD) credible intervals and MCMC convergence diagnostics (effective sample size, Gelman-Rubin R-hat, acceptance rate).

CUQIpy is an optional dependency: install it with the package extra pip install bssunfold[cuqi] (or pip install cuqipy). When cuqipy is missing, the package still imports normally (the flag bssunfold.CUQI_AVAILABLE is False) and unfold_cuqi raises an informative ImportError.

NumPy 2.4 note. Upstream cuqipy 1.5.1 declares numpy<=2.2.0, which conflicts with this package’s numpy>=2.4.1 floor. A maintained fork with the cap relaxed to numpy<2.5 and the NUTS int()-on-1-element-array issue fixed is available at Radiationsafety/CUQIpy (branch numpy2-support, dist version 1.5.2); it is validated to reproduce upstream results bit-for-bit on identical seeded chains. On NumPy >= 2.4 pre-install it with:

pip install "cuqipy @ git+https://github.com/Radiationsafety/CUQIpy@numpy2-support"

before pip install bssunfold[cuqi] (pip accepts the fork as satisfying cuqipy>=1.5.0). The uv-managed environment (lockfile and CI) resolves cuqipy from the fork automatically via [tool.uv.sources]. On NumPy <= 2.2 the plain PyPI cuqipy works unchanged.

Bayesian model#

Terms. Posterior distribution \(p(\boldsymbol{\theta} \mid b)\) — the probability density of the model parameters (here the log-spectrum) after the measurements are taken; it combines the likelihood \(p(b \mid \theta)\) with the prior \(p(\theta)\) via Bayes’ theorem. Credible interval — an interval containing a stated fraction (e.g. 95 %) of the posterior mass; unlike a frequentist confidence interval it makes a direct probabilistic statement about the spectrum itself (Gelman et al., 2013). HPD (highest posterior density) interval — the shortest credible interval at a given level; every point inside an HPD interval has higher posterior density than any point outside it.

The unknown spectrum is modelled on the log scale,

\[f_j = \exp(\theta_j), \qquad j = 1, \dots, n,\]

which enforces non-negativity by construction and linearises the multiplicative structure of the folded readings. The data model is a Gaussian likelihood on the (whitened, optionally relative-error scaled) readings,

\[b \sim \mathcal{N}\!\big(A\, f(\theta),\; \sigma^2 I\big), \qquad \sigma = \texttt{noise\_level} \cdot \|b\|,\]

and the log-spectrum carries a smoothness prior anchored on a data-driven center \(\mu\) (the non-negative least-squares solution of \(A x = b\), or the user-supplied initial_spectrum mapped to the log scale):

  • prior="gmrf" — a Gaussian Markov random field with a first- or second-order difference precision operator (\(\texttt{gmrf\_order}\) = 1 or 2; higher is smoother);

  • prior="ou" — a dense Ornstein-Uhlenbeck (exponentially decaying) correlation Gaussian prior with the correlation length lengthscale energy bins.

Samplers#

The posterior is explored with one of the CUQIpy samplers selected by sampler:

  • sampler="pcn" — preconditioned Crank-Nicolson (robust random-walk on the log scale; centred parameterisation \(\theta = \mu + t\));

  • sampler="cwmh" — component-wise Metropolis-Hastings (per-bin updates, small memory footprint);

  • sampler="mala" / sampler="ula" — (Metropolis-adjusted) Langevin algorithms using the analytic gradient of the log posterior; the sampler operates on a MAP-whitened parameterisation (Gauss-Newton MAP preconditioning) because the posterior curvature of the unfolding problem spans many orders of magnitude across directions. ula is experimental and may diverge on stiff problems;

  • sampler="nuts" — the No-U-Turn Sampler (native CUQIpy NUTS with dual-averaging step-size adaptation);

  • sampler="gibbs" / sampler="gibbs_nuts" — hierarchical HybridGibbs sampling where the GMRF smoothness precision \(\delta\) is inferred from the data through a conjugate Gamma hyperprior (\(\delta \sim \mathrm{Gamma}(\alpha, \beta)\), exact conjugate update each scan) while the spectral block is updated with pCN or NUTS respectively. The posterior draws of the hyperparameter are returned in cuqi_stats["delta_samples"].

Multi-chain runs (chains >= 2) are seeded independently and the between-chain agreement is quantified with the Gelman-Rubin statistic.

Diagnostics#

The result dictionary contains the usual unfolding fields (spectrum — posterior mean, spectrum_uncertainty — per-bin posterior standard deviation, spectrum_lower / spectrum_upper — HPD bounds at credible_level, residual / residual_norm) plus a cuqi_stats sub-dictionary with the sampling diagnostics:

  • mean, median, std — posterior summaries per energy bin;

  • hpd_lower / hpd_upper — HPD credible interval bounds;

  • ess — effective sample size of the (thinned) chains;

  • rhat — Gelman-Rubin potential scale reduction factor;

  • acc_rate — sampler acceptance rate(s);

  • theta_samples / samples — raw posterior draws;

  • delta_samples — posterior draws of the smoothness precision (hierarchical samplers only);

  • prior_center — the NNLS-anchored prior center.

Interpreting the posterior (CUQIpy interpretability)#

A point estimate from a classical solver answers “what spectrum fits the data?”; the CUQIpy posterior answers the interpretability questions a point estimate cannot:

  • How uncertain is each energy bin? — the per-bin posterior standard deviation and the HPD interval width. Bins whose HPD interval spans decades are not resolved by the sphere set; narrow intervals flag data-determined regions. This is the Bayesian counterpart of the local detector-sensitivity analysis of Interpreting Unfolding Results, but it is global: it marginalises over all parameter directions at once instead of perturbing one reading at a time (Tarantola, 2005).

  • Are the reported uncertainties trustworthy? — only if the MCMC chains have converged and are long enough. Check rhat and ess (below) before quoting any credible interval.

  • Which smoothness prior does the data support? — with the hierarchical samplers the posterior of the smoothness precision delta_samples shows which smoothness levels the measurements themselves favour, a data-driven regularisation choice rather than a hand-tuned \(\alpha\) (Gelman et al., 2013).

  • Is the non-negativity constraint distorting the result? — the log-scale parameterisation \(f = \exp(\theta)\) enforces positivity exactly, so the posterior mass can never cross zero; a posterior concentrated far from zero means the constraint is inactive, while a posterior piling up at the lower edge flags the same “artificial corner” that the shadow-price analysis of Interpreting Unfolding Results detects for QP solvers.

Diagnostics and how to read them#

  • Gelman-Rubin :math:`hat R` (cuqi_stats["rhat"]) compares the between-chain variance \(B\) with the within-chain variance \(W\) of \(m\) chains of length \(n\):

    \[\hat V = \frac{n-1}{n}\,W + \frac{B}{n}, \qquad \hat R = \sqrt{\hat V / W}.\]

    \(\hat R \approx 1\) means the chains agree — they sample the same posterior. \(\hat R > 1.1\) is a hard warning: the credible intervals are unreliable and the run needs more samples or a better tuned sampler (Gelman and Rubin, 1992). The rank-normalised split-\(\hat R\) computed by ArviZ, used here, is the modern standard with the stricter practical threshold \(\hat R < 1.01\) (Vehtari et al., 2021).

  • Effective sample size (ESS) (cuqi_stats["ess"]) converts \(N\) correlated MCMC draws into the number of independent draws that would carry the same information, \(\mathrm{ESS} = N / (1 + 2\sum_t \rho_t)\) with \(\rho_t\) the autocorrelation at lag \(t\) (Vehtari et al., 2021). An ESS below a few hundred means the reported posterior mean is itself noisy — increase n_samples or thin. ESS also quantifies why one sampler beats another: NUTS typically achieves a much higher ESS per draw than random-walk pCN on the strongly correlated unfolding posterior, at a higher cost per step (Hoffman and Gelman, 2014).

  • Acceptance rate (cuqi_stats["acc_rate"]) — the fraction of proposals accepted. For random-walk samplers the optimal regime is roughly 0.2-0.5 (Roberts et al., 1997); near 1 means tiny steps (slow exploration, low ESS), near 0 means the chain is stuck. NUTS adapts its step size to a target acceptance automatically (Hoffman and Gelman, 2014).

  • HPD width vs. bin energy — plot result["spectrum_upper"] - result["spectrum_lower"] on a log scale: a narrow plateau means the spectrum shape is identified where the sphere responses overlap; flaring wings mark the thermal and high-energy edges where the Bonner sphere set loses resolution (Thomas and Alevra, 2002).

Relation to the pyoptexplain interpretation#

Interpreting Unfolding Results (pyoptexplain) explains the optimisation problem — which constraints bind, what each detector is worth — by analysing the solved QP. CUQIpy explains the inference problem — how much of the spectrum is actually determined by the data. They are complementary views of the same underdetermined system:

Question

pyoptexplain (Interpreting Unfolding Results)

CUQIpy (unfold_cuqi)

Is the solution stable?

±1-5 % reading perturbation tests

Posterior std / HPD width per bin

Which detectors matter?

One-at-a-time detector importance

Global marginalisation over all directions

Is \(x \ge 0\) distorting?

Shadow prices / non-negativity relaxation

Posterior mass position near zero

Is the regularisation right?

\(\alpha\) sweep of point solutions

Hierarchical posterior of the smoothness precision

Convergence certificate

Solver status, KKT residuals

\(\hat R < 1.01\), ESS, acceptance rate

Example: quantify which bins are data-determined#

from bssunfold import Detector

det = Detector()
result = det.unfold_cuqi(
    readings,
    sampler="gibbs_nuts",
    prior="gmrf",
    hierarchical=True,
    n_samples=2000,
    n_burnin=1000,
    chains=2,
    credible_level=95.0,
    random_state=42,
)
stats = result["cuqi_stats"]

# 1) Convergence gate: only trust intervals after this passes
assert stats["rhat"] < 1.01 and stats["ess"] > 200

# 2) Per-bin relative uncertainty = interpretability map
rel_unc = result["spectrum_uncertainty"] / result["spectrum"]
# bins with rel_unc << 1 are data-determined;
# bins with rel_unc ~ 1 are prior-dominated.

# 3) HPD width per bin (log-scale view)
hpd_width = result["spectrum_upper"] - result["spectrum_lower"]

Example#

from bssunfold import Detector

det = Detector()                # default GSF response functions
result = det.unfold_cuqi(
    readings,                   # measured sphere readings
    sampler="gibbs_nuts",       # hierarchical Gibbs with a NUTS block
    prior="gmrf",               # GMRF log-spectrum prior, order 1
    gmrf_order=1,
    hierarchical=True,          # infer the smoothness precision
    n_samples=2000,
    n_burnin=1000,
    thin=1,
    chains=2,                   # >= 2 chains enable R-hat
    credible_level=95.0,
    random_state=42,
)

spec = result["spectrum"]       # posterior mean spectrum
lo, hi = result["spectrum_lower"], result["spectrum_upper"]
stats = result["cuqi_stats"]
print(stats["ess"], stats["rhat"], stats["acc_rate"])

A lower-level entry point solve_cuqi_bayesian() exposes the same model on the raw response-matrix/response-vector level for use outside the Detector workflow.

References#

  • CUQIpy: CUQIpy — Computational Uncertainty Quantification for Inverse Problems, DTU, CUQI-DTU/CUQIpy

  • Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., Rubin, D. B. (2013). Bayesian Data Analysis, 3rd ed. CRC Press. doi:10.1201/b16018

  • Gelman, A. & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences, Statistical Science 7(4), 457-472. doi:10.1214/ss/1177011136

  • Hoffman, M. D., Gelman, A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Machine Learning Research 15, 1593-1623.

  • Riis, N. A. B. et al. (2019). pCN sampling for Bayesian inverse problems (preconditioned Crank-Nicolson in CUQI/CUQIpy).

  • Roberts, G. O., Gelman, A., Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. Ann. Appl. Probab. 7, 110-120. doi:10.1214/aoap/1034625254

  • Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM. doi:10.1137/1.9780898717791

  • Thomas, D. J., Alevra, A. V. (2002). Bonner sphere spectrometers — a critical review. Nucl. Instrum. Meth. A 476, 12-20. doi:10.1016/S0168-9002(01)01379-1

  • Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved \(\hat R\) for assessing convergence of MCMC. Bayesian Analysis 16, 667-718. doi:10.1214/20-BA1221