Uncertainty & Error

Contents

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

assess(), mean_error(), quantile_error(), cvar_error()

Parameter uncertainty — the parameters were estimated from finite data

More data

parameter_uncertainty(), propagate(), profile_likelihood()

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 Process path, 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. The native kernels provide exact first derivatives but no second derivatives of their 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.

The statistic is estimated from a single sample of the attribute, so the Monte Carlo error reported is that of the same draw the estimate came from. Use required_samples() on the result to size a run for a target precision.

Parameters:

  • compiled – A CompiledModel.

  • attribute_name (str) – Which attribute to sample. It must be scalar valued. For a path valued attribute, sample it and pass one step’s draws to mean_error(), quantile_error() or cvar_error().

  • statistic (str) – "mean", "prob", "quantile" or "cvar". "prob" needs level as the threshold and estimates P(attribute <= level) as the mean of an indicator. "quantile" and "cvar" need level as the probability level.

  • level (float | None) – Threshold or probability level, as above.

  • n_samples (int | None) – Draws to use. Defaults to the model’s compiled n_samples.

  • params (Mapping[str, Any] | None) – Unconstrained parameter values by name, such as a fit’s raw_params.

  • key (Any) – 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 variation between runs.

  • n_batches (int) – Batch count for the mean and probability error (see mean_error()).

Returns:

An EstimateWithError, or a TailRiskEstimate for "cvar".

Raises:

ValueError – If statistic is unknown, level is missing where required, or the attribute is not scalar valued.

Return type:

EstimateWithError | TailRiskEstimate

Example

>>> est = assess(compiled, "load", "quantile", level=0.99, n_samples=20000)
>>> est.value, est.interval()
>>> required_samples(est, target_relative_error=0.01)
class prophys.engine.EstimateWithError(statistic, value, std_error, n_samples, effective_sample_size, autocorrelation_time, attribute=None, level=None)[source]#

A Monte Carlo estimate of a mean, probability or quantile with its sampling error.

Returned by assess(), mean_error() and quantile_error(). str() prints the estimate, its 95% interval, relative error, sample size, effective sample size and autocorrelation time.

Parameters:
statistic: str#

Which statistic this is, one of STATISTICS or the name given to mean_error().

value: float#

The estimate itself.

std_error: float#

Monte Carlo standard error of value.

n_samples: int#

Number of finite draws the estimate rests on.

effective_sample_size: float#

Draws’ worth of independent information, n / iact. Equal to n_samples for independent draws, lower for a correlated series. Always n_samples for a quantile, whose draws are taken as independent.

autocorrelation_time: float#

Integrated autocorrelation time (Geyer’s initial positive sequence), 1.0 for independent draws and always 1.0 for a quantile.

attribute: str | None = None#

Model attribute this estimates, when it came from assess().

level: float | None = None#

The q of a quantile estimate or the threshold of a prob estimate. None for a mean.

property relative_error: float#

Standard error as a fraction of the estimate. nan for an estimate that is zero or not finite, where no relative statement is meaningful.

interval(level=0.95)[source]#

Normal theory confidence interval (lower, upper) for the estimate at level.

Parameters:

level (float)

Return type:

tuple[float, float]

is_adequate(target_relative_error=0.05)[source]#

Whether the sample achieves a relative error at or below the target. False when the relative error is undefined.

Parameters:

target_relative_error (float)

Return type:

bool

class prophys.engine.TailRiskEstimate(alpha, var, cvar, std_error, n_tail, n_samples, attribute=None)[source]#

Value at Risk and Conditional Value at Risk estimates with the CVaR’s standard error.

Returned by assess() with statistic="cvar" and by cvar_error().

Parameters:
alpha: float#

Confidence level of the tail.

var: float#

Value at Risk, the alpha quantile.

cvar: float#

Conditional Value at Risk, var + mean((x - var)^+) / (1 - alpha), about the mean of the draws at or above the VaR.

std_error: float#

Influence function standard error of cvar, which includes the uncertainty of the VaR threshold itself.

n_tail: int#

Draws at or above the VaR.

n_samples: int#

Finite draws used in total.

attribute: str | None = None#

Model attribute this estimates, when it came from assess().

property relative_error: float#

Standard error as a fraction of the CVaR. nan for a CVaR that is zero or not finite.

interval(level=0.95)[source]#

Normal theory confidence interval (lower, upper) for the CVaR at level.

Parameters:

level (float)

Return type:

tuple[float, float]

is_adequate(target_relative_error=0.05)[source]#

Whether the CVaR’s relative error is at or below the target. False when undefined.

Parameters:

target_relative_error (float)

Return type:

bool

prophys.engine.mean_error(draws, *, n_batches=0, name='mean')[source]#

Sample mean of draws with its batch means standard error.

Batch means make the error valid for correlated draws (an MCMC chain, a time series) as well as independent ones. Correlation up to a block length raises the reported error, and the effective sample size comes from the integrated autocorrelation time. Non finite draws are dropped, and fewer than four finite draws fall back to sd / sqrt(n).

Parameters:

  • draws (Any) – One dimensional array of draws.

  • n_batches (int) – Number of blocks, clipped to between 2 and n / 2. The default (0) uses about sqrt(n) blocks of about sqrt(n) draws.

  • name (str) – Value of the result’s statistic field.

Returns:

An EstimateWithError.

Raises:

  • ValueError – If draws has more than one axis.

  • RuntimeError – A jax.errors.JaxRuntimeError from the native kernel if no draw is finite.

Return type:

EstimateWithError

prophys.engine.quantile_error(draws, q)[source]#

Sample q quantile of draws with its Maritz and 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. The draws are treated as independent and non finite draws are dropped. The quantile interpolates linearly, as numpy does.

Parameters:

  • draws (Any) – One dimensional array of independent draws.

  • q (float) – Probability level in [0, 1]. At q = 0 the standard error is nan.

Returns:

An EstimateWithError with level=q.

Raises:

  • ValueError – If draws has more than one axis.

  • RuntimeError – A jax.errors.JaxRuntimeError from the native kernel if q is outside [0, 1] or fewer than two draws are finite.

Return type:

EstimateWithError

prophys.engine.cvar_error(draws, alpha)[source]#

Upper tail CVaR of draws at level alpha, with its influence function error.

The CVaR takes the Rockafellar and Uryasev form VaR + mean((x - VaR)^+) / (1 - alpha). Non finite draws are dropped.

Parameters:

  • draws (Any) – One dimensional array of independent draws.

  • alpha (float) – Confidence level in [0, 1). The tail is the draws at or above the alpha quantile.

Returns:

A TailRiskEstimate with the VaR, the CVaR and its standard error.

Raises:

  • ValueError – If draws has more than one axis.

  • RuntimeError – A jax.errors.JaxRuntimeError from the native kernel if alpha is outside [0, 1) or fewer than two draws are finite.

Return type:

TailRiskEstimate

prophys.engine.required_samples(estimate, target_relative_error)[source]#

Number of draws needed to reach a target relative error, extrapolated from an 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:
Returns:

The required sample size, which is below n_samples when the target is already met. The estimate’s own n_samples when its relative error is undefined or zero.

Raises:

ValueError – If target_relative_error is not positive.

Return type:

int

prophys.engine.parameter_uncertainty(compiled, observations, *, result=None, raw_params=None, weights=None, prior_neg_log_prob=None, trainable=None, frozen=None)[source]#

Compute standard errors and correlations of a fit’s parameters (Laplace approximation).

The Hessian of the negative log likelihood the fit minimized is taken at the fitted point in unconstrained space. Its inverse, the covariance, is mapped into constrained space by the delta method (see laplace_covariance()). The weights, prior_neg_log_prob, trainable and frozen arguments must match the ones the fit used, since they define the objective whose curvature is measured.

Parameters excluded by trainable or frozen are held fixed at their fitted values. They were not estimated, so they carry no uncertainty and do not enter the Hessian.

Parameters:

  • compiled – The fitted model.

  • observations (Mapping[str, Any]) – Attribute name to observations, as passed to the fit.

  • result – The fit’s CalibrationResult.

  • raw_params (Mapping[str, Any] | None) – Alternatively, the fitted unconstrained parameters. Takes precedence over result.

  • weights (Mapping[str, float] | None) – As used in the fit.

  • prior_neg_log_prob (Callable[[Mapping[str, Any]], Any] | None) – As used in the fit.

  • trainable (Iterable[str] | None) – As used in the fit.

  • frozen (Iterable[str] | None) – As used in the fit.

Returns:

A ParameterUncertainty report. Array valued parameters contribute one entry per element, named "name[i,j]".

Raises:

ValueError – If neither result nor raw_params is given, no parameter is selected, or the Hessian is not positive definite (the fit has not reached a minimum, or a parameter is not identified by the data).

Return type:

ParameterUncertainty

Example

>>> fit = finetune(compiled, {"conc": measured})
>>> unc = parameter_uncertainty(compiled, {"conc": measured}, result=fit)
>>> unc.interval("emission_rate"), unc.collinear_pairs()
class prophys.engine.ParameterUncertainty(names, estimates, std_errors, covariance, correlation, n_observations, log_likelihood)[source]#

Standard errors, covariance and correlations of fitted parameters (Laplace approximation).

Returned by parameter_uncertainty() and laplace_covariance(). All quantities are in constrained space, the values the model graph uses and the user reads, obtained from the Hessian in unconstrained space 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]".

estimates: dict[str, float]#

Fitted value per name.

std_errors: dict[str, float]#

Standard error per name.

covariance: list[list[float]]#

Covariance matrix over names.

correlation: list[list[float]]#

Correlation matrix over names.

n_observations: int#

Informative observations the fit used, for reference when comparing fits. Whatever the caller passed to laplace_covariance().

log_likelihood: float#

the negated loss, prior term included when one was supplied.

Type:

Log likelihood at the optimum

interval(name, level=0.95)[source]#

Normal theory confidence interval (lower, upper) for one parameter.

Raises:

KeyError – If name is not one of names.

Parameters:
Return type:

tuple[float, float]

intervals(level=0.95)[source]#

Normal theory confidence intervals for every parameter, by name.

Parameters:

level (float)

Return type:

dict[str, tuple[float, float]]

correlation_between(a, b)[source]#

Correlation coefficient between the estimates of parameters a and b.

Raises:

ValueError – If a or b is not one of names.

Parameters:
Return type:

float

collinear_pairs(threshold=0.95)[source]#

Parameter pairs whose absolute correlation is at least threshold.

Returns (a, b, rho) tuples, strongest first. Any result means the data determine a combination of those parameters much better than either alone. The individual standard errors remain correct, but reading them as independent uncertainties overstates what the fit resolved.

Parameters:

threshold (float)

Return type:

list[tuple[str, str, float]]

property relative_errors: dict[str, float]#

Standard error as a fraction of each estimate, nan where the estimate is zero.

prophys.engine.laplace_covariance(loss_fn, raw_params, param_specs, *, n_observations=0)[source]#

Laplace (observed information) covariance of parameters at the minimum of a loss.

The lower level form of parameter_uncertainty(), for a loss built by hand. loss_fn is a negative log likelihood over unconstrained parameter values, as negative_log_likelihood() builds, and raw_params must be at its minimum. The Hessian there (from observed_information()) is the observed Fisher information. Its inverse is the parameter covariance, which is then mapped into constrained space by the delta method.

A Hessian that is not positive definite is a statement about the fit. At a genuine minimum it is positive definite, so a failed factorization means the optimizer stopped elsewhere (a saddle, a flat direction or an unconverged run), and any covariance derived there would be meaningless.

Parameters:

  • loss_fn (Callable[[Mapping[str, Any]], Any]) – Negative log likelihood of a dict of unconstrained values.

  • raw_params (Mapping[str, Any]) – The minimizing unconstrained values, by name. Every entry is treated as estimated.

  • param_specs (Mapping[str, Any]) – The matching Param objects, by name, such as compiled.param_specs.

  • n_observations (int) – Observation count recorded in the result.

Returns:

A ParameterUncertainty in constrained space.

Raises:

ValueError – If raw_params is empty, or the Hessian has non finite entries or is not positive definite.

Return type:

ParameterUncertainty

prophys.engine.propagate(uncertainty, fn)[source]#

Propagate parameter uncertainty into a derived quantity by the delta method.

fn is evaluated at the fitted parameters, and its Jacobian there 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 nonlinear derived quantity, sample the parameters from the covariance (sample_parameters()) and push the samples through fn instead.

Parameters:
Returns:

A DerivedUncertainty with the value, standard error and covariance of every component.

Return type:

DerivedUncertainty

Example

>>> derived = propagate(unc, lambda p: p["a"] * p["b"])
>>> derived.values[0], derived.std_errors[0], derived.interval()
class prophys.engine.DerivedUncertainty(values, std_errors, covariance)[source]#

Delta method standard errors of a quantity derived from the fitted parameters.

Returned by propagate(). A scalar quantity is stored as a single component, index 0.

Parameters:
values: list[float]#

The derived quantity, evaluated at the fitted parameters.

std_errors: list[float]#

Standard error of each component.

covariance: list[list[float]]#

Covariance matrix over the components.

interval(index=0, level=0.95)[source]#

Normal theory confidence interval (lower, upper) for component index.

Parameters:
Return type:

tuple[float, float]

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 log likelihood in one parameter and read off its likelihood ratio interval.

At each grid value the named parameter is frozen at that constrained value and finetune() optimizes every other parameter again, starting from the fitted point. The resulting curve is the profile log likelihood. The confidence interval is the range where the profile stays within half the chi square (1 degree of freedom) critical value of the log likelihood at the fit.

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. Every other parameter is optimized, including any the original fit held frozen.

name must be a scalar Param. Profiling one element of an array valued parameter would need the array split into separate leaves.

Parameters:

  • compiled – The fitted model.

  • observations (Mapping[str, Any]) – Attribute name to observations, as passed to the fit.

  • name (str) – The scalar parameter to profile.

  • grid (Sequence[float]) – Constrained values to fix it at. Should bracket the expected interval on both sides. For a bounded Param, values at or beyond a bound are moved just inside it.

  • result – The fit’s CalibrationResult.

  • raw_params (Mapping[str, Any] | None) – Alternatively, the fitted unconstrained parameters. Takes precedence over result.

  • weights (Mapping[str, float] | None) – As used in the fit.

  • prior_neg_log_prob (Callable[[Mapping[str, Any]], Any] | None) – As used in the fit.

  • level (float) – Confidence level: 0.68, 0.90, 0.95 or 0.99.

  • n_steps (int) – Optimizer steps per grid point.

  • learning_rate (float) – Adam learning rate per grid point.

Returns:

A ProfileLikelihood. A nan interval bound means the grid did not extend far enough on that side.

Raises:

  • ValueError – If level is unsupported, the parameter is not scalar, neither result nor raw_params is given, or a grid value is not positive for a softplus parameter.

  • KeyError – If the model has no parameter name.

Return type:

ProfileLikelihood

class prophys.engine.ProfileLikelihood(name, grid, log_likelihood, estimate, max_log_likelihood, level, interval)[source]#

Profile log likelihood of one parameter with its likelihood ratio interval.

Returned by profile_likelihood().

Parameters:
name: str#

The profiled parameter (constrained space).

grid: list[float]#

Values the parameter was fixed at.

log_likelihood: list[float]#

Log likelihood at each grid point, maximized over all other parameters.

estimate: float#

The parameter’s fitted value (constrained space).

max_log_likelihood: float#

Log likelihood at the fitted parameters.

level: float#

Confidence level of interval.

interval: tuple[float, float]#

where the profile, followed outward from its highest grid point, falls below max_log_likelihood by half the chi square (1 degree of freedom) critical value, linearly interpolated. A nan bound means the grid was too narrow to bracket the interval on that side.

Type:

Likelihood ratio confidence interval

prophys.engine.sample_parameters(uncertainty, n_draws, seed=0)[source]#

Draw parameter vectors from the Laplace (Gaussian) approximation of the fit.

Draws from N(estimate, covariance) of a ParameterUncertainty. These are the draws a predictive score (waic(), 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.inference instead and pass those draws.

Parameters:

  • uncertainty (ParameterUncertainty) – The Laplace approximation.

  • n_draws (int) – Number of draws.

  • seed (int) – NumPy random seed.

Returns:

Flattened parameter name to an array of n_draws constrained values.

Raises:

ValueError – If n_draws is not positive.

Return type:

dict[str, ndarray]

class prophys.engine.Report[source]#

Base class of every diagnostic result in prophys.engine.

Subclasses are frozen dataclasses. On construction every field goes through to_plain(), so a report never holds a live JAX array. str(report) gives a readable console summary, to_dict() a record ready for JSON, and fields() the (name, value) pairs.

to_dict()[source]#

Return this report as a dict ready for JSON serialization.

The dict holds every field plus the report’s class name under the key "report", so a serialized record identifies itself.

Return type:

dict[str, Any]

fields()[source]#

Iterate over the report’s (field name, value) pairs in declaration order.

Return type:

Iterator[tuple[str, Any]]

Contents