Model Validation

Contents

Model Validation#

Four questions about a fitted model, four groups of functions in prophys.engine.scoring.

Question

Function

Needs

How well does it fit, penalized for complexity?

information_criteria()

A point estimate

How well will it predict new data?

waic(), loo()

Parameter draws

Is its uncertainty the right size?

calibration(), crps()

Observations

Does it reproduce the features of the data?

posterior_predictive_check()

Observations

report_card() runs every applicable one in a single call and is the natural thing to attach to a ModelPackage.

Information criteria#

ic = prp.information_criteria(compiled, {"y": data}, result=result)
print(ic)

InformationCriteria:
  log-likelihood = -552.62  (2 parameters, 300 observations)
  AIC  = 1109.24
  AICc = 1109.28
  BIC  = 1116.65

The parameter count is the number actually estimated, so trainable and frozen must match the fit. The observation count excludes records marked missing, which contribute no likelihood term.

Prefer AICc to AIC whenever \(n/k\) is below about 40; it converges to AIC as \(n\) grows, so it is never the worse choice. BIC penalizes parameters more heavily at any realistic sample size and targets the model that generated the data rather than the one that predicts best.

No prior term is included: an information criterion is defined on the likelihood, and adding the regularizer to the value it is computed from would double-count it. Score a regularized or hierarchical fit with WAIC or PSIS-LOO instead.

Predictive scores#

WAIC and PSIS-LOO estimate the expected log predictive density on new data by integrating over a distribution of parameters, penalizing the effective flexibility the data actually saw rather than the parameter count. Both take parameter draws — from the Laplace approximation, or from a posterior:

u = prp.parameter_uncertainty(compiled, {"y": data}, result=result)
draws = prp.sample_parameters(u, 400)

print(prp.waic(compiled, "y", data, draws))
print(prp.loo(compiled, "y", data, draws))

PSIS-LOO additionally reports a per-observation Pareto shape \(k\) from its importance-weight tail fit. Values above 0.7 mark observations whose leave-one-out estimate is not reliable and should be obtained by actually refitting without them:

score = prp.loo(compiled, "y", data, draws)
print(score.unreliable_points)

Comparing models#

compare() ranks models and pairs each difference with the standard error of that difference:

ranking = prp.compare({"simple": score_a, "richer": score_b})
print(ranking)
ranking.is_distinguishable("richer")

Pairing matters. Two models scored on the same observations have correlated errors, so the standard error of their difference is smaller — often much smaller — than the two scores’ standard errors combined. All scores must therefore be over the same observations in the same order, which compare checks.

Calibration#

A model can have an excellent mean and still be useless for risk work by being systematically over- or under-confident. The check is the probability integral transform: if the predictive distribution \(F\) is correct, \(u = F(y)\) for a genuine observation is uniform on \((0, 1)\).

print(prp.calibration(compiled, "y", data, params=result.raw_params))

CalibrationReport:
  y: 300 observations, calibrated, unbiased
  PIT mean = 0.4976 (0.5 when calibrated), variance = 0.0843 (0.0833)
  KS = 0.0351, CvM = 0.0473, AD = 0.3339
  reliability index = 0.1200

The report names the failure mode directly:

  • dispersion — "under-dispersed" (too confident: the truth lands outside the predicted range too often), "over-dispersed" (too vague), or "calibrated", from the PIT variance against its calibrated value of \(1/12\).

  • bias — whether predictions sit systematically above or below the data, from the PIT mean against 0.5.

Three test statistics accompany them. Anderson-Darling weights the tails, where the Kolmogorov-Smirnov and Cramér-von Mises statistics are least sensitive; for a model that exists to answer tail questions it is the one to read. The reliability index is the rank histogram’s total absolute deviation from flat — an effect size with no null distribution attached, so it stays informative at sample sizes where any test rejects.

CRPS#

The continuous ranked probability score is a proper scoring rule for the whole predictive distribution, expressed in the observation’s own units:

\[\mathrm{CRPS}(F, y) = E|X - y| - \tfrac{1}{2} E|X - X'|\]

It stays finite where the observation lands in a region the model gave zero density (unlike a log-likelihood), rewards sharpness and calibration jointly, and reduces to absolute error for a deterministic forecast — so a probabilistic model and a point prediction can be compared on one scale.

print(prp.crps(compiled, "y", data, params=result.raw_params))

The default fair=True applies the ensemble-size bias correction; without it, scores computed at different n_ensemble are not comparable, because the uncorrected estimator measures the spread of the ensemble members rather than of the distribution behind them.

Posterior predictive checks#

Simulate replicate datasets from the fitted model and compare any summary statistic of the real data against the replicate distribution:

check = prp.posterior_predictive_check(
    compiled, "y", data, statistic="skew", params=result.raw_params
)
print(check.p_value, check.reproduces_the_data)

The p-value is two-sided and small when the model fails to reproduce the statistic. A large value is not evidence the model is right — only that this statistic does not distinguish it from the data.

Choose a statistic the likelihood does not fit directly. Checking the mean of a model whose mean parameter is free will pass by construction and says nothing; checking the skew of a Gaussian fitted to lognormal data fails immediately, which is the point. Built-in names are in TEST_STATISTICS; any callable on a 1-D array works too.

Passing param_draws makes each replicate use a different parameter draw, so the check propagates parameter uncertainty as well as predictive noise.

Report card#

card = prp.report_card(
    compiled, "y", data, result=result, param_draws=draws
)
for name, report in card.items():
    print(name)
    print(report)

The dict it returns is serializable and is what export_model() records under diagnostics, so the diagnostics travel with the exported distributions.

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

Compute AIC, AICc and BIC for a fit of a compiled model to observations.

The parameter count is the number of scalar values actually estimated, so trainable and frozen must match the fit. The observation count excludes records marked missing, which contribute no likelihood term and so should not inflate a sample size penalty.

No prior term is included. An information criterion is defined on the likelihood, and adding a penalty to the value it is computed from would count the regularization twice. Use waic() or loo() to score a regularized or hierarchical fit.

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) – The fitted unconstrained parameters, used instead of result when given.

  • weights (Mapping[str, float] | None) – Likelihood weight per attribute, as used in the fit.

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

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

Returns:

An InformationCriteria report.

Raises:

ValueError – If neither result nor raw_params is given.

Return type:

InformationCriteria

class prophys.engine.InformationCriteria(log_likelihood, n_parameters, n_observations, aic, aicc, bic)[source]#

Information criteria (AIC, AICc and BIC) of a fitted model, from information_criteria().

Parameters:
log_likelihood: float#

Maximized log likelihood of the fit (no prior term).

n_parameters: int#

Parameters actually estimated (frozen ones do not count).

n_observations: int#

Observations contributing a likelihood term (missing ones do not).

aic: float#

-2 log L + 2k. Lower is better.

aicc: float#

AIC with the small sample correction 2k(k+1)/(n-k-1). Prefer it to AIC whenever n/k is below about 40. It converges to AIC as n grows, so it is never the worse choice. inf when n <= k + 1, where no correction is defined.

bic: float#

-2 log L + k log n (inf without observations). Penalizes parameters more heavily than AIC at any realistic sample size, and targets the model that generated the data rather than the one that predicts best.

prophys.engine.waic(compiled, attribute_name, observations, param_draws)[source]#

Score predictive accuracy with the Widely Applicable Information Criterion (WAIC).

Estimates the expected log predictive density from the pointwise log likelihood across parameter draws, penalized by its variance. Cheaper than loo() and usually close to it. Where the two disagree materially, the Pareto diagnostic of loo() says which observations are responsible.

Parameters:
Returns:

A PredictiveScore with criterion="waic".

Return type:

PredictiveScore

prophys.engine.loo(compiled, attribute_name, observations, param_draws)[source]#

Estimate leave one out predictive accuracy by Pareto smoothed importance sampling (PSIS LOO).

Reports the expected log predictive density like waic(), plus a Pareto shape per observation that flags where the estimate itself is unreliable (see PredictiveScore.unreliable_points).

Parameters:
Returns:

A PredictiveScore with criterion="loo" and pareto_k.

Return type:

PredictiveScore

class prophys.engine.PredictiveScore(criterion, elpd, elpd_se, p_eff, deviance, n_observations, pointwise, pareto_k=None)[source]#

Expected log predictive density estimated by WAIC or PSIS LOO, from waic() or loo().

Parameters:
criterion: str#

"waic" or "loo".

elpd: float#

Estimated expected log pointwise predictive density. Higher is better.

elpd_se: float#

Standard error of elpd across observations.

p_eff: float#

Effective number of parameters.

deviance: float#

-2 elpd, the scale AIC and BIC are reported on. Lower is better.

n_observations: int#

Number of observations scored.

pointwise: list[float]#

Elpd contribution of each observation, from which compare() computes differences.

pareto_k: list[float] | None = None#

the fitted tail shape of the importance weights for each observation.

Type:

PSIS LOO only

property unreliable_points: list[int]#

Indices whose Pareto shape exceeds PARETO_K_THRESHOLD or is not finite. Empty for WAIC, which has no reliability diagnostic per observation.

prophys.engine.compare(scores)[source]#

Rank models by predictive score, with the standard error of each difference.

All scores must be over the same observations in the same order. The contributions of each observation are differenced pairwise, which is what makes the standard error of the difference meaningful. Scores computed on different data would give a number with no interpretation, so the observation counts are checked (only the counts, not the values).

Parameters:

scores (Mapping[str, PredictiveScore]) – Model name to its PredictiveScore. At least two, all from the same criterion.

Returns:

A ModelComparison, best model first.

Raises:

ValueError – If fewer than two scores are given, or they differ in criterion or observation count.

Return type:

ModelComparison

Example

>>> ranking = compare({"linear": loo(m1, "y", obs, d1), "quadratic": loo(m2, "y", obs, d2)})
>>> print(ranking)
class prophys.engine.ModelComparison(criterion, names, elpd, delta, delta_se)[source]#

Ranking of models by a predictive score, with each model’s difference from the best.

Returned by compare().

Parameters:
criterion: str#

The criterion the scores share, "waic" or "loo".

names: list[str]#

Model names, ordered best first.

elpd: dict[str, float]#

Each model’s elpd estimate.

delta: dict[str, float]#

Difference from the best model’s elpd (0 for the best, negative for the rest).

delta_se: dict[str, float]#

Standard error of that difference, from the paired contributions of each observation. Pairing matters. Two models scored on the same observations have correlated errors, so the standard error of the difference is smaller, often much smaller, than the standard errors of the two scores combined.

is_distinguishable(name, n_errors=2.0)[source]#

Whether name scores worse than the best by more than n_errors standard errors of the difference. False when that standard error is zero, as for the best model itself.

Parameters:
Return type:

bool

prophys.engine.calibration(compiled, attribute_name, observations, *, params=None, n_bins=10, n_ensemble=4096, key=None, seed=0)[source]#

Test whether the predictive distribution is calibrated against observations (PIT).

Computes the probability integral transform (PIT) of every exactly observed value and tests it for uniformity. Where the attribute’s distribution has a CDF, the PIT is computed exactly from it. Otherwise it falls back to ensemble ranks against n_ensemble predictive draws, with ties broken by uniform random draws so a discrete predictive distribution does not register as miscalibrated. Censored and missing records are skipped.

Parameters:

  • compiled – The model.

  • attribute_name (str) – The observed attribute.

  • observations (Any) – Its observations.

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

  • n_bins (int) – Bins of the rank histogram.

  • n_ensemble (int) – Predictive draws for the rank fallback.

  • key (Any) – PRNG key. Defaults to PRNGKey(seed).

  • seed (int) – Seed used when key is not given.

Returns:

A CalibrationReport.

Raises:

ValueError – If fewer than two values are exactly observed.

Return type:

CalibrationReport

class prophys.engine.CalibrationReport(attribute, n_observations, ks, cvm, ad, pit_mean, pit_variance, histogram, reliability_index)[source]#

Probability integral transform (PIT) calibration of a predictive distribution.

Returned by calibration(). dispersion and bias give a plain language verdict and the test statistics quantify it.

Parameters:
attribute: str#

The attribute scored.

n_observations: int#

Exactly observed values the PIT was computed from.

ks: float#

Kolmogorov Smirnov statistic against uniformity.

cvm: float#

Cramér von Mises statistic.

ad: float#

Anderson Darling statistic, which weights the tails.

pit_mean: float#

0.5 under calibration. A value away from it indicates bias.

pit_variance: float#

1/12 under calibration. A larger value indicates under dispersion.

histogram: list[float]#

Rank histogram counts.

reliability_index: float#

Total absolute deviation of the histogram from flat, in [0, 2).

property dispersion: str#

"calibrated", "under-dispersed" (too confident) or "over-dispersed" (too vague), from the PIT variance against its calibrated value of 1/12.

The 10% band is a reporting convention, not a test. It keeps a well calibrated model from being labelled by ordinary sampling noise. Read ks and ad for the formal statistics.

property bias: str#

"unbiased", "low" (PIT mean below 0.45, the predictions sit above the data) or "high" (PIT mean above 0.55, the predictions sit below the data).

prophys.engine.crps(compiled, attribute_name, observations, *, params=None, n_ensemble=4096, fair=True, key=None, seed=0)[source]#

Score the predictive distribution with the continuous ranked probability score (CRPS).

Scores the whole predictive distribution against each exactly observed value, in the observation’s units. Lower is better. The predictive distribution is represented by n_ensemble draws, and every observation is scored against the same draws. Censored and missing records are skipped.

Parameters:

  • compiled – The model.

  • attribute_name (str) – The observed attribute.

  • observations (Any) – Its observations.

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

  • n_ensemble (int) – Number of predictive draws.

  • fair (bool) – Apply the bias correction for ensemble size (default). Without it, scores computed at different n_ensemble are not comparable.

  • key (Any) – PRNG key. Defaults to PRNGKey(seed).

  • seed (int) – Seed used when key is not given.

Returns:

A CrpsScore.

Raises:

ValueError – If no value is exactly observed.

Return type:

CrpsScore

class prophys.engine.CrpsScore(attribute, mean, std_error, n_observations, n_ensemble, pointwise)[source]#

Continuous ranked probability score (CRPS) of a predictive distribution, from crps().

Parameters:
attribute: str#

The attribute scored.

mean: float#

Mean CRPS over observations, in the attribute’s own units. Lower is better, and zero only for a perfect deterministic forecast.

std_error: float#

Standard error of that mean across observations.

n_observations: int#

Exactly observed values scored.

n_ensemble: int#

Predictive draws the score was estimated from, shared by every observation.

pointwise: list[float]#

CRPS of each observation.

prophys.engine.posterior_predictive_check(compiled, attribute_name, observations, *, statistic='std', params=None, param_draws=None, n_replicates=500, key=None, seed=0)[source]#

Compare a summary statistic of the data against replicate datasets simulated from the model.

Each replicate is a dataset of the same size as observations. When param_draws is supplied, the replicates cycle through the parameter draws, so the check propagates parameter uncertainty as well as predictive noise. With only params, replicates vary by predictive noise alone and the check is correspondingly stricter.

statistic is a name from TEST_STATISTICS or any callable on a one dimensional array. Choose one the likelihood does not fit directly. A check on the mean of a model fitted by maximizing a likelihood whose mean parameter is free passes by construction and says nothing.

The returned p value is two sided and small when the model fails to reproduce the statistic (see PosteriorPredictiveCheck.p_value).

Parameters:

  • compiled – The model.

  • attribute_name (str) – The observed attribute.

  • observations (Any) – Its observations. Only exactly observed values are used.

  • statistic (str | Callable[[ndarray], float]) – "mean", "std", "min", "max", "median", "skew", "range", or a callable.

  • params (Mapping[str, Any] | None) – Unconstrained parameter values, used only when param_draws is not given.

  • param_draws (Mapping[str, Any] | None) – Constrained parameter draws (as for waic()), cycled over the replicates.

  • n_replicates (int) – Number of replicate datasets.

  • key (Any) – PRNG key. Defaults to PRNGKey(seed).

  • seed (int) – Seed used when key is not given.

Returns:

A PosteriorPredictiveCheck.

Raises:

  • ValueError – If fewer than two values are exactly observed.

  • KeyError – If statistic is an unknown name.

Return type:

PosteriorPredictiveCheck

class prophys.engine.PosteriorPredictiveCheck(attribute, statistic, observed, replicate_mean, replicate_std, p_value, n_replicates)[source]#

Posterior predictive check of one summary statistic, from posterior_predictive_check().

Parameters:
attribute: str#

The attribute checked.

statistic: str#

Name of the summary statistic (a callable’s __name__).

observed: float#

The statistic computed on the real observations.

replicate_mean: float#

Mean of the statistic over the replicate datasets.

replicate_std: float#

Standard deviation of the statistic over the replicate datasets.

p_value: float#

twice the smaller tail probability of the observed statistic under the replicate distribution.

Small values mean the model does not reproduce this feature of the data. 1.0 means the observed value sits at the replicate median. A large value is not evidence that the model is right, only that this particular statistic does not distinguish it from the data.

Type:

Two sided posterior predictive p value

n_replicates: int#

Number of replicate datasets simulated.

property reproduces_the_data: bool#

Whether p_value >= 0.05, that is the observed statistic falls inside the central 95% of the replicate distribution.

prophys.engine.report_card(compiled, attribute_name, observations, *, result=None, raw_params=None, param_draws=None, statistics=('std', 'max'), n_ensemble=4096, seed=0)[source]#

Run every applicable goodness of fit diagnostic for one fitted attribute.

Intended as the default thing to run after a fit and to attach to a ModelPackage (pass it as diagnostics to export_model()). The individual functions give finer control. The information criteria count every model parameter as estimated, so call information_criteria() directly for a fit that froze some.

Parameters:

  • compiled – The fitted model.

  • attribute_name (str) – The observed attribute.

  • observations (Any) – Its observations.

  • result – The fit’s CalibrationResult.

  • raw_params (Mapping[str, Any] | None) – The fitted unconstrained parameters, used instead of result when given.

  • param_draws (Mapping[str, Any] | None) – Optional constrained parameter draws. They enable WAIC and PSIS LOO and make the predictive checks integrate over them.

  • statistics (Sequence[str]) – Statistics for the posterior predictive checks.

  • n_ensemble (int) – Predictive draws for calibration and CRPS.

  • seed (int) – PRNG seed for every Monte Carlo step.

Returns:

"information_criteria" (when fitted parameters are given), "calibration", "crps", one "ppc_<statistic>" per entry of statistics, and "waic" and "loo" when param_draws is given.

Return type:

Report name to report

prophys.engine.pointwise_log_likelihood(compiled, attribute_name, observations, param_draws)[source]#

Log likelihood of every observation under every parameter draw.

This is the [n_draws, n_obs] matrix that waic() and loo() are computed from.

Parameters:

  • compiled – The model.

  • attribute_name (str) – The observed attribute.

  • observations (Any) – Its observations (an array or Observations record).

  • param_draws (Mapping[str, Any]) – Parameter name to an array of n_draws constrained scalar values, as returned by sample_parameters() or found in Posterior.draws. Draws outside a bounded parameter’s range are clamped just inside it, and draws of a softplus parameter are clamped to at least 1e-6.

Returns:

NumPy array of shape (n_draws, n_obs).

Raises:

ValueError – If param_draws is empty or its arrays differ in length.

Return type:

ndarray

Kernel correctness#

The native kernels themselves are validated independently of any model:

  • Against independent references, for values and gradients: scipy, OpenTURNS, shapely, scikit-learn, networkx, arviz, SALib, scoringrules and the rainflow package, closed forms evaluated in mpmath at 40–50 significant digits, complex-step differentiation of complex-capable references, and automatic differentiation of independently written JAX formulas.

  • CPU against GPU, bit for bit. Every GPU kernel runs on an emulated device (several threads, in a shuffled order) and must reproduce the CPU kernel’s output exactly, at both precisions, across batch and broadcast layouts; a parity suite repeats the comparison on CUDA hardware.

Contents