Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 30 additions & 5 deletions deeptab/metrics/distributional.py
Original file line number Diff line number Diff line change
Expand Up @@ -191,12 +191,31 @@ def __repr__(self) -> str:
return f"LogScore(distribution={self._nll.distribution!r})"


def _lognormal_crps(y: np.ndarray, loc: np.ndarray, scale: np.ndarray) -> np.ndarray:
"""Closed-form CRPS of a log-normal predictive distribution (Baran & Lerch, 2015).

``loc`` and ``scale`` are the mean and standard deviation of ``log(X)``.
Observations ``y <= 0`` lie below the support, where the CRPS reduces to
``E[X] - y - E|X - X'| / 2``.
"""
from scipy.stats import norm

mean = np.exp(loc + scale**2 / 2.0)
half_spread = mean * (2.0 * norm.cdf(scale / np.sqrt(2.0)) - 1.0)
positive = y > 0
z = (np.log(np.where(positive, y, 1.0)) - loc) / scale
inside = y * (2.0 * norm.cdf(z) - 1.0) - 2.0 * mean * (norm.cdf(z - scale) + norm.cdf(scale / np.sqrt(2.0)) - 1.0)
below = mean - y - half_spread
return np.where(positive, inside, below)


class CRPS(DeepTabMetric):
"""Continuous Ranked Probability Score (CRPS) for univariate distributions.

Uses vectorised Gaussian ``properscoring`` routines for supported
parameter layouts when available. Otherwise returns a mean absolute
error fallback without estimating predictive spread from observed targets.
Uses the closed-form log-normal CRPS for ``lognormal`` and, when
available, vectorised Gaussian ``properscoring`` routines for the other
supported parameter layouts. Otherwise returns a mean absolute error
fallback without estimating predictive spread from observed targets.

Expected ``y_pred`` format (2-D array, columns are distribution parameters):

Expand All @@ -210,7 +229,9 @@ class CRPS(DeepTabMetric):
With ``properscoring`` installed, the exact Gaussian CRPS is computed for
``normal``. Student-T and Johnson SU use a Gaussian location/scale
approximation, and mixtures use a moment-matched Gaussian. Without that
dependency, all families fall back to mean absolute error.
dependency, these families fall back to mean absolute error. ``lognormal``
always uses the exact closed-form log-normal CRPS, where ``loc`` and
``scale`` are the mean and standard deviation of ``log(y)``.

Notes
-----
Expand Down Expand Up @@ -245,7 +266,11 @@ def __call__(self, y_true: np.ndarray, y_pred: np.ndarray) -> float:
y_true = np.asarray(y_true, dtype=float).ravel()
y_pred = np.asarray(y_pred, dtype=float)

if self.family in ("normal", "lognormal", "studentt", "johnsonsu"):
if self.family == "lognormal":
loc = _col(y_pred, self.loc_col if self.loc_col is not None else 0)
scale = np.clip(_col(y_pred, self.scale_col if self.scale_col is not None else 1), 1e-9, None)
return float(np.mean(_lognormal_crps(y_true, loc, scale)))
if self.family in ("normal", "studentt", "johnsonsu"):
default_columns = {"studentt": (1, 2), "johnsonsu": (2, 3)}
loc_col, scale_col = default_columns.get(self.family, (0, 1))
loc = _col(y_pred, self.loc_col if self.loc_col is not None else loc_col)
Expand Down
20 changes: 11 additions & 9 deletions docs/api/metrics/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -198,14 +198,16 @@ Distributional / LSS Metrics

.. warning::

``CRPS`` uses predicted location and scale for normal and lognormal layouts,
Gaussian approximations for Student's t and Johnson SU, and a moment-matched
Gaussian for mixtures. Other family layouts return mean absolute error of
the first predicted column, even when ``properscoring`` is installed.
That column must contain a predicted mean, not an arbitrary raw family
parameter. Missing ``properscoring`` also produces an MAE fallback for the
supported layouts. MAE does not assess predictive uncertainty; use a
family-specific metric or raw-output NLL when that distinction matters.
``CRPS`` is exact for the normal layout and for the lognormal layout, whose
location and scale are the mean and standard deviation of ``log(y)``. It
uses Gaussian approximations for Student's t and Johnson SU, and a
moment-matched Gaussian for mixtures. Other family layouts return mean
absolute error of the first predicted column, even when ``properscoring``
is installed. That column must contain a predicted mean, not an arbitrary
raw family parameter. Missing ``properscoring`` also produces an MAE
fallback for the Gaussian-based layouts. MAE does not assess predictive
uncertainty; use a family-specific metric or raw-output NLL when that
distinction matters.

.. list-table::
:header-rows: 1
Expand All @@ -230,7 +232,7 @@ Distributional / LSS Metrics
- ``crps``
- ``False``
- ``False``
- Gaussian CRPS for supported layouts; otherwise a documented MAE fallback
- Exact for normal and lognormal; Gaussian approximation for other supported layouts; otherwise a documented MAE fallback
* - :class:`IntervalScore`
- ``interval_score``
- ``False``
Expand Down
17 changes: 17 additions & 0 deletions tests/metrics/test_distributional_metrics.py
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,23 @@ def test_crps_lower_for_better_predictions(self):
bad = np.column_stack([rng.normal(0, 1, 200), np.ones(200) * 2.0])
assert CRPS(family="normal")(y_true, good) < CRPS(family="normal")(y_true, bad)

def test_lognormal_crps_matches_numerical_integration(self):
from scipy import integrate, stats

loc = np.log([100.0, 50.0, 20.0, 5.0])
scale = np.array([0.1, 0.2, 0.3, 0.5])
targets = np.array([105.0, 45.0, 22.0, -1.0])

def reference(target, mean_log, sd_log):
cdf = stats.lognorm(sd_log, scale=np.exp(mean_log)).cdf
above = integrate.quad(lambda x: (1.0 - cdf(x)) ** 2, max(target, 0.0), np.inf, limit=200)[0]
below = integrate.quad(lambda x: cdf(x) ** 2, 0.0, target, limit=200)[0] if target > 0 else -target
return above + below

expected = np.mean([reference(*row) for row in zip(targets, loc, scale, strict=True)])
result = CRPS(family="lognormal")(targets, np.column_stack([loc, scale]))
assert result == pytest.approx(expected, rel=1e-6)

def test_poisson_deviance_nonneg(self, count_data):
y_true, y_pred = count_data
assert PoissonDeviance()(y_true, y_pred) >= 0.0
Expand Down
Loading