Posterior Inference

calibrate() returns a point estimate; Uncertainty & Error puts a Gaussian around it. prophys.inference gives the posterior itself, so a skewed, bounded or multi-modal posterior is reported as it is rather than as the best-fitting ellipse around its mode.

Requires the inference extra:

pip install "prophys[inference]"

The bridge is deliberately thin. A prophys model already supplies everything a gradient-based sampler needs — a differentiable log-likelihood over an unconstrained parameter vector — so the work is exposing that as a NumPyro potential function, declaring priors, and packaging the draws in the shape the rest of the engine consumes.

Sampling

import prophys as prp
from prophys.inference import sample_posterior

fit = prp.calibrate(compiled, "y", data, n_steps=2000)
posterior = sample_posterior(
    compiled,
    {"y": data},
    init_params=fit.raw_params,   # start at the mode: shorter warmup
    n_warmup=500,
    n_samples=1000,
    n_chains=2,
)
print(posterior.summary())

PosteriorSummary:
  2000 draws over 2 chain(s), 0 divergence(s)
  mu = 2.02293 +/- 0.1022  [1.83093, 2.21944]  r_hat 1.001, ess 1000
  sigma = 1.45369 +/- 0.07293  [1.31567, 1.60417]  r_hat 0.999, ess 688

Sampling happens in unconstrained space, where the geometry is smooth and unbounded and NUTS performs well; the draws are mapped back through Param.forward before being reported, so a positive scale parameter’s posterior never contains a negative value. Both forms are available: draws (constrained) and raw_draws.

Convergence

summary() reports the three things that decide whether the draws describe the posterior at all:

• r_hat — split Gelman-Rubin, comparing chains and each chain’s two halves, so a single chain that drifted over its run is detected. Above about 1.01 the chains have not mixed.

• ess — effective sample size per parameter. A few hundred suffices for a posterior mean; a tail quantile needs more.

• divergences — NUTS transitions that could not follow the posterior’s geometry. Any at all mean the draws may be biased in a way more sampling does not fix.

converged requires all three to be healthy. A single-chain run reports False whatever the draws look like: with one chain there is no between-chain comparison to make.

Chains started from init_params are jittered apart from each other. Identical starting points can agree for many draws without having explored anything, which is precisely the failure r_hat exists to detect.

Priors

Priors are declared per parameter, on the unconstrained value — the scale the sampler works on. For a Param with bounds or a softplus transform, that scale is the pre-image of its constrained range.

posterior = sample_posterior(
    compiled,
    {"y": data},
    priors={
        "mu": ("normal", {"loc": 0.0, "scale": 5.0}),
        "sigma": ("cauchy", {"loc": 0.0, "scale": 2.5}),
    },
)

Named priors are normal, cauchy, student_t, uniform and laplace; anything else is passed as a NumPyro distribution instance. Parameters without an entry get a wide Normal(0, 10): sampling needs a proper prior to have a normalizable target, and a flat improper one leaves an unbounded parameter free to wander where the likelihood is flat.

Feeding the rest of the engine

draws is the same mapping sample_parameters() produces, so every predictive score and posterior-predictive function accepts it directly:

prp.waic(compiled, "y", data, posterior.thin(400))
prp.loo(compiled, "y", data, posterior.thin(400))
prp.posterior_predictive_check(
    compiled, "y", data, statistic="max", param_draws=posterior.thin(400)
)

thin() takes a random subset rather than every k-th draw: with a residual periodicity in the chain, systematic thinning can land repeatedly on the same phase of it.

Embedding in a larger NumPyro program

potential_function() exposes the model’s negative log-likelihood on its own, for use as one factor among several — inside a hierarchical model whose upper levels are declared in NumPyro directly, or alongside another likelihood term:

from prophys.inference import potential_function
import numpyro

potential = potential_function(compiled, {"y": data})

def model():
    params = {"mu": numpyro.sample("mu", dist.Normal(0.0, 5.0)), ...}
    numpyro.factor("prophys", -potential(params))

It is the same objective finetune() minimizes, so a posterior sampled from it and a point fit describe one model rather than two.

prophys.inference.sample_posterior(compiled, observations, *, priors=None, init_params=None, weights=None, trainable=None, frozen=None, n_warmup=500, n_samples=1000, n_chains=2, seed=0, target_accept=0.8, progress_bar=False)

Sample the posterior over compiled’s parameters with NUTS.

Parameters

priors:

Per-parameter priors over the unconstrained value, as NumPyro distributions or (name, kwargs) pairs (see _PRIOR_FACTORIES). Parameters without an entry get a wide Normal(0, 10). Priors are declared on the unconstrained scale because that is where the sampler works; for a Param with bounds or a softplus transform, that scale is the pre-image of its constrained range, so a prior there is a prior on the transformed parameter.

init_params:

Unconstrained starting values, typically a point fit’s raw_params. Starting from the mode shortens warmup considerably and makes an unconverged run easier to recognize, since the chains begin somewhere the posterior actually has mass.

n_chains:

Two or more, so r_hat can detect chains that disagree. A single chain gives draws but no way to tell whether they describe the posterior.

Returns a Posterior whose draws plug directly into prophys.engine.scoring.waic(), prophys.engine.scoring.loo() and prophys.engine.scoring.posterior_predictive_check().

Parameters:

• observations (Mapping[str, Any])

• priors (Mapping[str, Any] | None)

• init_params (Mapping[str, Any] | None)

• weights (Mapping[str, float] | None)

• trainable (Iterable[str] | None)

• frozen (Iterable[str] | None)

• n_warmup (int)

• n_samples (int)

• n_chains (int)

• seed (int)

• target_accept (float)

• progress_bar (bool)

Return type:

Posterior

class prophys.inference.Posterior(draws, raw_draws, extra=<factory>, n_chains=1)

A posterior sample over a compiled model’s parameters.

Not a Report: this carries the full draws, which are data rather than a summary. Call summary() for the report, or pass draws straight to prophys.engine.scoring.waic() and friends.

Parameters:

• draws (dict[str, ndarray])

• raw_draws (dict[str, ndarray])

• extra (dict[str, Any])

• n_chains (int)

draws: dict[str, ndarray]

Constrained-space posterior draws, one array of shape (n_draws,) per parameter. This is the mapping every predictive score and posterior-predictive function in the engine accepts.

raw_draws: dict[str, ndarray]

The same draws in unconstrained (sampler) space.

extra: dict[str, Any]

divergences, tree depth, step size, and so on.

Type:

Sampler diagnostics as returned by NumPyro

summary()

Per-parameter means, credible intervals, r_hat and effective sample sizes.

Return type:

PosteriorSummary

thin(n, seed=0)

n draws sampled without replacement, for a scoring pass that does not need the whole chain.

Drawn at random rather than by taking every k-th draw: with a residual periodicity in the chain, systematic thinning can land repeatedly on the same phase of it, while a random subset cannot.

Parameters:

• n (int)

• seed (int)

Return type:

dict[str, ndarray]

class prophys.inference.PosteriorSummary(names, mean, std, quantiles, r_hat, ess, n_draws, n_chains, divergences)

Per-parameter summary of a posterior sample.

Parameters:

• names (list[str])

• mean (dict[str, float])

• std (dict[str, float])

• quantiles (dict[str, list[float]])

• r_hat (dict[str, float])

• ess (dict[str, float])

• n_draws (int)

• n_chains (int)

• divergences (int)

quantiles: dict[str, list[float]]

The 2.5%, 25%, 50%, 75% and 97.5% posterior quantiles per parameter.

r_hat: dict[str, float]

Gelman-Rubin potential scale reduction across chains. Values above about 1.01 mean the chains have not mixed and the draws do not yet describe one distribution. nan with a single chain, which cannot detect this at all.

ess: dict[str, float]

Effective sample size per parameter, accounting for autocorrelation within the chains. A few hundred is enough for a posterior mean; a tail quantile needs more.

n_draws: int

Post-warmup draws, summed across chains.

divergences: int

NUTS transitions that diverged. Any at all mean the sampler could not follow the posterior’s geometry somewhere, and the draws may be biased in a way no amount of further sampling fixes.

credible_interval(name)

The 95% central credible interval for one parameter.

Parameters:

name (str)

Return type:

tuple[float, float]

property converged: bool

Whether every parameter’s r_hat is below 1.01 and no transition diverged. A single-chain run cannot establish this and reports False.

prophys.inference.potential_function(compiled, observations, *, weights=None, trainable=None, frozen=None)

The model’s negative log-likelihood as a NumPyro potential function.

Takes a dict of unconstrained parameter values and returns a scalar. This is the same objective finetune minimizes, so a posterior sampled from it and a point fit are describing one model, not two.

Exposed separately from sample_posterior() so the model can be embedded in a larger NumPyro program — as one factor among several, or inside a hierarchical model whose upper levels are declared in NumPyro directly.

Parameters:

• observations (Mapping[str, Any])

• weights (Mapping[str, float] | None)

• trainable (Iterable[str] | None)

• frozen (Iterable[str] | None)

Return type:

Callable[[Mapping[str, Any]], Any]