diff --git a/docs/source/user_guide/concepts/beam.rst b/docs/source/user_guide/concepts/beam.rst index f9185379a..761f72a3e 100644 --- a/docs/source/user_guide/concepts/beam.rst +++ b/docs/source/user_guide/concepts/beam.rst @@ -138,19 +138,36 @@ three dimensions. Each dimension answers a different operational question: **Configurable metrics.** Metrics are supplied as provider objects. BEAM includes providers such as ``RMAEProvider`` (Relative Mean Absolute Error) and ``RCRPSProvider`` -(Relative Continuous Ranked Probability Score) for probabilistic evaluation. You can -implement custom providers to add domain-specific metrics. +(Relative Continuous Ranked Probability Score) for probabilistic evaluation. Interval +metrics such as ``RCSProvider`` and +``RIQDProvider`` help assess quantile calibration and sharpness for symmetric +quantile ranges. You can implement custom providers to add domain-specific metrics. .. code-block:: python from openstef_beam.evaluation import EvaluationConfig, EvaluationPipeline - from openstef_beam.evaluation.metric_providers import RMAEProvider, RCRPSProvider + from openstef_beam.evaluation.metric_providers import ( + RCRPSProvider, + RCSProvider, + RMAEProvider, + RIQDProvider, + ) eval_pipeline = EvaluationPipeline( config=EvaluationConfig(), quantiles=forecaster.quantiles, - window_metric_providers=[RMAEProvider(), RCRPSProvider()], - global_metric_providers=[RMAEProvider(), RCRPSProvider()], + window_metric_providers=[ + RMAEProvider(), + RCRPSProvider(), + RCSProvider(), + RIQDProvider(), + ], + global_metric_providers=[ + RMAEProvider(), + RCRPSProvider(), + RCSProvider(), + RIQDProvider(), + ], ) The output is an :class:`~openstef_beam.evaluation.EvaluationReport` containing per-subset @@ -232,4 +249,4 @@ benchmarks (hundreds of targets, multiple models) remain tractable. - :ref:`concept_models` for the forecasting models that BEAM evaluates. - :ref:`concept_metalearning` for how BEAM results inform model selection decisions. - :doc:`/user_guide/guides/backtesting_tutorial` for a hands-on walkthrough of setting up and running a backtest. - - :doc:`/api/beam` for the full openstef-beam API reference. \ No newline at end of file + - :doc:`/api/beam` for the full openstef-beam API reference. diff --git a/docs/source/user_guide/guides/probabilistic_forecasting.rst b/docs/source/user_guide/guides/probabilistic_forecasting.rst index 2f552dfcc..584abadf2 100644 --- a/docs/source/user_guide/guides/probabilistic_forecasting.rst +++ b/docs/source/user_guide/guides/probabilistic_forecasting.rst @@ -253,6 +253,7 @@ Calibration quality can be assessed by comparing expected vs. observed quantile Key metrics for probabilistic forecast quality include: - **Calibration error**: the difference between expected and observed coverage per quantile +- **Regression Coverage Score (RCS)**: the fraction of actual values inside a prediction interval such as P10-P90 - **Sharpness**: the width of prediction intervals (narrower is better, given proper calibration) - **Pinball loss**: the proper scoring rule for quantile forecasts, penalizing both miscalibration and lack of sharpness @@ -263,4 +264,4 @@ See :doc:`/user_guide/guides/backtesting_tutorial` for how to evaluate forecast - :doc:`/user_guide/guides/forecasting` for the overall forecasting workflow (fitting, predicting, model selection). - :doc:`/user_guide/concepts/models` for understanding how different model types compare. - :doc:`/user_guide/guides/backtesting_tutorial` for evaluating forecast performance systematically. - - :doc:`/user_guide/guides/reliability_fallback` for operational concerns like fallback behavior when data is missing. \ No newline at end of file + - :doc:`/user_guide/guides/reliability_fallback` for operational concerns like fallback behavior when data is missing. diff --git a/packages/openstef-beam/src/openstef_beam/evaluation/evaluation_helper.py b/packages/openstef-beam/src/openstef_beam/evaluation/evaluation_helper.py new file mode 100644 index 000000000..99a479128 --- /dev/null +++ b/packages/openstef-beam/src/openstef_beam/evaluation/evaluation_helper.py @@ -0,0 +1,82 @@ +# SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project +# +# SPDX-License-Identifier: MPL-2.0 + +"""Helper functions shared between evaluation metric providers.""" + +from collections.abc import Callable, Sequence +from typing import Any + +import numpy as np +import numpy.typing as npt + +from openstef_beam.evaluation.models.subset import QuantileMetricsDict +from openstef_core.types import Quantile + +type SymmetricQuantileMetric = Callable[..., float] + + +def compute_symmetric_quantile_metrics( + y_true: npt.NDArray[np.floating], + y_pred: npt.NDArray[np.floating], + quantiles: Sequence[Quantile], + metric_name: str, + metric: SymmetricQuantileMetric, + *, + selected_quantiles: Sequence[Quantile] | None = None, + **metric_kwargs: Any, +) -> QuantileMetricsDict: + """Compute metrics for quantiles that have a complementary quantile. + + For each selected quantile, the helper finds its complementary quantile + (``1 - q``), orders the corresponding predictions as lower and upper bounds, + and computes the supplied interval metric. Median quantiles and quantiles + without a complementary counterpart are skipped. + + Args: + y_true: True values with shape (num_samples,). + y_pred: Predicted values with shape (num_samples, num_quantiles). + quantiles: Quantiles used for prediction, in the same order as y_pred columns. + metric_name: Name under which to store the computed metric. + metric: Callable that computes a metric from true values and interval bounds. + selected_quantiles: Optional subset of quantiles to compute metrics for. + metric_kwargs: Additional keyword arguments passed to the metric callable. + + Returns: + QuantileMetricsDict containing metric values for matching quantile pairs. + """ + quantile_indices = {quantile: index for index, quantile in enumerate(quantiles)} + metrics: QuantileMetricsDict = {} + + for quantile, quantile_index in quantile_indices.items(): + if selected_quantiles is not None and quantile not in selected_quantiles: + continue + + complementary_quantile = quantile.complementary() + if quantile == complementary_quantile: + continue + + complementary_index = quantile_indices.get(complementary_quantile) + if complementary_index is None: + continue + + if quantile < complementary_quantile: + lower_pred = y_pred[:, quantile_index] + upper_pred = y_pred[:, complementary_index] + else: + lower_pred = y_pred[:, complementary_index] + upper_pred = y_pred[:, quantile_index] + + metrics[quantile] = { + metric_name: metric( + y_true=y_true, + y_pred_lower_q=lower_pred, + y_pred_upper_q=upper_pred, + **metric_kwargs, + ) + } + + return metrics + + +__all__ = ["compute_symmetric_quantile_metrics"] diff --git a/packages/openstef-beam/src/openstef_beam/evaluation/metric_providers.py b/packages/openstef-beam/src/openstef_beam/evaluation/metric_providers.py index 384dd6e54..3ca28af47 100644 --- a/packages/openstef-beam/src/openstef_beam/evaluation/metric_providers.py +++ b/packages/openstef-beam/src/openstef_beam/evaluation/metric_providers.py @@ -16,6 +16,7 @@ import pandas as pd from pydantic import Field +from openstef_beam.evaluation.evaluation_helper import compute_symmetric_quantile_metrics from openstef_beam.evaluation.models.subset import MetricsDict, QuantileMetricsDict from openstef_beam.metrics import ( completeness, @@ -28,6 +29,7 @@ precision_recall, r2, rcrps, + rcs, relative_pinball_loss, riqd, rmae, @@ -629,8 +631,6 @@ class RIQDProvider(MetricProvider): def metric_names(self) -> frozenset[str]: return frozenset({"rIQD"}) - median_quantile: Quantile = Quantile(0.5) - measurement_range_lower_q: Quantile = Field( default=Quantile(0.05), description="Lower quantile bound for measurement range normalization.", @@ -661,42 +661,60 @@ def compute_probabilistic( Returns: QuantileMetricsDict containing rIQD metrics for each processable quantile. """ - metrics: QuantileMetricsDict = {} - - for i, quantile in enumerate(quantiles): - if self.quantiles is not None and quantile not in self.quantiles: - continue + return compute_symmetric_quantile_metrics( + y_true=y_true, + y_pred=y_pred, + quantiles=quantiles, + selected_quantiles=self.quantiles, + metric_name="rIQD", + metric=riqd, + measurement_range_lower_q=self.measurement_range_lower_q, + measurement_range_upper_q=self.measurement_range_upper_q, + ) - symmetric_quantile = 1.0 - quantile - if np.isclose(quantile, symmetric_quantile, atol=1e-6): - continue # skip if same quantile (e.g., 0.5) +class RCSProvider(MetricProvider): + """Provides Regression Coverage Score metrics. - symmetric_indices = np.nonzero(np.isclose(quantiles, symmetric_quantile, atol=1e-6))[0] + Measures the fraction of observed values inside symmetric prediction + intervals, such as P10-P90. For each quantile, finds its symmetric + counterpart and computes RCS between them. + """ - if len(symmetric_indices) == 0: - continue # no symmetric quantile found, skip + @property + @override + def metric_names(self) -> frozenset[str]: + return frozenset({"RCS"}) - symmetric_idx = symmetric_indices[0] + @override + def compute_probabilistic( + self, + y_true: npt.NDArray[np.floating], + y_pred: npt.NDArray[np.floating], + quantiles: list[Quantile], + ) -> QuantileMetricsDict: + """Compute RCS for each quantile by finding its symmetric counterpart. - if quantile < self.median_quantile: - lower_pred = y_pred[:, i] - upper_pred = y_pred[:, symmetric_idx] - else: - lower_pred = y_pred[:, symmetric_idx] - upper_pred = y_pred[:, i] + For each quantile q, finds the symmetric quantile (1-q) and computes + RCS between them. Only processes quantiles for which a symmetric + counterpart is available. - metrics[quantile] = { - "rIQD": riqd( - y_true=y_true, - y_pred_lower_q=lower_pred, - y_pred_upper_q=upper_pred, - measurement_range_lower_q=self.measurement_range_lower_q, - measurement_range_upper_q=self.measurement_range_upper_q, - ) - } + Args: + y_true: True values, 1D array of shape (num_samples,). + y_pred: Predicted values, 2D array of shape (num_samples, num_quantiles). + quantiles: Quantiles used for prediction, sequence of length (num_quantiles,). - return metrics + Returns: + QuantileMetricsDict containing RCS metrics for each processable quantile. + """ + return compute_symmetric_quantile_metrics( + y_true=y_true, + y_pred=y_pred, + quantiles=quantiles, + selected_quantiles=self.quantiles, + metric_name="RCS", + metric=rcs, + ) class RelativePinballLossProvider(MetricProvider): @@ -748,6 +766,7 @@ def compute_deterministic( "PeakMetricProvider", "R2Provider", "RCRPSProvider", + "RCSProvider", "RIQDProvider", "RMAEPeakHoursProvider", "RMAEProvider", diff --git a/packages/openstef-beam/src/openstef_beam/metrics/__init__.py b/packages/openstef-beam/src/openstef_beam/metrics/__init__.py index 54b8728f1..d0691a5c3 100644 --- a/packages/openstef-beam/src/openstef_beam/metrics/__init__.py +++ b/packages/openstef-beam/src/openstef_beam/metrics/__init__.py @@ -28,6 +28,7 @@ pinball_loss, precision_recall, r2, + rcs, relative_pinball_loss, riqd, rmae, @@ -54,6 +55,7 @@ "precision_recall", "r2", "rcrps", + "rcs", "relative_pinball_loss", "riqd", "rmae", diff --git a/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py b/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py index 33040074d..defcc6040 100644 --- a/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py +++ b/packages/openstef-beam/src/openstef_beam/metrics/metrics_deterministic.py @@ -504,6 +504,44 @@ def riqd( return float(riqd) +def rcs( + y_true: npt.NDArray[np.floating], + y_pred_lower_q: npt.NDArray[np.floating], + y_pred_upper_q: npt.NDArray[np.floating], +) -> float: + """Calculate the Regression Coverage Score (RCS) for prediction intervals. + + RCS measures the fraction of observed values that fall within the predicted + lower and upper quantile bounds. A calibrated 90% prediction interval should + have an RCS close to 0.9. + + Args: + y_true: Observed values with shape (num_samples,). + y_pred_lower_q: Predicted values of lower quantile with shape (num_samples,). + y_pred_upper_q: Predicted values of upper quantile with shape (num_samples,). + + Returns: + The fraction of observations inside the interval. Values closer to the + nominal interval coverage indicate better calibration. + + Example: + Evaluate coverage of a P10-P90 prediction interval + + >>> import numpy as np + >>> y_true = np.array([100, 120, 110, 130]) + >>> y_pred_lower_q = np.array([90, 115, 105, 125]) + >>> y_pred_upper_q = np.array([110, 125, 115, 128]) + >>> rcs(y_true, y_pred_lower_q, y_pred_upper_q) + 0.75 + + Note: + The interval bounds are inclusive: observations equal to either bound + are counted as covered. + """ + is_in_interval = (y_true >= y_pred_lower_q) & (y_true <= y_pred_upper_q) + return float(np.mean(is_in_interval)) + + def r2( y_true: npt.NDArray[np.floating], y_pred: npt.NDArray[np.floating], diff --git a/packages/openstef-beam/tests/unit/evaluation/test_evaluation_helper.py b/packages/openstef-beam/tests/unit/evaluation/test_evaluation_helper.py new file mode 100644 index 000000000..bd0cdf534 --- /dev/null +++ b/packages/openstef-beam/tests/unit/evaluation/test_evaluation_helper.py @@ -0,0 +1,95 @@ +# SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project +# +# SPDX-License-Identifier: MPL-2.0 + +import numpy as np + +from openstef_beam.evaluation.evaluation_helper import compute_symmetric_quantile_metrics +from openstef_core.types import Quantile + + +def test_compute_symmetric_quantile_metrics_uses_complementary_quantiles() -> None: + y_true = np.array([10.0, 20.0, 30.0]) + y_pred = np.array( + [ + [1.0, 5.0, 9.0], + [2.0, 6.0, 10.0], + [3.0, 7.0, 11.0], + ] + ) + quantiles = [Quantile(0.1), Quantile(0.5), Quantile(0.9)] + + def metric(y_true: np.ndarray, y_pred_lower_q: np.ndarray, y_pred_upper_q: np.ndarray) -> float: + assert y_true.tolist() == [10.0, 20.0, 30.0] + return float(np.mean(y_pred_upper_q - y_pred_lower_q)) + + result = compute_symmetric_quantile_metrics( + y_true=y_true, + y_pred=y_pred, + quantiles=quantiles, + metric_name="spread", + metric=metric, + ) + + assert result == { + Quantile(0.1): {"spread": 8.0}, + Quantile(0.9): {"spread": 8.0}, + } + + +def test_compute_symmetric_quantile_metrics_skips_missing_pairs_and_selected_quantiles() -> None: + y_true = np.array([10.0, 20.0, 30.0]) + y_pred = np.array( + [ + [1.0, 3.0, 7.0], + [2.0, 4.0, 8.0], + [3.0, 5.0, 9.0], + ] + ) + quantiles = [Quantile(0.1), Quantile(0.3), Quantile(0.7)] + + result = compute_symmetric_quantile_metrics( + y_true=y_true, + y_pred=y_pred, + quantiles=quantiles, + selected_quantiles=[Quantile(0.7)], + metric_name="spread", + metric=lambda y_true, y_pred_lower_q, y_pred_upper_q: float(np.mean(y_pred_upper_q - y_pred_lower_q)), + ) + + assert result == {Quantile(0.7): {"spread": 4.0}} + + +def test_compute_symmetric_quantile_metrics_passes_metric_kwargs() -> None: + y_true = np.array([10.0, 20.0, 30.0]) + y_pred = np.array( + [ + [1.0, 9.0], + [2.0, 10.0], + [3.0, 11.0], + ] + ) + quantiles = [Quantile(0.1), Quantile(0.9)] + + def metric( + y_true: np.ndarray, + y_pred_lower_q: np.ndarray, + y_pred_upper_q: np.ndarray, + multiplier: float, + ) -> float: + assert y_true.tolist() == [10.0, 20.0, 30.0] + return float(np.mean(y_pred_upper_q - y_pred_lower_q) * multiplier) + + result = compute_symmetric_quantile_metrics( + y_true=y_true, + y_pred=y_pred, + quantiles=quantiles, + metric_name="weighted_spread", + metric=metric, + multiplier=0.5, + ) + + assert result == { + Quantile(0.1): {"weighted_spread": 4.0}, + Quantile(0.9): {"weighted_spread": 4.0}, + } diff --git a/packages/openstef-beam/tests/unit/evaluation/test_metric_provider.py b/packages/openstef-beam/tests/unit/evaluation/test_metric_provider.py index afd3acca6..a072bbde5 100644 --- a/packages/openstef-beam/tests/unit/evaluation/test_metric_provider.py +++ b/packages/openstef-beam/tests/unit/evaluation/test_metric_provider.py @@ -8,7 +8,11 @@ import pandas as pd import pytest -from openstef_beam.evaluation.metric_providers import RIQDProvider, RMAEPeakHoursProvider +from openstef_beam.evaluation.metric_providers import ( + RCSProvider, + RIQDProvider, + RMAEPeakHoursProvider, +) from openstef_core.datasets import ForecastDataset from openstef_core.types import Quantile @@ -121,3 +125,71 @@ def test_riqd_provider_symmetric_quantile_logic( # Verify that skipped quantiles are not in the results for skipped_quantile in expected_skipped: assert Quantile(skipped_quantile) not in result + + +@pytest.mark.parametrize( + ("quantiles", "expected_pairs", "expected_skipped"), + [ + # Standard symmetric quantiles + ([0.1, 0.5, 0.9], [(0.1, 0.9), (0.9, 0.1)], [0.5]), + # More quantiles with multiple symmetric pairs + ([0.05, 0.25, 0.5, 0.75, 0.95], [(0.05, 0.95), (0.25, 0.75), (0.75, 0.25), (0.95, 0.05)], [0.5]), + # Missing symmetric counterpart + ([0.1, 0.5, 0.8], [], [0.1, 0.5, 0.8]), + # Single quantile (median) + ([0.5], [], [0.5]), + # Asymmetric quantiles + ([0.1, 0.3, 0.7], [(0.3, 0.7), (0.7, 0.3)], [0.1]), + ], + ids=["standard_symmetric", "multiple_pairs", "missing_counterpart", "single_median", "asymmetric_partial"], +) +def test_rcs_provider_symmetric_quantile_logic( + quantiles: list[float], expected_pairs: list[tuple[float, float]], expected_skipped: list[float] +) -> None: + """Test that RCSProvider correctly identifies and processes symmetric quantile pairs.""" + # Arrange + provider = RCSProvider() + + # Create test data + start_time = datetime.fromisoformat("2025-01-01T00:00:00") + times = [start_time + timedelta(hours=i) for i in range(24)] + index = pd.DatetimeIndex(times) + + # Create predictions with specified quantiles + quantile_data = {} + for i, q in enumerate(quantiles): + quantile_data[f"quantile_P{int(q * 100):02d}"] = [i * 10 + j for j in range(24)] + + subset = ForecastDataset( + data=pd.DataFrame( + data={ + **quantile_data, + "horizon": timedelta(hours=24), + "load": range(24), + }, + index=index, + ), + target_column="load", + sample_interval=timedelta(hours=1), + ) + + # Act + with patch("openstef_beam.evaluation.metric_providers.rcs", return_value=0.75) as mock_rcs: + result = provider(subset) + + # Assert + # Check that RCS was called for each expected quantile pair + assert mock_rcs.call_count == len(expected_pairs) + + # Verify that results contain metrics for quantiles with symmetric counterparts + expected_result_quantiles = {Quantile(pair[0]) for pair in expected_pairs} + assert set(result.keys()) == expected_result_quantiles + + # Verify each result contains the RCS metric + for quantile in expected_result_quantiles: + assert "RCS" in result[quantile] + assert result[quantile]["RCS"] == 0.75 + + # Verify that skipped quantiles are not in the results + for skipped_quantile in expected_skipped: + assert Quantile(skipped_quantile) not in result diff --git a/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py b/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py index 9bb1e72e6..884c1cde5 100644 --- a/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py +++ b/packages/openstef-beam/tests/unit/metrics/test_metrics_deterministic.py @@ -16,6 +16,7 @@ pinball_loss, precision_recall, r2, + rcs, relative_pinball_loss, riqd, rmae, @@ -460,6 +461,22 @@ def test_riqd_returns_nan_when_inputs_empty() -> None: assert np.isnan(result) +def test_rcs() -> None: + """Test the regression coverage score with a sample prediction interval.""" + y_true = np.array([100.0, 120.0, 110.0, 130.0]) + y_pred_lower_q = np.array([90.0, 115.0, 105.0, 125.0]) + y_pred_upper_q = np.array([110.0, 125.0, 115.0, 128.0]) + + result = rcs( + y_true=y_true, + y_pred_lower_q=y_pred_lower_q, + y_pred_upper_q=y_pred_upper_q, + ) + + assert isinstance(result, float) + assert result == 0.75 + + @pytest.mark.parametrize( ( "y_true",