Reinvented everywhere, packaged nowhere: pastry’s opening

No maintained, general-purpose Python package takes N heteroscedastic per-object likelihoods over K bins, or posterior samples drawn under an interim prior, and returns population fractions with calibrated uncertainty plus smoothness, selection and interim-prior extensions, even though at least ten fields have independently built exactly this model. Astronomers call it hierarchical n(z), colour–magnitude or occurrence-rate inference; particle physicists call its point estimate iterative Bayesian unfolding; statisticians call it the NPMLE or g-modelling; machine-learning researchers call it EMQ or Bayesian quantification; geneticists call it mixed-stock analysis. The resulting code is domain-locked (gwpopulation, RAIL, QuaPy, RooUnfold), point-estimate-only (mixsqp, REBayes, extreme deconvolution), confined to R or Julia (ashr, deconvolveR, Empirikos.jl, rubias, StarFormationHistories.jl), or paper code that was never packaged or is no longer maintained (CHIPPR, Starlight, exopop, PyFBU and every neural deconvolution repository surveyed). The one maintained Python implementation of the exact baseline in astronomy, RAIL’s VarInfStackSummarizer, is mean-field variational Bayes that the HSC Y3 team says under-covers by about 20% (Rau et al. 2023). Pastry’s opportunity is therefore engineering and statistical hygiene rather than a new model: explicit input types for likelihoods, interim priors and selection; exact posterior sampling by default, with EM and variational modes labelled as approximations; overdispersion and smoothness priors, because a single Dirichlet propagates only shot noise; and simulation-based calibration shipped as a user-facing tool. The evidence favours an emcee/scikit-learn-style layered API that returns plain arrays and ArviZ-compatible draws, a NumPy/SciPy core with JAX and flows as optional extras, compression rather than brute force at N ≥ 10⁸, and a PyPI name other than pastry, which is taken. The main unknowns are empirical: nobody has published a coverage comparison of Gibbs, NUTS, variational and EM-plus-Laplace inference for this model as per-object likelihoods broaden, and no neural deconvolution method has been coverage-tested, so pastry will have to generate that evidence itself.

At least ten fields reinvented the same mixture-weights model

Strip away the vocabulary and every tool in this survey targets one likelihood, p({d} | ρ) = Π_i Σ_k L_ik ρ_k, where L_ik is the probability of object i’s data if its true value lies in bin k and ρ is the population histogram on the simplex. Hogg, Myers & Bovy set out the general recipe in 2010: take per-object likelihood functions, or posterior samples drawn under a known uninformative prior, and treat the population fit as “a kind of heteroscedastic deconvolution”, with photometric redshifts named as an application (Hogg, Myers & Bovy 2010). On fixed bins, maximising this likelihood is a finite-dimensional convex problem for any per-object likelihoods, the observation of Koenker & Mizera on which modern NPMLE solvers rest (Soloff, Guntuboyina & Sen), and each evaluation costs one N×K matrix–vector product. Adding a Dirichlet prior on ρ and latent bin labels z_i gives the hierarchical Dirichlet–multinomial model that pastry implements.

Astronomy alone has built it at least four times. Leistedt, Mortlock & Peiris put a Dirichlet prior on photo-z histogram heights and sampled it with a two-step Gibbs sampler on binned likelihoods, drawing 10⁴ samples in about ten minutes on a laptop (Leistedt et al. 2016). Leistedt & Hogg ran the same sampler over 1.4 million Gaia TGAS stars, a model with “4 million parameters in total, most of which are marginalized out”, by precomputing each star’s likelihood in each fixed colour–magnitude kernel (Starlight). Foreman-Mackey, Hogg & Morton fitted piecewise-constant exoplanet occurrence rates with a completeness function, interim-prior importance weights and a Gaussian-process smoothness prior (exopop ms.tex). The LVK’s 2016 BinnedCBCMasses Stan program takes per-event weight matrices, which are L_ik by another name, under a GP prior on log-rates with marginalised selection uncertainty (population.stan), and PixelPop and gppop now run that structure on 10⁴ to 8×10⁴ bins on GPUs (PixelPop; gppop).

Outside astronomy the same mathematics appears under other names. Iterative Bayesian unfolding in particle physics is the EM algorithm for these weights, identical to Richardson–Lucy and Shepp–Vardi (Kuusela & Panaretos), and the OmniFold authors prove that their unbinned per-event iteration is a generalised EM with the true values as latent variables (Andreassen et al.); binned unfolding is pastry’s input with one shared row per reconstructed bin, L_ik = R_b(i),k. In machine-learning quantification, QuaPy’s EMQ code divides classifier posteriors by the training prevalence and averages them, which is EM with interim-prior correction (QuaPy aggregative.py), and its BayesianMAPLS samples the same likelihood under a Dirichlet prior with NUTS (QuaPy _bayesian.py). Fisheries geneticists estimate stock proportions by EM (ONCOR, SPAM) or Dirichlet Gibbs (BAYES, rubias) from per-fish likelihoods under each reference population (ONCOR manual; rubias infer_mixture.R), microbial source tracking runs the Dirichlet Gibbs sampler on count tables (SourceTracker2), and RNA-seq quantifiers such as salmon run the same EM, plus a variational variant, over equivalence classes of reads (salmon). A 2023 comparison of ordinal quantifiers found algorithms “proposed by authors from very different research fields, such as data mining and astrophysics, who were unaware of each others’ developments” (Bunse et al.).

One identity ties the point estimators together and gives pastry a teaching device. From the optimality conditions on the simplex, the maximum-likelihood ρ satisfies ρ_k = (1/N) Σ_i r_ik, where r_ik ∝ L_ik ρ_k is object i’s posterior over bins under the population ρ; the NPMLE is the fixed point of “stack the per-object posteriors computed under the current population” (Soloff et al.; mixsqp manual). By direct algebra, stacking posteriors computed under an interim prior is exactly one EM step from that prior, IBU with early stopping is a truncated EM, and a Dir(α) MAP adds α_k − 1 pseudo-counts. That explains why stacking is biased unless data or prior are “perfectly informative” (Malz 2021) and why HEP treats the iteration count as a hidden regulariser (Cowan 2002). It also hands RAIL and QuaPy users a one-sentence migration story: pastry runs to convergence, with uncertainty, the iteration they already perform once.

Pastry concept	Astronomy	GW populations	HEP unfolding	EB statistics	ML quantification	Genetics and genomics
Per-object likelihood L_ik	binned photo-z likelihood; per-star kernel likelihood	per-event weights: PE samples ÷ PE prior	response-matrix row of the reconstructed bin	likelihood matrix L	classifier posterior ÷ training prevalence	genotype likelihood per baseline stock; read compatibility
Point estimate of ρ	stacking (one EM step); EM	maximum hyper-likelihood	IBU, D’Agostini, Richardson–Lucy	NPMLE (Kiefer–Wolfowitz)	EMQ, SLD, MLLS	conditional ML (SPAM, ONCOR); EM/VBEM (salmon)
Posterior over ρ	Dirichlet Gibbs; mean-field VB (RAIL)	NUTS over binned log-rates	Fully Bayesian Unfolding	mixVBEM (ashr, variational)	BayesianMAPLS, BayesianKDEy	BAYES, rubias, SourceTracker Gibbs
Interim-prior correction	CHIPPR; importance weights (FM+14)	prior column divided out	“simulation-prior dependence”	none found	division by π_train	none found
Selection	completeness; Poisson process	injections, P_det	efficiency and acceptance	truncation (deconvolveR tg)	none found	none found
Smoothness	logistic GP, GP, random walk	ICAR, GP, P-splines	Tikhonov, GP unfolding	log-spline g (Efron)	ordinal Tikhonov	none found

The convergence has a practical payoff: each field’s benchmarks and validation designs transfer directly. HEP contributes coverage studies of regularised estimators (Kuusela & Panaretos), empirical Bayes contributes canonical datasets bundled with Empirikos.jl and REBayes (Empirikos datasets), quantification contributes the LeQua challenges, including an ordinal task (LeQua 2024), and Leistedt & Hogg’s TGAS run is a ready-made astronomy benchmark (Starlight).

Maintained tools stop at point estimates, one domain, or R

Across the census the pattern is consistent. Maintained tools restrict their inputs to one field’s data products, stop at point estimates, or live in R and Julia; the tools that implement the full hierarchical model with honest uncertainty are paper code. The subsections below compare the fields, and a census table closes the section.

Astronomy’s exact samplers are stale and its maintained tool is approximate

For photometric redshifts, LSST DESC’s RAIL is the only maintained, pip-installable, general tool that turns per-galaxy p(z) into n(z) samples (RAIL estimation docs). Its only hierarchical option, VarInfStackSummarizer, is coordinate-ascent mean-field VB for a flat Dir(1) prior that evaluates each PDF at grid points rather than integrating over bins, applies no interim-prior correction and holds all data in memory (var_inf.py). Its originators reported that the scheme “yielded underestimated error bars, likely due to the restricted covariance structure of the Dirichlet” (Rau et al. 2022), and HSC Y3 put the under-coverage at about 20%, from “numerical experience” rather than a documented test (Rau et al. 2023). RAIL’s own example notebook warns that it “tends to get a little wonky if you use too many bins” and passes an niter argument that the source never reads, because the configuration key is n_iter (Sampled_Summarizers.ipynb); the package’s summarizer tests are end-to-end smoke tests with no accuracy or coverage assertions (rail_base tests). Rubin’s first data-preview photo-z paper still reports stacked n(z) without an uncertainty band (Zhang et al. 2025).

The exact hierarchical samplers exist but are frozen. CHIPPR, the only photo-z code that divides out the interim prior, calls itself “stale and should be understood to be a prototype” and ships a Python 2.7 demo (aimalz/chippr); Leistedt’s Gibbs sampler in Photoz-tools was last pushed in 2020 (ixkael/Photoz-tools), BHM-nz in 2023 (Harry45/BHM-nz) and frankenz in 2018 (joshspeagle/frankenz). A 2026 KiDS application of the Leistedt sampler still subsampled to 10⁵ galaxies because running on all of them “would be too slow” (Kyriacou et al. 2026). The survey pipelines (DES SOMPZ with three-step Dirichlet resampling, KiDS SOM calibration, HSC’s logistic-GP model) obtain more realistic errors, but they work from colour cells, calibration samples and clustering rather than generic likelihood matrices (Sánchez et al. 2020; Myles et al. 2021). Downstream, the interfaces are standardising on qp ensembles, CosmoSIS realisation files for hyperrank, and nz_prior’s compressed shift, PCA and GP priors (qp; nz_multirank; nz_prior), which fixes pastry’s export targets.

Stellar populations tell the same story. Leistedt & Hogg’s colour–magnitude model is pastry’s baseline almost exactly, but its Starlight code was last committed in 2017 and never released (Starlight). The extreme-deconvolution (XD) family that dominates practice assumes Gaussian per-object noise and returns a maximum-likelihood mixture with no posterior uncertainty. Bovy’s reference C code is maintained, with CI on Python 3.9–3.14, but is not on PyPI (extreme-deconvolution); astroML’s XDGMM lacks projection matrices (xdeconv.py); and only pygmmis, last released in 2022, handles truncation, and only with a known completeness function (pygmmis). Anderson et al. note that using XD as a prior “is technically using the data twice” (photoParallax paper.tex). The best-engineered code for bin weights with uncertainties is in Julia: StarFormationHistories.jl fits Hess-diagram counts with a Poisson likelihood ratio, offers MAP with inverse-Hessian errors plus HMC and MCMC, ships CI, coverage reports and documentation, and in v2.0.0 (October 2026) added per-sample derived quantities that “account for the correlations between time bins” and background templates that stop contaminants biasing fits (CHANGELOG). PANCAKE, by contrast, fits Poisson Hess diagrams with a mean-squared-error loss and clipped gradient descent (Zheng et al. 2025).

Exoplanet occurrence studies produced pastry’s entire extension set a decade ago, but only as paper code. Foreman-Mackey, Hogg & Morton combined an inhomogeneous Poisson process with completeness, importance-sampled posteriors divided by the interim prior, a piecewise-constant rate density and a GP smoothness prior with marginalised hyperparameters, sampled by elliptical slice sampling (exopop ms.tex); the repository was “written specifically for our paper” and last touched in 2014 (dfm/exopop). The Kepler DR25 occurrence code is “not intended to be used by the community” (DR25-occurrence-public), and a 2026 hierarchical-histogram script silently replaces log(0) by zero, dropping objects with no support (sblunt/eccentricities). The reusable astronomy packages are narrow. PosteriorStacker fits a flat-Dirichlet histogram to one-dimensional posterior samples drawn under a flat prior, through a command-line script, with inconsistent GPL/AGPL licensing (PosteriorStacker). MRExo builds pastry’s N×K matrix internally but returns simplex MLEs with bootstrap errors and has no tests (mrexo). dftools, the most general luminosity-function tool, with arbitrary per-object PDFs and selection, is R code last committed in 2021, and its Python port never implemented non-parametric models (dftools; pydftools).

Gravitational-wave codes are mature, GPU-ready and GW-shaped

Gravitational-wave population inference is the best-engineered neighbour. gwpopulation (MIT, JOSS 2025, NumPy/CuPy/JAX backends) is the LVK reference implementation of the importance-reweighted hierarchical likelihood: it divides a prior column out of each event’s posterior samples, estimates detection efficiency from injections, and returns per-sample Monte Carlo variances (hyperpe.py; vt.py). PixelPop runs NUTS over 10⁴ ICAR-coupled bins on GPUs (PixelPop), and gppop runs a Kronecker-factorised GP over about 8×10⁴ bins on four GPUs, taking a per-bin weight matrix, i.e. L_ik, as input (gppop). The whole ecosystem is shaped by GW data, however. Inputs are posterior-sample tables plus injection sets with total_generated and analysis_time; gwpopulation is coupled to bilby and wcosmo and switches array backends through global state (backend.py); GWInferno, icarogw and gppop are not on PyPI; and test suites range from about a dozen modules in gwpopulation to 54 lines in FIGARO and none in gppop (gwpopulation; FIGARO; gppop). The generic bilby.hyper likelihood is domain-neutral but has no selection term and no Monte Carlo variance (bilby 2.8.2). What transfers is gwpopulation’s plain-function model contract, the popsummary result schema of hyperposterior samples, reweighted per-object samples and rates on grids (popsummary), and the Monte Carlo accounting discussed below.

HEP unfolding has lost its lightweight Bayesian package

Binned unfolding is pastry’s model with a single shared likelihood row per reconstructed bin, and production still runs on C++. RooUnfold (BSD-3, v3.1.0, on conda-forge but not PyPI) bundles IBU, SVD, Tikhonov, GP and TUnfold methods (RooUnfold), and TUnfold’s last standalone release dates from 2019 (TUnfold). The Python tools froze. PyUnfold’s last release was in 2018, and it returns per-bin errors rather than a covariance matrix (unfold.py). PyFBU, the one package that sampled a posterior with nuisance parameters, is “no longer actively maintained” and needs Theano, PyMC3 3.5 and numpy<1.20 (pyFBU/fbu), so the 2025 Neural Posterior Unfolding authors re-implemented Fully Bayesian Unfolding in TensorFlow notebooks to “provide an updated version” (Torales Acosta et al.). Since 2020, method development has moved to OmniFold and more than fifteen generative or discriminative unbinned variants, nearly all single-paper code; the 2025 practitioners’ guide reports ensembles of 4–10 networks (100 for ATLAS), 50–100 bootstrap runs and 500 to 10,000 GPU hours per analysis (Canelli et al.). No HEP tool accepts a heteroscedastic per-object likelihood matrix: binned tools assume one response matrix, and unbinned tools assume paired simulation (Andreassen et al.).

Statistics owns the solvers and the interval estimates, outside Python

Empirical-Bayes deconvolution has the most mature algorithms, and almost all of them live in R or Julia. mixsqp solves the mixture-weights MLE for any N×K likelihood matrix, accepting log input, and achieved at least a 100-fold speed-up over the MOSEK interior-point solver at n ≈ 10⁶ (Kim et al. 2020). REBayes requires the commercial MOSEK solver (REBayes). ashr accepts any likelihood matrix with a Dirichlet prior and, through mixVBEM, is the only EB package surveyed that treats ρ as random, by mean-field VB (ashr manual). deconvolveR returns delta-method standard errors and bias for a penalised log-spline g (deconv.R), and Empirikos.jl offers F-localisation and AMARI confidence intervals for functionals of g (Empirikos.jl). Python’s entire offering is npeb, an alpha last released in 2024 with no tests, print-based progress and an n×n kernel by default (GLMixture.py), and mixsqpx, an alpha JAX port restricted to Python 3.12 (mixsqpx). Python researchers fall back on rpy2, R and MOSEK, as the replication code for Chen’s Econometrica paper does (close-replication). ebnm’s authors name the cost of this fragmentation: “each existing package has a different interface and outputs, which hinders comparisons” (Willwerscheid et al.).

Quantification and genetics ship the Bayesian model for classifiers and genotypes

QuaPy (BSD-3, v0.2.3, October 2026) has the most complete uncertainty toolkit built around this model in Python. Its BayesianMAPLS samples Σ_i logsumexp_k(log p(k|x_i) + log π_test,k − log π_train,k) under a Dirichlet prior with NumPyro NUTS, every confidence-aware method returns an (estimate, region) pair from predict_conf, and calibrate_temperature widens over-tight posteriors by choosing a likelihood temperature against validation coverage (QuaPy _bayesian.py; QuaPy confidence.py). Its inputs, however, are features plus a scikit-learn classifier; no documented entry point takes a likelihood matrix, and none adds smoothness or selection (QuaPy). In fisheries genetics, rubias runs Dirichlet–multinomial Gibbs with parametric-bootstrap bias correction and baseline resampling, and returns traces, per-fish posteriors and posterior-predictive totals (rubias infer_mixture.R); CRAN archived it in February 2026 because the maintainer’s email was undeliverable (CRAN rubias), and a GitHub search found no Python stock-identification tool at all (GitHub search). The only Python member of this family, SourceTracker2, has not been released to PyPI since 2016 (PyPI sourcetracker).

Neural deconvolution exists only as single-paper code

Normalising-flow and diffusion deconvolution is a research frontier without a library. Dockhorn et al.’s flow-deconvolution code has no licence and a single 2020 commit (density-deconvolution); Neural Empirical Bayes is MIT-licensed notebooks last touched in 2022 (NEB); Yan & Sanders’s 2025 Gaia flows have neither packaging nor tests (Denoising-Milky-Way-with-NF); Neural-g, which trains a softmax-output network on exactly pastry’s gridded objective, is an R wrapper around PyTorch without tests (neuralG); and DiEM’s diffusion-prior EM needs a Slurm cluster (diem). None of them reports calibrated uncertainty on the learned density (Vandegar et al.; Yan & Sanders 2025). The infrastructure underneath is healthy: flowjax (MIT, v19.1.2, October 2026) accepts arbitrary training losses and interoperates with NumPyro, zuko serves PyTorch, and NumPyro 0.22 supplies NUTS (flowjax; zuko; NumPyro).

Census at a glance

The table compares the tools closest to pastry’s scope on the attributes that matter for the specification. “Generic input” means the tool accepts an arbitrary per-object likelihood matrix or posterior samples with a declared interim prior.

Tool (field)	Stack, licence, status (Oct 2026)	Per-object input	Engine and population uncertainty	Selection · smoothness · interim prior	Role for pastry
RAIL VarInfStack (photo-z)	Python, MIT; pz-rail-base 2.0.11, Sept 2026	qp PDFs evaluated at grid points	mean-field CAVI; Dirichlet draws that under-cover	none · none · none	drop-in replacement target; RAIL plug-in
Photoz-tools / BHM-nz (photo-z)	Python notebooks, MIT; last push 2020 / 2023	binned template likelihoods	exact two-block Gibbs; posterior samples	selection variable · generalised Dirichlet proposed, not built · n/a	reference implementation of the baseline
chippr (photo-z)	Python 2, MIT; “stale” prototype, 2022	gridded posteriors + interim prior	emcee; posterior samples	none · squared-exponential hyperprior · yes	only photo-z interim-prior code
Starlight (Gaia CMD)	Python/Cython, MIT; last commit 2017	precomputed per-kernel likelihoods	Gibbs over 1.4 M stars; posterior samples	none · none · n/a	TGAS benchmark
extreme-deconvolution, pygmmis (astro densities)	C + Python; XD maintained, not on PyPI; pygmmis 1.2.3, 2022	Gaussian mean and covariance only	EM maximum likelihood; no uncertainty	pygmmis: known completeness · covariance regulariser · no	Gaussian-noise baseline to beat
StarFormationHistories.jl (Hess diagrams)	Julia, MIT; v2.0.0, Oct 2026; CI, codecov	binned counts + template matrix	Poisson likelihood ratio; MAP + Hessian, HMC, MCMC	via templates · hierarchical metallicity relations · n/a	engineering exemplar; response-matrix input
exopop (occurrence rates)	Python, MIT; 2014 paper code	posterior samples	Poisson process, elliptical slice; samples	completeness · GP · yes	the full extension set, in paper form
PosteriorStacker (astronomy)	Python, GPL/AGPL (inconsistent); 0.6.4, Sept 2026	1-D samples under a flat prior	flat-Dirichlet histogram via UltraNest; samples + evidence	none · none · flat prior only	closest PyPI analogue; too narrow
MRExo (exoplanets)	Python, GPL-3; 1.1.5, 2023; no tests	values + asymmetric errors, mapped to N×K	simplex MLE; bootstrap	none · basis degree by CV/AIC · n/a	same internal structure
gwpopulation (GW)	Python, MIT; 1.3.1, Dec 2025; NumPy/CuPy/JAX	samples with prior column + injections	any bilby sampler; samples + Monte Carlo variance	injections · spline nodes · yes	likelihood contract to copy
PixelPop / gppop (GW)	JAX/NumPyro, MIT / PyTorch + Pyro, GPL; PixelPop 0.3.4, Oct 2026; gppop not on PyPI, no tests	samples mapped to per-bin weights (= L_ik)	GPU NUTS over 10⁴ (ICAR) to 8×10⁴ (Kronecker GP) bins	injections · ICAR / GP · yes	proof that smooth binned priors scale
RooUnfold (HEP)	C++/ROOT, BSD-3; 3.1.0, conda-forge only	one shared response matrix	IBU, SVD, Tikhonov, GP; propagated errors, toys	efficiency · Tikhonov · simulation prior	HEP benchmark; response-matrix adapter
PyUnfold / PyFBU (HEP)	Python; frozen 2018 / unmaintained (Theano)	response matrix	IBU with per-bin errors / NUTS posterior	efficiency · spline or Tikhonov · n/a	vacant niche pastry can fill
mixsqp, REBayes (EB)	R/C++; mixsqp 2023; REBayes 2026, needs MOSEK	any N×K likelihood matrix	convex NPMLE; no uncertainty	none · post-hoc kernel smoothing only · no	golden-value references
ashr (EB)	R; CRAN 2.2-63	any matrix + Dirichlet prior	MAP, or mean-field VB (mixVBEM)	none · unimodal dictionary · no	only EB tool treating ρ as random
deconvolveR, Empirikos.jl (EB)	R / Julia; both active in 2026	P/Q/y matrices; likelihood types	log-spline + delta-method errors; NPMLE + F-localisation and AMARI intervals	truncation · log-spline / convex prior classes · no	the only interval estimates for g
npeb, mixsqpx (EB)	Python; alpha, 2024 / 2026 (Python 3.12 only)	any kernel matrix	cvxpy/MOSEK; JAX mix-SQP; no uncertainty	none	the entire Python EB offering
QuaPy (ML quantification)	Python, BSD-3; 0.2.3, Oct 2026; 16 test modules	classifier outputs, not L	EM; NUTS; bootstrap regions; temperature calibration	none · none · ÷ training prior	uncertainty API to copy
rubias (fisheries genetics)	R/C++, CC0; CRAN-archived Feb 2026	genotypes vs reference baselines	Gibbs, bias-corrected bootstrap, baseline resampling; traces	none · none · n/a	prior and aggregation lessons
salmon / kallisto (RNA-seq)	Rust / C++, BSD	reads in equivalence classes	EM/VBEM; bootstrap or Gibbs replicates	none	scaling by deduplication
density-deconvolution, NEB, Denoising-MW-NF, neuralG (neural)	PyTorch research code, 2020–2026; several unlicensed; no tests	Gaussian noise, Monte Carlo samples or grids	IWAE or Monte Carlo marginal likelihood; point estimates	selection only by masking Monte Carlo draws · none · implicit flat prior	reimplement; do not depend

frankenz’s Gibbs sampler is exact; its pairwise sampler and flexible reference are not

frankenz holds pastry’s direct ancestors: a collapsed “population” sampler and a data-augmentation “hierarchical” sampler for n(z) (last push October 2018; one PyPI release, 0.1.5; joshspeagle/frankenz). The audit read samplers.py and tutorials 5 and 6 at commit b753e64, ran the code unmodified on NumPy 2.5 and SciPy 1.18, and tested it on freshly built mocks, because the tutorials’ pickled mock is not in the repository. The references were exact quadrature for three bins, Geweke “getting it right” tests, and ablations of a re-implementation that reproduces frankenz bit for bit. The full write-up, with every table and script, is in the research notes (frankenz_audit.md).

The hierarchical sampler is correct. Its object step and conjugate Dirichlet step form an exact two-block Gibbs sampler. Against the exact posteriors it gives an rms z-score of 0.94 over 132 statistics, it matches an independent reference at K=15, N=2000 within |z| ≤ 2.0, and it passes the Geweke test. The two samplers provably target the same marginal posterior for the bin fractions, so any difference between tutorials 5 and 6 comes from the population sampler alone.

The population sampler does not leave the posterior invariant. It makes Metropolis moves along pairs of bins with a proposal width that depends on the current state, but its acceptance test omits the Hastings factor that such a width requires, and it holds the width fixed across three inner steps. Either defect alone breaks the sampler; restoring both recovers the exact posterior. The consequences are large:

• any bin whose posterior reaches zero is driven to about 1e-13 within 2–15 saved samples at the tutorial settings, with errors up to z ≈ 200, and its mass moves to the other bins (up to 1.2 posterior SD);

• with no bin near zero, every posterior SD is still 8–18% too large;

• in the Geweke test the chain collapses onto a vertex of the simplex.

Sparse bins are exactly where this bites: distribution tails, edge bins and bins under broad likelihoods. It fits tutorial 5’s remark that the result is “a little off, especially near zero”, which the notebook attributed to the prior.

The “flexible reference” step belongs to no joint model. It is not the full conditional of any distribution in which the fractions given the counts are Dirichlet: a compatibility check that must be 0 gives 5.33. It fails the Geweke test with max |z| = 70, against 1.7 for an exact replacement step. In practice it discards an uncontrolled share of even a representative reference sample, inflating SDs by a factor of 1.1–1.9, and the share is set by the ratio of catalogue to reference size rather than by any stated trust in the reference. A coherent replacement is either a power prior, Dir(α + λm) with λ in [0, 1], or a two-level model, ρ ~ Dir(c ρ_ref) with ρ_ref ~ Dir(β + m) and c fixed or inferred.

Wrong inputs cause the largest errors. Both samplers need likelihoods, but the tutorials and docstrings call the inputs “PDFs”.

• Posteriors computed under a non-uniform interim prior bias the fractions by up to 12.6 posterior SD and cut 68%/95% interval coverage to 0.29/0.52. Dividing out the interim prior removes the bias.

• Tutorial 5’s normalisation over objects instead of bins (cell 25) would bias results by up to 50 SD if implemented as written.

• frankenz’s bin integral, used where the model needs the bin average, biases results by up to 28 SD once bins have unequal widths.

Software defects compound the statistical ones.

• run_mcmc never resumes from the last sample: both samplers restart from the stacked PDFs.

• reset() leaves old samples in place.

• The population sampler crashes with a numpy.random.Generator.

• Bins with zero stacked mass, 37 of 80 in the tutorial set-up, are frozen at exactly zero.

• No input is validated.

Mixing is the remaining limit. With a flat prior, the Gibbs sampler’s median autocorrelation time grows from 7 to 74, 162 and 180 sweeps as the likelihood width grows from 0.5 to 1, 2 and 4 bin widths. A simple HMC in log-ratio coordinates stays at 1–14, giving 1.4–10 times the per-bin effective samples per second of a vectorised Gibbs sampler and about 200 times that of the original code at tutorial scale. The tutorials’ hierarchical runs carry only about 120 effective samples per bin. Vectorising frankenz’s per-object Python loop alone, which takes 99% of each sweep, speeds a sweep up 6–9 times.

Audit finding	pastry v0.1 response
Pairwise sampler invalid	Not ported. NUTS on isometric log-ratio coordinates is the default; every engine is tested against exact K=3 posteriors.
Gibbs sampler exact but slow to mix for broad likelihoods	Kept, vectorised, with Rao-Blackwellised per-object outputs; R-hat, ESS and MCSE are reported by default.
Flexible reference incoherent	Not ported. Dirichlet.from_counts(m, strength=λ) gives the power prior, and a logistic-normal prior centred on the log reference histogram lets the data decide how far to depart from it; a hierarchical Dirichlet with inferred concentration (spec §5.4) is planned.
Posteriors used as likelihoods	Likelihoods.from_posteriors requires the interim prior and divides it out; a slow regression test shows the bias of skipping this.
Bin integral vs bin average	Likelihoods stores bin averages and takes an explicit representation; unequal-bin tests cover both.
Frozen zero-support bins, exact-zero Dirichlet draws	Interior initialisation; log-space Gamma draws; a test checks the prior-driven mass α_k/(A + N).
Resume, reset, RNG, validation	Resumable engines (run(a); run(b) equals run(a + b)), Generator-only randomness with spawned chain streams, strict input checks.
Geweke tests need power	The suite includes a deliberately broken kernel that the Geweke test must reject.

Honest error bars fail in six documented ways

“Realistic uncertainties” is where the landscape is weakest. The engines that ship differ less in their point estimates, most of which converge to the same maximum, than in what their error bars mean, and the literature documents six distinct ways those error bars go wrong. The table summarises the engines; the paragraphs that follow give the evidence.

Engine	Ships in	What its uncertainty means	Documented weakness
EM, IBU, EMQ	RooUnfold, PyUnfold, QuaPy, ONCOR, salmon	nothing natively; propagated errors or bootstrap added on	slow near the optimum; a zeroed bin never recovers; the iteration count acts as an unstated regulariser
Convex NPMLE solvers	mixsqp, REBayes, npeb	none; a sparse, discrete ĝ	intervals collapse onto support points; MOSEK dependence
Mean-field VB	RAIL VarInfStack, ashr mixVBEM, salmon VBEM	a Dirichlet approximation to p(ρ)	under-covers, worst when per-object likelihoods overlap
Bootstrap or replicas	OmniFold, QUnfold, QuaPy, salmon	variability of the estimator	no coverage guarantee under prior shift or regularisation bias
Gibbs data augmentation	Leistedt et al., Starlight, rubias, SourceTracker	exact posterior	mixing slows as the fraction of missing information grows
NUTS on a transformed simplex	PixelPop, gppop, QuaPy’s Bayesian methods, PyFBU	exact posterior under any differentiable prior	geometry depends on the simplex transform; long tails for sparse α
Delta method; F-localisation	deconvolveR; Empirikos.jl	standard errors; frequentist intervals for functionals of g	the delta method ignores bias; F-localisation is conservative
Flows trained on a marginal likelihood	research code only	none (point estimates)	estimator bias of order 1/K; selection destabilises training

Approximate inference narrows intervals. Mean-field VB factorises the bin labels from ρ and so ignores the posterior coupling that broad likelihoods create. Rau et al. found that the Dirichlet scheme “yielded underestimated error bars, likely due to the restricted covariance structure of the Dirichlet that makes an independence assumption between neighboring bins” (Rau et al. 2022); HSC Y3 put the under-coverage at about 20% (Rau et al. 2023); and the general result for mixture weights is that VB intervals are “unrealistically narrow, especially if the components of the mixture model are not well separated” (Wang & Titterington 2005), which is the broad-likelihood regime. Point estimators fare worse. In ebnm’s evaluation the NPMLE “tended to have much poorer coverage” because its discrete prior “can greatly underestimate uncertainty” (Willwerscheid et al.); Efron’s log-spline intervals were the shortest in Ignatiadis & Wager’s simulations but “do not achieve nominal coverage, while all other methods do” (Ignatiadis & Wager); bootstrapped maximum-likelihood intervals for luminosity functions “can be too narrow” (Kelly, Fan & Vestergaard); and quantification’s bootstrap regions “cannot provide valid coverage guarantees under prior probability shift” (Moreo & Salvati, LQ 2025).

Regularisation bias breaks coverage. Kuusela & Panaretos showed that in unfolding “neither Bayesian credible intervals nor standard bootstrap confidence intervals succeed in achieving good frequentist coverage … due to the inherent bias”, and they recovered nearly nominal coverage only with an iteratively bias-corrected bootstrap (Kuusela & Panaretos). Coarse bins are a regulariser too: the response “depends on the unknown shape of the true spectrum”, wide bins introduce “a non-trivial bias”, and the recommended fix is to unfold on fine bins and then construct intervals for linear functionals of them (Stanley, Patil & Kuusela). Rau et al. found the same in photo-z work, where oversmoothing with coarse bins “can lead to posteriors that can be too narrow” (Rau et al. 2022). Pastry should therefore bin finely by default, report coarse summaries as functionals computed per posterior draw, and never let an EM iteration count stand in for a prior.

A single Dirichlet sees only shot noise. The hierarchical Dirichlet–multinomial treats objects as independent draws from one ρ, so its posterior width reflects counting statistics alone. In redshift calibration, plain Dirichlet resampling “is only propagating shot noise” while sample variance can be “an order of magnitude larger”, and both bootstrap and plain Dirichlet intervals fail to cover the truth where a variance-inflated three-step Dirichlet (3sDir) gives unit-normal pulls (Sánchez et al. 2020); DES Y3 notes that its Dirichlet step “only reproduces the variance expected from Poisson noise” (Myles et al. 2021). Other systematics dwarf counting noise: correcting Kepler candidates for reliability moves the Earth-analogue rate for GK stars from 0.034 to 0.015 (Bryson et al. 2020), and isochrone systematics in colour–magnitude fitting run two to four times the random errors (Zheng et al. 2025). At N ≥ 10⁷ the posterior width in well-identified directions shrinks as N^(−1/2) while these terms stay fixed, so overdispersion and likelihood-calibration hooks matter more than sampler precision.

Wrong inputs propagate silently. Feeding posteriors where likelihoods are expected bakes the interim prior into the answer. CHIPPR found a misspecified interim prior “disastrous” for every estimator it tested (Malz & Hogg); Malz concluded that “stacking must be abandoned” for future surveys (Malz 2021); Euclid’s plain stacked posteriors miss the accuracy requirement by five to ten times before PIT debiasing (Ilbert et al. 2021); and Rubin’s stacked machine-learning PDFs “show the imprint of the training sample” (Zhang et al. 2025). In GW, a parameter-estimation prior that vanishes at χ_p = 0 and 1 produced population features that “are not necessarily astrophysical in origin” (GWTC-4.0). Miscalibrated likelihoods bias results too. Leistedt & Hogg declined to inflate parallax errors because “overestimated errors typically cause excessive deconvolution” (Starlight); the EM quantifier assumes calibrated posteriors, “which is not true of modern neural networks” (Alexandari et al.); and QuaPy now ships likelihood tempering tuned on validation coverage, an implicit admission that classifier-derived likelihoods are miscalibrated (QuaPy _bayesian.py).

Sample-based likelihoods add Monte Carlo error. When per-object information arrives as posterior samples, each likelihood term is an importance-sampling estimate. The LVK now discards hyperparameter points where the total log-likelihood variance exceeds 1 (GWTC-4.0). Heinzel & Vitale show that the per-integral effective-sample-size cuts used earlier “are not necessarily safe”, that the plug-in likelihood is biased but correctable, and that information loss should stay below Ê ≲ 0.2 bits, a diagnostic their population-error package computes (Heinzel & Vitale; population-error). Talbot & Golomb call selection Monte Carlo “the limiting factor” for resolving narrow population features (Talbot & Golomb). The cuts are not innocent: in GWTC-4, changing them altered the inferred χ_p distribution, and for B-spline models the variance limit “was the main driver of smoothness” after the smoothing hyperparameter “consistently railed against prior boundaries” (GWTC-4.0). Likelihoods supplied on the same grid as the bins carry no such noise, an advantage pastry’s documentation should state plainly.

Mixing, identifiability and numerics bite at scale. Gibbs data augmentation and EM both slow down as the fraction of missing information grows (Tan 2023), which is exactly the broad-likelihood regime of faint galaxies and noisy parallaxes. In related categorical models data augmentation mixes dramatically worse than marginal Metropolis (Johndrow et al.), yet no published benchmark compares Gibbs, NUTS and VB for this model. Broad likelihoods also make adjacent bins strongly anti-correlated, and the naive histogram estimator “fails in the regime of strong inter-bin correlations or quasi-empty bins” (Leistedt et al. 2016); deconvolution converges at “very slow nonparametric rates” (Efron 2016); and “distinct source distributions may result in the same distribution over observations” (Vandegar et al.). Small numerical choices matter as well: a 2026 occurrence-rate script replaces log(0) by zero and silently drops objects without support (sblunt/eccentricities), and StarFormationHistories.jl had to remove an epsilon clamp that made “the gradient inconsistent with the objective” (CHANGELOG).

Together these failures define what “realistic uncertainty” must mean in pastry’s specification: exact posterior sampling as the reference, approximations labelled with their known bias, an overdispersion layer beyond the multinomial, input types that make interim priors and calibration explicit, Monte Carlo variance tracked whenever samples enter, posterior functionals reported per draw rather than as sums of marginals, and coverage demonstrated by simulation rather than asserted.

Every planned extension has precedent but no shared home

Each extension on pastry’s roadmap (generalised-Dirichlet and smoothness priors, selection effects, interim-prior correction, flow densities with marginal or joint training) has been implemented somewhere, usually once, in one field’s code. The subsections below collect what those implementations teach.

Conjugate priors add structure along an ordering, not smoothness

The Dirichlet’s covariances are all negative and fixed by the means plus a single precision (SciPy dirichlet). Its conjugate generalisations relax that rigidity while keeping Gibbs and EM updates closed-form given latent counts. The generalised Dirichlet has 2(K−1) parameters and “a more general covariance structure” (Generalized Dirichlet distribution), and a Dirichlet-tree “can give an independent variance to each p_k”, with positive correlation among leaves that share ancestors (Minka). Leistedt et al. noted in 2016 that a generalised-Dirichlet prior could encode inter-bin correlations but did “not consider this case” (Leistedt et al. 2016), and pastry will have to implement both families itself: the only Python generalised Dirichlet found is efax’s exponential-family form, which has no sampler (efax), while R’s MGLM supplies generalised-Dirichlet-multinomial fits for cross-checks (MGLM). Both families encode structure only through a bin ordering or a tree, so they suit hierarchical groupings and calibrated variances but cannot say “neighbouring bins should be similar”; smoothness needs the next family.

The concentration matters as much as the family. With three observations in 1,000 cells, the Jeffreys prior leaves the observed cells only 0.008 of the posterior mass, and Berger, Bernardo & Sun recommend α ≈ 1/K for sparse tables (Berger, Bernardo & Sun); broad likelihoods shrink the effective counts further. A symmetric Dirichlet over fine bins also biases coarse aggregates when groups contain unequal numbers of bins, the reporting-unit bias rubias corrects (rubias), so tying α to bin width as a base measure keeps re-binning prior-consistent. ashr’s device of encoding the prior as pseudo-observations cannot represent α_k < 1 (mix_opt.R), so sparse priors need samplers rather than penalised solvers.

Smoothness priors scale when their precision matrix is sparse

Smoothness lives in Gaussian priors on log-ratios (logistic-normal, random walks, ICAR, P-splines, GPs), which are not conjugate to the multinomial and therefore need gradient-based samplers on an unconstrained sum-to-zero vector, or Pólya–gamma stick-breaking augmentation where Gibbs is required (Linderman, Johnson & Adams). Every field has used them: random-walk priors on log heights for eccentricities (Hogg, Myers & Bovy 2010), a GP with marginalised hyperparameters for occurrence rates (exopop ms.tex), logistic GPs for n(z) (Rau, Wilson & Mandelbaum), P-spline difference penalties for GW masses (GWInferno), Tikhonov curvature with SURE-based strength selection in TUnfold (TUnfold.h), first- and second-difference penalties in ordinal quantification (qunfold likelihood.py) and Efron’s penalised log-spline g (deconv.R). Scaling depends on structure. A dense GP costs O(K³) per factorisation (Riihimäki & Vehtari), whereas PixelPop’s sparse ICAR runs NUTS over 10⁴ bins and gppop’s Kronecker GP over 8×10⁴ (PixelPop; gppop), and NumPyro already ships CAR, ZeroSumNormal and Hilbert-space GP components (NumPyro distributions).

The hard part is the smoothing strength. GWTC-4 fixed its B-spline penalty by hand after it railed against prior boundaries (GWTC-4.0); Hogg et al. called learning the smoothing strength and bin number jointly “beyond our scope” (Hogg, Myers & Bovy 2010); and PixelPop’s improper ICAR prior needs an explicit initialisation (PixelPop). Pastry should marginalise scales analytically where conjugate, as PixelPop’s Gamma reparametrisation does, report their posteriors, and offer empirical-Bayes fits of the scale, as Kuusela & Panaretos (Kuusela & Panaretos) and GP unfolding (Bozson, Cowan & Spanò) do, together with a prior-dominance diagnostic. Overdispersion belongs in the same module: a 3sDir-style concentration rescaling, per-field layers ρ_f ~ Dir(λρ̄), or a logistic-normal with a simulation-derived covariance each address the shot-noise-only failure described above (Sánchez et al. 2020).

Selection effects force a choice between two likelihoods

Mandel, Farr & Gair identify two classic mistakes: treating detection as independent of the data, which puts an extra detection probability in each object’s term and drops the normaliser, and including the normaliser while keeping the detection probability in the numerator (Mandel, Farr & Gair). Which form is right depends on the selection model. If selection is a cut on the same observed data that enter L_ik, it appears only in the normaliser Σ_k Q_k ρ_k; if detection is a separate random event whose efficiency depends on true properties, as in Kepler injection–recovery tests, Q_k also multiplies each object’s term (exopop ms.tex; sblunt/eccentricities). The Poisson-rate likelihood used for occurrence rates, Hess diagrams and GW populations reduces to the normalised-fraction form only under a scale-invariant prior on the total, since “This correspondence only holds with a 1/N scale-invariant prior” (Mandel, Farr & Gair). HEP models acceptance as Poisson thinning, f* = εf (Kuusela & Panaretos); deconvolveR reports both the selected and the thinning-corrected prior for zero-truncated counts (deconv.R); pygmmis needs completeness over all of feature space plus a callback predicting noise where nothing was observed (pygmmis); and for binned GW models the injection-estimated detection efficiency is linear in ρ, with closed-form Monte Carlo variance per bin (gppop). Selection functions should be expressed in observed attributes and defined “beyond the domain in q where we actually have cataloged objects” (Rix et al. 2021), and providers such as gaiaunlimited exist to supply them (gaiaunlimited). None of the quantification or stock-identification documentation reviewed mentions a selection term (QuaPy; rubias). Pastry should make users declare their selection model, offer both fraction and rate parameterisations, and return the intrinsic and selected populations side by side.

Interim-prior correction is a column rescaling with known failure modes

The mixture likelihood is invariant to rescaling each object’s row of L (mixsqp manual), so per-object evidences and normalisations can be dropped, but anything that varies across bins within a row cannot: the interim prior, bin widths and selection. Posteriors produced under an interim prior π⁰ therefore enter as L_ik ∝ p_ik/π⁰_k. GW codes apply this correction through a prior column (gwpopulation hyperpe.py), EMQ through training prevalences (QuaPy aggregative.py), and CHIPPR is the only photo-z code that applies it at all (Malz & Hogg). For posterior samples, each object’s term becomes an importance-sampling average with documented caveats. Foreman-Mackey et al. warn that it “could be a high variance estimator”, depending on the number of samples and the interim prior, whose support must be broader than the population’s (exopop ms.tex); Hogg et al. found that 50 samples per star were “only slightly worse” than 10⁵, for a 2,000-fold speed-up (Hogg, Myers & Bovy 2010); and GWTC-4’s χ_p episode shows what happens when the interim prior vanishes where the population has mass (GWTC-4.0). When the interim prior is implicit, as for machine-learning photo-z, quantification estimates it as the mean held-out posterior (Alexandari et al.; QuaPy aggregative.py). No photo-z n(z) code found takes per-galaxy posterior samples directly, and qp can only histogram them (qp basic usage), so sample inputs with interim-prior reweighting would be new in that domain. Pastry should distinguish likelihood inputs from posterior + interim_prior inputs, refuse posteriors that arrive without a declared prior, report per-object effective sample sizes and support mismatches, and include deliberately misspecified interim priors in its calibration tests.

Flows earn their cost only for samples in three or more dimensions

On gridded inputs, any smooth parameterisation ρ(θ) (splines, a logistic GP, a softmax-output network, or the bin masses of a bounded flow) plugs into the exact likelihood, which is precisely what Neural-g does with a softmax-output network (Wang et al. 2024); no IWAE or Monte Carlo is needed. The IWAE-style marginal objective is needed only when per-object information comes as samples or as continuous likelihoods in several dimensions, and there its biases are documented. The log-mean-exp estimator is biased low by O(1/K) in the number of samples K, and a Russian-roulette (SUMO-type) estimator removes the bias at small K; in Neural Empirical Bayes tests, direct Monte Carlo with K = 1024 beat both the ELBO and IWAE (Vandegar et al.). Dockhorn et al. found that “the variational bound is not tight, and biases us towards worse prior models”, with high run-to-run variance (Dockhorn et al.). The selection term −N ln P̂_det carries the opposite bias, so the full objective is no longer a bound, and Ruhe et al. dropped selection entirely because “Optimization was highly unstable when including the selection function” (Ruhe et al.). The most practical published recipe is Yan & Sanders’s: draw fresh samples from each star’s likelihood every epoch, pre-train on the noisy points, then deconvolve. The useful K plateaued between 4 and 128 depending on sample size, and training took 4 h at K = 32 and 16 h at K = 128 on one T4 GPU, against about 1 h for GPU extreme deconvolution (Yan & Sanders 2025).

Joint hierarchical sampling is cheap in pastry’s discrete case, because it is the conjugate Gibbs sampler, whereas GW analyses needed roughly 500-dimensional HMC to sample event parameters and hyperparameters jointly (Mancarella & Gerosa). For flows the joint route becomes Monte Carlo EM, alternating per-object draws under the current density with refits, which DiEM uses for diffusion priors while warning that “empirical Bayes is an ill-posed problem” and that a conservative initial prior is preferable (Rozet et al.). Uncertainty on learned densities remains unsolved: Ruhe et al. ran HMC over only 90 flow weights and judged it “not scalable” (Ruhe et al.), and amortised neural posteriors are often overconfident (Hermans et al.). Two traps deserve explicit documentation: pzflow’s convolve_errs option is noise augmentation rather than deconvolution (pzflow), and flows extrapolate beyond the support of samples and injections, so population claims should be restricted to supported regions (Iampieri & Mastrogiovanni). For backends, flowjax and zuko are the maintained MIT options (flowjax; zuko). Meta’s flow_matching carries a non-commercial CC BY-NC 4.0 licence (flow_matching LICENSE), nflows has not been released since 2020 (PyPI nflows), and conditional flow-matching losses can plateau while physical fidelity is still improving (ScatterPrism). Simulation-based approaches that need a forward simulator of each object’s data, such as neural posterior estimation or optimal-transport population fits, complement a likelihood-based core rather than replace it (Leyde et al.; Alsing et al. 2024).

Response matrices and point estimates widen the audience cheaply

Users without per-object PDFs can still be served. A response-matrix adapter turns HEP histograms into pastry inputs and restores the posterior sampling with nuisance parameters that PyFBU once provided (PyFBU), provided the response’s own Monte Carlo noise is handled honestly, since estimated response matrices “can inadvertently regularize the problem” (Zhu et al.). BayesianCC learns a confusion matrix jointly with the prevalences from labelled counts (Ziegler & Czyż), and StarFormationHistories.jl fits Hess-diagram counts as Poisson draws from a template matrix (SFH.jl docs). Epidemiology shows why the simplex matters here: plug-in prevalence corrections leave the simplex, and epiR recommends a Bayesian approach whenever apparent prevalence falls below one minus specificity (epiR manual), a case a simplex-constrained posterior handles automatically.

Specification lessons: layered API, exact defaults, built-in calibration

Copy emcee’s plain inputs and scikit-learn’s contracts, in four layers

emcee won astronomers because it works “with ‘black box’ models instead of structured graphical models” (Foreman-Mackey et al. 2019) and requires “hand-tuning of only 1 or 2 parameters” (Foreman-Mackey et al. 2013); its run_mcmc call, sample() generator and HDF5 backend with the same accessors as the in-memory one remain the model for run control (emcee backends). Its weaknesses are equally instructive: it seeds from NumPy’s global legacy random state (ensemble.py), and its most-discussed issues concern initialisation, autocorrelation estimates and parallelism (emcee issues). scikit-learn contributes the estimator contract, with hyperparameters stored unmodified in __init__ and “no logic, not even input validation” there, fit returning self, and learned attributes ending in underscores (scikit-learn developer guide), along with its principles of consistency, inspection, non-proliferation of classes and sensible defaults (Buitinck et al.). Full compliance is a poor fit, because pastry’s input is a likelihood matrix rather than a feature matrix and check_estimator would need many expected failures; a scikit-learn-style estimator is the right target. dynesty 3.0 replaced piles of tuning keywords with configured strategy objects (dynesty #495); nautilus declares priors as objects and resumes automatically from a checkpoint path (nautilus checkpointing); UltraNest reports which stopping criterion was met and why (UltraNest tutorial); and tinygp shows the value of exposing a model as plain functions that other inference libraries can call (Why tinygp?).

These lessons compose into four layers: a one-call pastry.fit(...) with documented defaults; estimator classes that accept configured prior, selection and sampler objects; a generator and checkpointing layer for run control, with callbacks in the style of scikit-learn’s accepted callback proposal (SLEP023); and array-level functions (log-density, gradient, EM step) plus base classes for other libraries. A sketch consistent with these exemplars, not a settled design:

import pastry

res = pastry.fit(logL, bins=edges, rng=0)                 # layer 0: one call, documented defaults
res.summary(); res.plot(); dt = res.to_datatree()

model = pastry.HistogramModel(                             # layer 1: estimator + strategy objects
    prior=pastry.priors.ICAR(scale="marginalise"),
    selection=pastry.Selection(Q, mode="true-property"),
    sampler=pastry.samplers.NUTS(), rng=0)
model.fit(pastry.inputs.from_samples(samples, bins=edges, interim_prior=pi0))
model.fractions_, model.posterior_

for state in model.sample(checkpoint="run.h5"): ...       # layer 2: run control
lp = pastry.logdensity(rho, logL)                          # layer 3: plain functions for NumPyro, BlackJAX, emcee

Inputs should centre on a canonical log-likelihood matrix, dense, sparse or chunked, with optional per-object weights, as mixsqp accepts log input (mixsqp manual). Adapters should cover gridded PDFs integrated over bins (RAIL evaluates them at grid points instead; var_inf.py), ragged posterior samples with an interim prior, padded without leaking into the counts as PixelPop learned after downsampling “discards samples from well-measured events” (PixelPop ragged-samples migration), Gaussian (x, σ) pairs, response matrices with counts, and qp ensembles (qp). Smooth densities plug in through a small protocol, bin masses for gridded inputs and a log-density for sample inputs, following the distribution objects of flowjax and zuko (flowjax; zuko).

Outputs should be a plain draws × K array plus an xarray.DataTree that ArviZ reads natively (Martin et al. 2026), per-object membership posteriors, posterior-predictive counts, derived quantities computed per draw, and region objects with coverage helpers in QuaPy’s style (QuaPy confidence.py). Diagnostics belong on every result: R̂ and ESS, the KKT gap of any MAP solution, the fraction of missing information, an Efron-style ratio of prior to data information (deconv.R), and flags for bins with negligible data support, which GWTC-5 hatches when they hold under 1% of samples (GWTC-5.0). Covariances on the simplex are rank-deficient, so comparison utilities should work in K − 1 log-ratio coordinates (Martinez Lopez et al.), and draws of several summaries should be exported jointly so correlations survive, as the PDG review urges; following its advice to publish what makes forward folding possible, pastry should also save the likelihood matrix or its generator next to the draws (PDG statistics review). Exports should target qp, CosmoSIS realisation files and nz_prior for cosmology (nz_prior; nz_multirank) and a popsummary-like HDF5 layout (popsummary). Randomness should flow through one rng= keyword per SPEC 7, rejecting legacy RandomState objects as PyMC does (SPEC 7; pymc.sample). The surveyed tools also show what to avoid: global backend switches (gwpopulation backend.py), mutable attribute-style configuration (PyFBU), commercial solver requirements (Willwerscheid et al.), and print-based progress with import-time warnings (GLMixture.py).

Default to exact sampling, label approximations, compress before scaling

The inference stack follows from the comparisons above. A MAP solution by SQUAREM-accelerated EM or an active-set Newton method gives a fast first look; first-order methods remained at least 0.35 log-likelihood units short after 1,000 iterations on a problem mix-SQP solved (Kim et al. 2020), and every MAP result should carry a KKT certificate, because EM can never revive a bin it has zeroed. The default posterior engine should be NUTS on a pastry-owned isometric-log-ratio simplex, the softmax of an orthonormal sum-to-zero vector that Stan adopted in version 2.37 (Stan transforms), rather than whatever a backend defaults to: NumPyro uses stick-breaking, PyMC a centred log-ratio and TFP an additive log-ratio (NumPyro distributions; PyMC Dirichlet; TFP SoftmaxCentered), and Stan switched without publishing a benchmark (Stan 2.37 release). Conjugate Gibbs serves as the exact cross-check and as the engine for sparse priors with α < 1, where unconstrained parameterisations develop long tails and Langevin schemes are biased near the boundary (Baker et al.). Mean-field VB stays available but is labelled overconfident, and evidences come from nested sampling at small K or from Chib’s estimator applied to conjugate Gibbs output (Chib 1995). Because the fraction of missing information governs both EM’s speed and Gibbs mixing (Tan 2023), estimating it at the MAP lets pastry recommend Gibbs or NUTS automatically.

At pastry’s target scales memory, not arithmetic, is the binding constraint. The Gaia TGAS likelihood matrix for 1.36 million stars on a 100×100 grid alone occupied about 109 GB, and the dual augmented-Lagrangian solver needed 472 minutes on a 32-core, 256 GB workstation (Zhang et al. 2024); by simple arithmetic, a dense float32 matrix at N = 10⁸ and K = 10³ needs 400 GB. The remedies already exist elsewhere. Objects with identical rows can be collapsed into weighted classes, as salmon does with equivalence classes (salmon) and Empirikos does with multiplicities (Empirikos NPMLE.jl); likelihoods can be stored sparsely above a log threshold, as npeb does (GLMixture.py); low rank can be exploited in one dimension, which mix-SQP relies on and which fails in two (Zhang et al. 2024); and per-object likelihoods should be computed once over fixed bins, Leistedt & Hogg’s route to 1.4 million stars (Starlight). Minibatch stochastic-gradient MCMC is a last resort. The core should depend only on NumPy and SciPy, with JAX samplers and flows as optional extras; TensorFlow Probability is a poor dependency, with no stable release since November 2024 (PyPI tensorflow-probability).

Ship the validation suite that no current tool has

Most tools surveyed test that code runs, not that intervals cover. RAIL’s summarizer tests are smoke tests (rail_base tests), gppop’s test directory is empty and FIGARO’s suite is 54 lines (gppop; FIGARO), and deconvolveR and npeb have no tests (deconvolveR; NPEB). A credible suite has layers. Exact cases run on every commit: one-hot likelihoods must return exactly Dir(α + counts), two-bin posteriors can be checked against quadrature, the EM fixed point against brute-force optimisation, gradients against finite differences, and rescaling any object’s row must leave the posterior unchanged. Gibbs code gets Geweke’s “getting it right” test, which compares forward simulation with a chain that alternates MCMC updates and fresh data (LIPS blog). Simulation-based calibration should rank not only each ρ_k but also the joint log-likelihood, because parameter ranks alone “could never detect large classes of problems including when the posterior is equal to the prior” (Modrák et al.). It should ship as a user-facing validation module driven by the user’s own per-object noise model, since the ArviZ team’s simuk works only with PyMC, Bambi or NumPyro models (simuk).

Coverage suites can borrow designs wholesale: Ignatiadis & Wager’s spiky and smooth priors with pointwise and simultaneous coverage (Ignatiadis & Wager), quantification’s artificial-prevalence protocols (QuaPy), HEP closure and stress tests (Kuusela & Panaretos), posterior-predictive checks of the observed data in the spirit of HEP’s preference for folding over unfolding (Cousins, May & Sun), and golden-value comparisons against mixsqp, REBayes and Empirikos on their bundled datasets (Empirikos datasets). Statistical tests should use dynesty’s tolerance of four times the run’s own Monte Carlo error, with an environment variable for looping over seeds (dynesty test_gau.py; dynesty utils.py), slow suites behind markers, and sequential tests rather than blanket reruns (Gandy & Scott). Every run should report rank-normalised R̂ and bulk and tail ESS (Vehtari et al.), which ArviZ’s NumPy-only array layer computes (arviz-stats). The centrepiece should be the benchmark nobody has published: Gibbs against NUTS, VB, EM-plus-Laplace, bootstrap and flows, with coverage tracked as likelihood breadth and missing information grow, under correct and deliberately misspecified interim priors and selection.

Package for longevity: SPEC 0, open history, optional extras

The PyPI name pastry belongs to an unrelated 2018 API wrapper for Chef, while pastry-inference was free on 22 October 2026 (PyPI pastry). Under SPEC 0, support for Python 3.12 and NumPy 2.2 ends in the fourth quarter of 2026 (SPEC 0), so a first release should require Python ≥ 3.12, or ≥ 3.13 if it ships after that quarter, together with current NumPy and SciPy. The Scientific Python Development Guide recommends hatchling with an SPDX licence identifier, ruff, strict pytest settings that turn warnings into errors, nox, uv-based CI and Trusted Publishing (SPDG packaging; SPDG pytest), and Sphinx with MyST-NB can execute every tutorial on each documentation build, so tutorials double as integration tests (MyST-NB). MkDocs has had no release since August 2024 and Material for MkDocs is in maintenance mode (PyPI mkdocs; pydevtools handbook). JOSS now rejects without review any submission lacking six months of public, incremental development history and evidence of research use, and it requires an AI-usage disclosure (JOSS submission rules), so the repository should be public early and paired with a science application.

Durability is the recurring failure across fields. CRAN archived rubias because the maintainer’s email bounced (CRAN rubias), PyFBU died with Theano (pyFBU/fbu), several of the neural repositories carry no licence (density-deconvolution), and emcee kept its core maintainable by moving MPI pools and parallel tempering into other packages because “there aren’t enough resources to maintain” them (emcee upgrade notes). Pastry should declare an OSI licence, keep releases multi-maintainer, follow NumPy’s policy of at least one year between deprecation and removal (NEP 23), require ArviZ 1.x, since version 1.0 removed InferenceData and changed the default interval probability to 0.89 (ArviZ migration guide), and cut scope as deliberately as emcee did.

The table condenses the specification decisions that the evidence supports.

Decision	Recommendation	Evidence
Distribution name	publish as pastry-inference, import as pastry	PyPI pastry
Canonical input	log-likelihood matrix (dense, sparse or chunked) with weights; adapters for bin-integrated PDFs, ragged samples + interim prior, Gaussian (x, σ), response matrix + counts, qp ensembles	mixsqp; PixelPop migration; RooUnfold
Interim prior	mandatory for posterior inputs; per-object ESS and support checks	Malz & Hogg; exopop ms.tex; GWTC-4.0
Selection	declared model (cut on observed data, or detection on true properties); fraction and rate parameterisations; intrinsic and selected outputs	Mandel, Farr & Gair; Kuusela & Panaretos
Priors	Dirichlet with explicit α presets (1, ½, 1/K, proportional to bin width); generalised Dirichlet and Dirichlet-tree; ICAR, random-walk, HSGP and Kronecker GP on sum-to-zero log-ratios; an overdispersion layer	Berger, Bernardo & Sun; rubias; Sánchez et al.; PixelPop
Inference	MAP with KKT certificate; NUTS on an isometric-log-ratio simplex by default; conjugate Gibbs as cross-check; VB labelled approximate; sampler recommendation from missing information	Kim et al.; Stan transforms; Rau et al. 2023
Scale	deduplicate rows, sparse storage, chunked streaming; never materialise dense L by default	Zhang et al.; salmon
Outputs	draws × K array and ArviZ DataTree; memberships; predictive counts; diagnostics; qp, CosmoSIS and popsummary-style exports	ArviZ 1.x; popsummary; nz_prior
Randomness	single rng= keyword, spawned per chain	SPEC 7
Validation	exact cases, Geweke test, SBC with joint log-likelihood, coverage suite, golden values against R and Julia solvers, published engine benchmark	Modrák et al.; Ignatiadis & Wager
Packaging	Python ≥ 3.12 (≥ 3.13 after 2026-Q4); NumPy/SciPy core; [jax], [flows], [plot], [io] extras; hatchling, ruff, nox, uv, Trusted Publishing; Sphinx + MyST-NB	SPEC 0; SPDG

Conclusion

The usual story of scientific software is inverted here: the statistics converged long ago across fields, while the engineering fragmented and decayed. That changes what pastry has to prove. Its model needs no defence, since it is the conjugate core of photo-z, colour–magnitude, occurrence-rate, unfolding, empirical-Bayes, quantification and stock-identification analyses alike; its uncertainties do, because every maintained alternative either reports none or reports intervals its own authors call too narrow. The highest-value deliverable after the core is therefore the cross-engine coverage benchmark that no field has published. That one study would justify pastry’s default engine, supply the research-use evidence JOSS now demands, and give RAIL, QuaPy and HEP users a concrete reason to switch, while the stacking-equals-one-EM-step identity gives them a familiar way in.

The second implication concerns where the error budget will sit. At the sample sizes pastry targets, counting noise becomes negligible next to sample variance, miscalibrated per-object PDFs, implicit interim priors and selection, so pastry’s competitive edge is the ability to express and test those systematics, through overdispersion layers, calibration and tempering hooks, declared selection models and calibration runs fed deliberately wrong inputs, rather than raw sampler speed. Neural extensions follow the same logic. On gridded inputs a flexible density needs no IWAE machinery and inherits exact posterior uncertainty, so marginal-likelihood flow training belongs in a later release, reserved for sample-only, multi-dimensional problems, where its documented biases and missing uncertainty quantification make it a research feature rather than a default.