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:

.. list-table::
   :header-rows: 1
   :widths: 25 35 40

   * - Source
     - Reduced by
     - Reported by
   * - **Monte-Carlo error** — the statistic was estimated from a finite
       draw
     - More samples
     - :func:`~prophys.engine.assess`, :func:`~prophys.engine.mean_error`,
       :func:`~prophys.engine.quantile_error`,
       :func:`~prophys.engine.cvar_error`
   * - **Parameter uncertainty** — the parameters were estimated from
       finite data
     - More data
     - :func:`~prophys.engine.parameter_uncertainty`,
       :func:`~prophys.engine.propagate`,
       :func:`~prophys.engine.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
------------------

:meth:`CompiledModel.assess <prophys.engine.CompiledModel.assess>` returns
the estimate together with the sampling error of the same draw:

.. code-block:: python

   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 :class:`~prophys.engine.EstimateWithError` (or a
:class:`~prophys.engine.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
  :math:`\sqrt{q(1-q)/n}/f(x_q)` form requires :math:`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
~~~~~~~~~~~~~~~~~~~~~

:func:`~prophys.engine.required_samples` converts an achieved error into the
draws a target precision would take, on the :math:`1/\sqrt{n}` scaling:

.. code-block:: python

   estimate = compiled.assess("loss", "cvar", level=0.99)
   prp.required_samples(estimate, target_relative_error=0.01)

Parameter uncertainty
----------------------

After a fit, :func:`~prophys.engine.parameter_uncertainty` inverts the
observed information — the Hessian of the same negative log-likelihood the
fit minimized — at the optimum:

.. code-block:: python

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

.. code-block:: python

   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,
:meth:`~prophys.engine.ParameterUncertainty.collinear_pairs` lists the
parameter pairs the data constrain only in combination:

.. code-block:: python

   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. :func:`~prophys.engine.propagate` carries the
covariance through the Jacobian of any function of the parameters:

.. code-block:: python

   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
:func:`~prophys.engine.sample_parameters` and push the draws through
instead.

Profile likelihood
~~~~~~~~~~~~~~~~~~~

The Laplace interval assumes the log-likelihood is quadratic near its
maximum. :func:`~prophys.engine.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.

.. code-block:: python

   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
:doc:`inference`.

.. autofunction:: prophys.engine.assess
.. autoclass:: prophys.engine.EstimateWithError
   :members:
.. autoclass:: prophys.engine.TailRiskEstimate
   :members:
.. autofunction:: prophys.engine.mean_error
.. autofunction:: prophys.engine.quantile_error
.. autofunction:: prophys.engine.cvar_error
.. autofunction:: prophys.engine.required_samples

.. autofunction:: prophys.engine.parameter_uncertainty
.. autoclass:: prophys.engine.ParameterUncertainty
   :members:
.. autofunction:: prophys.engine.laplace_covariance
.. autofunction:: prophys.engine.propagate
.. autoclass:: prophys.engine.DerivedUncertainty
   :members:
.. autofunction:: prophys.engine.profile_likelihood
.. autoclass:: prophys.engine.ProfileLikelihood
   :members:
.. autofunction:: prophys.engine.sample_parameters
.. autoclass:: prophys.engine.Report
   :members: