Uncertainty & Error
Every number a model reports carries uncertainty from two distinct sources, and prophys reports them separately because they are reduced by different things:
Source |
Reduced by |
Reported by |
|---|---|---|
Monte-Carlo error — the statistic was estimated from a finite draw |
More samples |
|
Parameter uncertainty — the parameters were estimated from finite data |
More data |
A model can have negligible Monte-Carlo error and large parameter uncertainty, or the reverse. Adding samples never reduces the second; adding data never reduces the first.
Monte-Carlo error
CompiledModel.assess returns
the estimate together with the sampling error of the same draw:
import prophys as prp
compiled = model.compile(n_samples=20_000)
compiled.assess("loss") # the mean
compiled.assess("loss", "prob", level=2.0) # P(loss <= 2)
compiled.assess("loss", "quantile", level=0.99) # the 99th percentile
compiled.assess("loss", "cvar", level=0.95) # tail risk
Each returns an EstimateWithError (or a
TailRiskEstimate for CVaR):
EstimateWithError:
quantile[loss] at 0.99 = 10.3263 +/- 0.7063
95% interval = [8.94206, 11.7106]
relative error = 6.840%
n_samples = 20000, ess = 20000.0 (iact 1.00)
Three estimators sit behind these, one per shape of statistic:
A mean (and an exceedance probability, which is the mean of an indicator) uses non-overlapping batch means plus an initial-positive-sequence autocorrelation time, so correlated draws — a
Processpath, an MCMC chain — report a larger error and a lower effective sample size than independent ones.A quantile uses the Maritz-Jarrett order-statistic estimator, which needs no density estimate at the quantile. The classical \(\sqrt{q(1-q)/n}/f(x_q)\) form requires \(f\) precisely where the sample is thinnest.
A CVaR uses the influence function of the Rockafellar-Uryasev functional, which includes the uncertainty of the VaR threshold. The standard error of the mean of the tail draws omits that term and is correspondingly optimistic.
Sample-size planning
required_samples() converts an achieved error into the
draws a target precision would take, on the \(1/\sqrt{n}\) scaling:
estimate = compiled.assess("loss", "cvar", level=0.99)
prp.required_samples(estimate, target_relative_error=0.01)
Parameter uncertainty
After a fit, parameter_uncertainty() inverts the
observed information — the Hessian of the same negative log-likelihood the
fit minimized — at the optimum:
result = prp.calibrate(compiled, "y", data, n_steps=3000)
u = prp.parameter_uncertainty(compiled, {"y": data}, result=result)
print(u)
ParameterUncertainty:
fitted on 500 observations, log-likelihood -1054.32
mu = 2.92676 +/- 0.09951 [2.73172, 3.1218]
sigma = 1.99021 +/- 0.07036 [1.8523, 2.12811]
Standard errors are reported in constrained space — the values the model graph uses and the reader interprets — obtained from the unconstrained-space Hessian by the delta method, so a positive scale parameter’s interval never straddles zero.
The Hessian is central-differenced from the exact analytic gradient. A
native kernel is registered as a jax.custom_vjp whose backward pass is
itself a plain FFI call with no rule of its own, so a model graph is
differentiable exactly once; differencing an exact gradient leaves a
relative accuracy of about \(\varepsilon^{2/3}\), far below the sampling
uncertainty being reported.
Unidentified parameters
At a genuine minimum the Hessian is positive definite. When it is not,
parameter_uncertainty refuses rather than returning numbers:
ValueError: The Hessian at the reported optimum is not positive definite,
so no covariance exists there. The optimizer has not reached a minimum, or
some parameter (or combination) is not identified by the data.
For a model that is identified but only weakly,
collinear_pairs() lists the
parameter pairs the data constrain only in combination:
for a, b, rho in u.collinear_pairs():
print(f"{a} and {b} are correlated at {rho:+.3f}")
Propagating into a derived quantity
Parameter uncertainty rarely matters on its own; what matters is its effect
on the quantity being reported. propagate() carries the
covariance through the Jacobian of any function of the parameters:
design_value = prp.propagate(u, lambda p: p["mu"] + 1.645 * p["sigma"])
print(design_value.interval(0.95))
This is a first-order propagation: exact for a linear function, and an
approximation whose quality falls with the curvature of the function over
the parameter uncertainty. For a strongly non-linear derived quantity,
sample the parameters with
sample_parameters() and push the draws through
instead.
Profile likelihood
The Laplace interval assumes the log-likelihood is quadratic near its
maximum. profile_likelihood() makes no such assumption:
it fixes one parameter across a grid, re-optimizes every other parameter at
each point, and reads the interval off the likelihood ratio.
import numpy as np
profile = prp.profile_likelihood(
compiled, {"y": data}, "sigma", np.linspace(1.5, 2.5, 21), result=result
)
print(profile.interval)
For a location parameter, where the likelihood is exactly quadratic, the two intervals coincide; for a bounded or skewed parameter the profile interval is asymmetric and the Laplace one is not. It costs one optimization per grid point.
Beyond the Gaussian approximation
Both routes above summarize the posterior by a Gaussian around its mode. For the posterior itself — skewed, bounded, or multi-modal as it may be — see Posterior Inference.
- prophys.engine.assess(compiled, attribute_name, statistic='mean', *, level=None, n_samples=None, params=None, key=None, n_batches=0)[source]
Estimate a statistic of a model attribute together with its Monte-Carlo error.
Parameters
- compiled:
A CompiledModel.
- attribute_name:
Which attribute to sample.
- statistic:
One of
STATISTICS."prob"needs level as the threshold (it estimatesP(attribute <= level));"quantile"and"cvar"need level as the probability level.- n_samples:
Draws to use; defaults to the model’s compiled n_samples.
- key:
PRNG key; defaults to the model’s compiled seed, which makes the assessment reproducible but ties repeated calls to one sample. Pass distinct keys to see the between-run variation directly.
All four statistics are estimated from a single sample of the attribute, so the Monte-Carlo error reported is that of the same draw the estimate came from.
- class prophys.engine.EstimateWithError(statistic, value, std_error, n_samples, effective_sample_size, autocorrelation_time, attribute=None, level=None)[source]
A Monte-Carlo estimate with its sampling error.
- Parameters:
- effective_sample_size: float
Draws’ worth of independent information, n / iact. Equal to n_samples for independent draws, lower for a correlated series.
- level: float | None = None
The q or alpha of a quantile/CVaR estimate, or the threshold of a prob estimate.
- property relative_error: float
Standard error as a fraction of the estimate. nan for an estimate of zero, where no relative statement is meaningful.
- class prophys.engine.TailRiskEstimate(alpha, var, cvar, std_error, n_tail, n_samples, attribute=None)[source]
Value-at-Risk and Conditional Value-at-Risk with the CVaR’s error.
- Parameters:
- prophys.engine.mean_error(draws, *, n_batches=0, name='mean')[source]
Sample mean of draws with its batch-means standard error.
n_batches selects the number of blocks; the default (0) uses about sqrt(n) blocks of about sqrt(n) draws. Correlation up to a block length is reflected in the reported error.
- Parameters:
- Return type:
- prophys.engine.quantile_error(draws, q)[source]
Sample q-quantile of draws with its Maritz-Jarrett standard error.
The error is an order-statistic estimate, so it stays defined in the far tail where a density-plug-in estimate would need the density precisely where the sample is thinnest.
- Parameters:
- Return type:
- prophys.engine.cvar_error(draws, alpha)[source]
CVaR of draws at alpha, with its influence-function error.
- Parameters:
- Return type:
- prophys.engine.required_samples(estimate, target_relative_error)[source]
Draws needed to reach target_relative_error, extrapolated from an achieved estimate.
Monte-Carlo error falls as 1/sqrt(n), so the requirement scales as the square of the ratio of achieved to target error. The result assumes the estimator’s variance does not itself change with n, which holds for a mean and is approximate for a tail quantile, where the far tail is better resolved at larger n than the extrapolation credits.
- Parameters:
estimate (EstimateWithError | TailRiskEstimate)
target_relative_error (float)
- Return type:
- prophys.engine.parameter_uncertainty(compiled, observations, *, result=None, raw_params=None, weights=None, prior_neg_log_prob=None, trainable=None, frozen=None)[source]
Standard errors for a fit of compiled to observations.
Pass the
CalibrationResultreturned by the fit as result (or the raw parameter dict as raw_params); the Hessian is taken there. The weights/prior_neg_log_prob/ trainable/frozen arguments must match the ones the fit used, since they define the objective whose curvature is being measured.Parameters excluded by trainable/frozen are held fixed: they were not estimated, so they carry no uncertainty and must not enter the Hessian, where they would contribute rows of an unrelated curvature.
- class prophys.engine.ParameterUncertainty(names, estimates, std_errors, covariance, correlation, n_observations, log_likelihood)[source]
Standard errors and correlations for a fitted parameter vector.
All quantities are in constrained space — the values the model graph uses and the user reads — obtained from the unconstrained-space Hessian by the delta method.
- Parameters:
- names: list[str]
Flattened parameter names, in the order the matrices index them. An array-valued Param contributes one entry per element, named
"route[2,0]".
- log_likelihood: float
Log-likelihood at the optimum (the negated loss, prior term included when one was supplied).
- collinear_pairs(threshold=0.95)[source]
Parameter pairs correlated beyond threshold.
A non-empty result means the data determine a combination of those parameters much better than either one alone. The individual standard errors remain correct, but reporting them as independent uncertainties overstates what the fit resolved.
- prophys.engine.laplace_covariance(loss_fn, raw_params, param_specs, *, n_observations=0)[source]
Observed-information covariance of raw_params under loss_fn.
loss_fn is a negative log-likelihood over unconstrained parameter values — exactly what
prophys.engine.calibrate.negative_log_likelihood()builds — and raw_params must be at its minimum. The Hessian there is the observed Fisher information; its inverse is the parameter covariance.Raises ValueError if the Hessian is not positive definite. That is a statement about the fit, not a numerical inconvenience: at a genuine minimum the Hessian is positive definite, so a failed factorization means the optimizer stopped somewhere that is not one — a saddle, a flat direction, or an unconverged run — and any covariance derived from it would be meaningless rather than merely imprecise.
- prophys.engine.propagate(uncertainty, fn)[source]
Delta-method uncertainty of fn evaluated at the fitted parameters.
fn takes a mapping of flattened parameter name to constrained value and returns a scalar or 1-D array. Its Jacobian at the estimate carries the parameter covariance into the derived quantity:
Cov(g) = J Cov(theta) J^T.This is a first-order propagation: it is exact for a linear fn and an approximation whose quality falls with the curvature of fn over the parameter uncertainty. For a strongly non-linear derived quantity, sample the parameters from the covariance and push the samples through fn instead.
- Parameters:
- Return type:
- class prophys.engine.DerivedUncertainty(values, std_errors, covariance)[source]
Delta-method uncertainty of a quantity derived from the parameters.
- prophys.engine.profile_likelihood(compiled, observations, name, grid, *, result=None, raw_params=None, weights=None, prior_neg_log_prob=None, level=0.95, n_steps=200, learning_rate=0.01)[source]
Profile the likelihood in one parameter over grid.
At each grid value the named parameter is held fixed (in constrained space) and every other parameter is re-optimized; the resulting curve is the profile log-likelihood. The confidence interval is the range where the profile stays within half the chi-square(1) critical value of its maximum.
Unlike the Laplace interval this makes no quadratic assumption, so it reports the asymmetric interval a bounded or skewed parameter actually has. It costs one optimization per grid point.
name must be a scalar Param; profiling one element of an array-valued parameter is not supported, since fixing a single element while re-optimizing the rest of the same array requires splitting the array into separate leaves.
- class prophys.engine.ProfileLikelihood(name, grid, log_likelihood, estimate, max_log_likelihood, level, interval)[source]
A one-parameter profile of the log-likelihood.
- Parameters:
- prophys.engine.sample_parameters(uncertainty, n_draws, seed=0)[source]
Draw parameter vectors from the Laplace (Gaussian) approximation.
Returns a mapping of parameter name to an array of n_draws constrained values, drawn from
N(estimate, covariance). These are the draws a predictive score (prophys.engine.scoring.waic(),prophys.engine.scoring.loo()) integrates over when no MCMC posterior is available.The draws inherit the Laplace approximation’s assumptions: the posterior is taken to be Gaussian in the reported (constrained) parameters, so a draw can fall outside a parameter’s admissible range when the estimate sits within a standard error or two of a bound. For a posterior that respects the constraints exactly, sample with
prophys.inferenceinstead and pass those draws.