GammaResidual¶

class GammaResidual(BaseConformityScorer):
    r"""Gamma residual scorer using relative prediction errors.

    Computes conformity scores as the signed relative error, normalised
    by the predicted value:

    $$s = \frac{y - \hat{y}}{\hat{y} + \epsilon}$$

    This scorer is useful when the scale of the target variable varies
    over time, because the conformity scores are relative to the prediction
    magnitude. The ``epsilon`` parameter prevents division by zero when
    predictions are near zero.

    Parameters
    ----------
    epsilon : float, default=1e-8
        Small constant added to the denominator to prevent division by
        zero.

    See Also
    --------
    - [`AbsoluteGammaResidual`][yohou.metrics.conformity.AbsoluteGammaResidual] : Symmetric variant using absolute relative errors.
    - [`Residual`][yohou.metrics.conformity.Residual] : Scale-dependent signed residual scorer.

    Examples
    --------
    >>> import polars as pl
    >>> from datetime import date
    >>> from yohou.metrics.conformity import GammaResidual
    >>> scorer = GammaResidual(epsilon=1e-8).fit(
    ...     pl.DataFrame({"time": [date(2020, 1, 1), date(2020, 1, 2)], "y": [1.0, 2.0]})
    ... )
    >>> y_truth = pl.DataFrame({"time": [date(2020, 1, 3)], "y": [10.0]})
    >>> y_pred = pl.DataFrame({"time": [date(2020, 1, 3)], "y": [8.0]})
    >>> scores = scorer.score(y_truth, y_pred)
    >>> round(scores.drop("time").to_series().item(), 4)
    0.25

    """

    _parameter_constraints: dict = {
        **BaseConformityScorer._parameter_constraints,
        "epsilon": [Interval(numbers.Real, 0, None, closed="neither")],
    }

    def __init__(
        self,
        epsilon: float = 1e-8,
        groups: list[str] | dict[str, float] | None = None,
        components: list[str] | dict[str, float] | None = None,
    ) -> None:
        super().__init__(
            groups=groups,
            components=components,
        )
        self.epsilon = epsilon

    def __sklearn_tags__(self):
        """Get the tags for this estimator."""
        tags = super().__sklearn_tags__()
        assert tags.scorer_tags is not None
        tags.scorer_tags.multiplicative = True
        return tags

    def score(self, y_truth: pl.DataFrame, y_pred: pl.DataFrame, /, **score_params) -> pl.DataFrame:
        """Compute gamma (relative) residual conformity scores.

        Parameters
        ----------
        y_truth : pl.DataFrame
            True target values with "time" column.

        y_pred : pl.DataFrame
            Predicted values with "time" column.

        Returns
        -------
        pl.DataFrame
            Relative conformity scores (y_truth - y_pred) / (y_pred + epsilon) with "time" column preserved.

        """
        check_is_fitted(self, ["_is_fitted"])

        # Filter out scorer from score_params to avoid conflict with explicit scorer=self
        score_params_filtered = {k: v for k, v in score_params.items() if k != "scorer"}

        # Validate and align (time dropped, returned as context)
        y_truth, y_pred, context = validate_scorer_data(
            self,
            y_truth,
            y_pred,
            **score_params_filtered,
        )

        # Compute scores and reconstruct with time
        scores_values = (y_truth - y_pred) / (y_pred + self.epsilon)
        scores = pl.DataFrame({"time": context.time_values}).hstack(scores_values)

        return scores

    def inverse_score(
        self, y_pred: pl.DataFrame, conformity_scores: pl.DataFrame, coverage_rate: float
    ) -> pl.DataFrame:
        """Construct prediction intervals from gamma conformity scores.

        Parameters
        ----------
        y_pred : pl.DataFrame
            Point predictions, optionally with "time" column.
        conformity_scores : pl.DataFrame
            Computed conformity scores from calibration set, optionally with "time" column.
        coverage_rate : float
            Desired coverage probability (e.g., 0.9 for 90% intervals).

        Returns
        -------
        pl.DataFrame
            Prediction intervals with lower and upper bounds, and time columns if input had them.

        """
        check_is_fitted(self, ["_is_fitted"])

        # Validate and align
        y_pred, conformity_scores, context = validate_scorer_data(
            self, y_true=None, y_pred=y_pred, scores=conformity_scores, inverse=True
        )

        # Compute quantiles
        lower_q, upper_q = self._compute_asymmetric_quantiles(conformity_scores, coverage_rate)

        # Determine explicit float types for Polars
        lower_q, upper_q = float(lower_q), float(upper_q)

        # Reconstruct y
        denom = y_pred + self.epsilon
        lower_bound = y_pred + lower_q * denom
        upper_bound = y_pred + upper_q * denom

        y_pred_interval = self._format_y_pred_interval(lower_bound, upper_bound, coverage_rate)

        return pl.DataFrame({"time": context.time_values}).hstack(y_pred_interval)