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]

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

The parameter count is the number actually estimated, so trainable/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 term to the value it is computed from would double-count the regularization. Use waic() or loo() to score a regularized or hierarchical fit.

Parameters:
Return type:

InformationCriteria

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

AIC, AICc and BIC for a fit.

Parameters:
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 this 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. 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]

Widely Applicable Information Criterion over param_draws.

Estimates the expected log predictive density from the pointwise log-likelihood’s mean and variance across draws. Cheaper than loo() and usually close to it; where the two disagree materially, loo’s Pareto diagnostic says which observations are responsible.

Parameters:
Return type:

PredictiveScore

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

Leave-one-out cross-validation by Pareto-smoothed importance sampling.

Reports the same elpd as waic() plus a per-observation Pareto shape; see PredictiveScore.unreliable_points.

Parameters:
Return type:

PredictiveScore

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

An expected-log-predictive-density estimate: WAIC or PSIS-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.

pointwise: list[float]

Per-observation elpd contributions, which model comparison differences are computed from.

pareto_k: list[float] | None = None

the fitted importance-weight tail shape per observation.

Type:

PSIS-LOO only

property unreliable_points: list[int]

Indices whose Pareto shape exceeds PARETO_K_THRESHOLD (or could not be fitted at all). Empty for WAIC, which has no per-observation reliability diagnostic.

prophys.engine.compare(scores)[source]

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

All scores must be over the same observations in the same order; the per-observation contributions are differenced pairwise, which is what makes the standard error of the difference meaningful. Comparing scores computed on different data would produce a number with no interpretation, so the observation counts are checked.

Parameters:

scores (Mapping[str, PredictiveScore])

Return type:

ModelComparison

class prophys.engine.ModelComparison(criterion, names, elpd, delta, delta_se)[source]

A ranking of models by a predictive score, with pairwise differences against the best.

Parameters:
names: list[str]

Model names, ordered best first.

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 per-observation contributions. 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.

Parameters:
Return type:

bool

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

PIT calibration of the predictive distribution against observations.

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 fresh uniform draws so a discrete predictive distribution does not register as miscalibrated.

Parameters:
Return type:

CalibrationReport

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

Probability-integral-transform calibration of a predictive distribution.

Parameters:
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; away from it indicates bias.

pit_variance: float

1/12 under calibration; above 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/ad for the formal statistics.

property bias: str

"unbiased", "low" (predictions sit above the data) or "high", from the PIT mean against 0.5.

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

Continuous ranked probability score against observations.

Scores the whole predictive distribution, in the observation’s units. fair=True (the default) applies the ensemble-size bias correction, without which scores computed at different n_ensemble are not comparable.

Parameters:
Return type:

CrpsScore

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

Continuous ranked probability score for a predictive distribution.

Parameters:
mean: float

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

std_error: float

Standard error of that mean across observations.

n_ensemble: int

Draws per observation the score was estimated from.

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 fitted model.

Each replicate is a dataset of the same size as observations. When param_draws is supplied, each replicate uses a different parameter draw, 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 1-D 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 will pass 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:
Return type:

PosteriorPredictiveCheck

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

A summary statistic of the data against its replicate distribution.

Parameters:
observed: float

The statistic computed on the real observations.

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

property reproduces_the_data: bool

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

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]

Every applicable diagnostic for one fitted attribute, in one call.

Returns a dict of report name to report: information criteria, calibration, CRPS, one posterior predictive check per entry in statistics, and — when param_draws is supplied — WAIC and PSIS-LOO. Intended as the default thing to run after a fit and to attach to a ModelPackage; the individual functions give finer control.

Parameters:
Return type:

dict[str, Report]

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

[n_draws, n_obs] log-likelihood matrix, the input both predictive scores are computed from.

param_draws maps parameter name to an array of n_draws constrained values — what prophys.engine.uncertainty.sample_parameters() returns, or a posterior from prophys.inference.

Parameters:
Return type:

ndarray