Posterior Inference

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)[source]#

Sample the posterior over a compiled model’s parameters with NUTS.

NUTS explores the exact log density of the model in unconstrained space, and the draws are mapped back through Param.forward, so a skewed, bounded or multimodal posterior is reported as it is. Chains run sequentially. Requires the inference extra (pip install "prophys[inference]").

Parameters:

compiled – The model.
observations (Mapping[str, Any]) – Attribute name to observations.
priors (Mapping[str, Any] | None) – Priors per parameter over the unconstrained value, as NumPyro distributions or (name, kwargs) pairs with name one of "normal", "cauchy", "student_t", "uniform", "laplace" (for example ("normal", {"loc": 0.0, "scale": 2.0})). Parameters without an entry get a wide Normal(0, 10). For a Param with bounds or a softplus transform, the unconstrained scale is the preimage of its constrained range, so a prior there is a prior on the transformed parameter.
init_params (Mapping[str, Any] | None) – Unconstrained starting values, typically a point fit’s raw_params. Starting from the mode shortens warmup and makes an unconverged run easier to recognize. With several chains the starts are jittered.
weights (Mapping[str, float] | None) – Optional likelihood weights per attribute.
trainable (Iterable[str] | None) – Names or glob patterns of the parameters to sample.
frozen (Iterable[str] | None) – Names or glob patterns of parameters not sampled. Parameters not sampled are held at their declared initial values (not at init_params).
n_warmup (int) – Warmup (adaptation) iterations per chain.
n_samples (int) – Retained draws per chain.
n_chains (int) – 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.
seed (int) – PRNG seed.
target_accept (float) – NUTS target acceptance probability. Raise it (for example to 0.95) when the summary reports divergences.
progress_bar (bool) – Show NumPyro’s progress bar.

Returns:

A Posterior whose draws plug directly into waic(), loo() and posterior_predictive_check().

Raises:

ImportError – If NumPyro is not installed.
ValueError – If no parameter is selected or n_chains < 1.
KeyError – If a prior names a parameter that is not sampled.

Return type:

Posterior

Example

>>> fit = finetune(compiled, {"y": y_obs})
>>> post = sample_posterior(compiled, {"y": y_obs}, init_params=fit.raw_params)
>>> print(post.summary())

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

Posterior draws over a compiled model’s parameters, from sample_posterior().

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]#: Posterior draws in constrained space, one array of shape (n_draws,) per scalar parameter ((n_draws, *shape) for an array parameter), chains concatenated. 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]#: Extra sampler fields returned by NumPyro, keyed by name. Only "diverging" (one flag per draw) is requested.

n_chains: int = 1#: Number of chains the draws came from.

summary()[source]#

Summarize the draws: means, credible intervals, r_hat and effective sample sizes per parameter.

Returns:: A PosteriorSummary report.
Return type:: PosteriorSummary

thin(n, seed=0)[source]#

Return a random subset of n draws, for a scoring pass that does not need the whole chain.

Draws are sampled without replacement, and all draws are returned when n is at least the number available.

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)[source]#

Posterior means, credible intervals and convergence diagnostics per parameter.

Returned by Posterior.summary(). Check converged before reading the numbers.

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)

names: list[str]#: Parameter names, sorted.

mean: dict[str, float]#: Posterior mean per parameter (constrained space).

std: dict[str, float]#: Posterior standard deviation per parameter.

quantiles: dict[str, list[float]]#: The 2.5%, 25%, 50%, 75% and 97.5% posterior quantiles per parameter.

r_hat: dict[str, float]#: Split potential scale reduction of Gelman and Rubin, comparing the halves of every chain. Values above about 1.01 mean the chains have not mixed and the draws do not yet describe one distribution. With a single chain it compares only the two halves of that chain, and it is nan with fewer than four draws per chain.

ess: dict[str, float]#: Effective sample size per parameter, accounting for autocorrelation within the chains. A few hundred is enough for a posterior mean, and a tail quantile needs more.

n_draws: int#: Draws after warmup, summed across chains.

n_chains: int#: Number of 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)[source]#

Return the 95% central credible interval (2.5%, 97.5%) of 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 run with a single chain cannot establish this and reports False.

prophys.inference.potential_function(compiled, observations, *, weights=None, trainable=None, frozen=None)[source]#

Return the model’s negative log likelihood as a NumPyro potential function.

This is the same objective finetune() minimizes (without a prior term), so a posterior sampled from it and a point fit describe 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:

compiled – The model.
observations (Mapping[str, Any]) – Attribute name to observations.
weights (Mapping[str, float] | None) – Optional likelihood weights per attribute.
trainable (Iterable[str] | None) – Names or glob patterns of the parameters the function takes. The others are held at their declared initial values.
frozen (Iterable[str] | None) – Names or glob patterns of parameters held at their declared initial values.

Returns:

A function of a dict of unconstrained values of the selected parameters, returning the scalar negative log likelihood.

Return type:

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

Example

>>> potential = potential_function(compiled, {"y": y_obs})
>>> numpyro.factor("prophys_likelihood", -potential({"mu": mu, "log_sigma": s}))

Posterior Inference

Contents

Posterior Inference#

Sampling#

Convergence#

Priors#

Feeding the rest of the engine#

Embedding in a larger NumPyro program#