feat: add regression coverage score metric

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.
+   - :doc:`/api/beam` for the full openstef-beam API reference.

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.
+   - :doc:`/user_guide/guides/reliability_fallback` for operational concerns like fallback behavior when data is missing.

packages/openstef-beam/src/openstef_beam/evaluation/evaluation_helper.py

@@ -0,0 +1,82 @@
+# SPDX-FileCopyrightText: 2025 Contributors to the OpenSTEF project <openstef@lfenergy.org>
+#
+# 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"]

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",

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",

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],