Model Validation
==================

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

.. list-table::
   :header-rows: 1
   :widths: 30 30 40

   * - Question
     - Function
     - Needs
   * - How well does it fit, penalized for complexity?
     - :func:`~prophys.engine.information_criteria`
     - A point estimate
   * - How well will it predict new data?
     - :func:`~prophys.engine.waic`, :func:`~prophys.engine.loo`
     - Parameter draws
   * - Is its uncertainty the right size?
     - :func:`~prophys.engine.calibration`, :func:`~prophys.engine.crps`
     - Observations
   * - Does it reproduce the features of the data?
     - :func:`~prophys.engine.posterior_predictive_check`
     - Observations

:func:`~prophys.engine.report_card` runs every applicable one in a single
call and is the natural thing to attach to a
:class:`~prophys.export.ModelPackage`.

Information criteria
---------------------

.. code-block:: python

   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 :math:`n/k` is below about 40; it converges
to AIC as :math:`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:

.. code-block:: python

   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 :math:`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:

.. code-block:: python

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

Comparing models
~~~~~~~~~~~~~~~~~

:func:`~prophys.engine.compare` ranks models and pairs each difference with
*the standard error of that difference*:

.. code-block:: python

   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 :math:`F` is correct,
:math:`u = F(y)` for a genuine observation is uniform on :math:`(0, 1)`.

.. code-block:: python

   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:

* :attr:`~prophys.engine.CalibrationReport.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 :math:`1/12`.
* :attr:`~prophys.engine.CalibrationReport.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:

.. math::

   \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.

.. code-block:: python

   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:

.. code-block:: python

   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
:data:`~prophys.engine.scoring.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
------------

.. code-block:: python

   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
:func:`~prophys.export.export_model` records under ``diagnostics``, so the
diagnostics travel with the exported distributions.

.. autofunction:: prophys.engine.information_criteria
.. autoclass:: prophys.engine.InformationCriteria
   :members:
.. autofunction:: prophys.engine.waic
.. autofunction:: prophys.engine.loo
.. autoclass:: prophys.engine.PredictiveScore
   :members:
.. autofunction:: prophys.engine.compare
.. autoclass:: prophys.engine.ModelComparison
   :members:
.. autofunction:: prophys.engine.calibration
.. autoclass:: prophys.engine.CalibrationReport
   :members:
.. autofunction:: prophys.engine.crps
.. autoclass:: prophys.engine.CrpsScore
   :members:
.. autofunction:: prophys.engine.posterior_predictive_check
.. autoclass:: prophys.engine.PosteriorPredictiveCheck
   :members:
.. autofunction:: prophys.engine.report_card
.. autofunction:: prophys.engine.pointwise_log_likelihood