Metrics

skforecast.metrics.mean_absolute_scaled_error

mean_absolute_scaled_error(y_true, y_pred, y_train)

Mean Absolute Scaled Error (MASE)

MASE is a scale-independent error metric that measures the accuracy of a forecast. It is the mean absolute error of the forecast divided by the mean absolute error of a naive forecast in the training set. The naive forecast is the one obtained by shifting the time series by one period. If y_train is a list of numpy arrays or pandas Series, it is considered that each element is the true value of the target variable in the training set for each time series. In this case, the naive forecast is calculated for each time series separately.

Parameters:

Name	Type	Description	Default
y_true	pandas Series, numpy ndarray	True values of the target variable.	required
y_pred	pandas Series, numpy ndarray	Predicted values of the target variable.	required
y_train	list, pandas Series, numpy ndarray	True values of the target variable in the training set. If list, it is consider that each element is the true value of the target variable in the training set for each time series.	required

Returns:

Name	Type	Description
mase	float	MASE value.

Examples:

import numpy as np
from skforecast.metrics import mean_absolute_scaled_error

y_true = np.array([3.0, 5.0, 2.5, 7.0])
y_pred = np.array([2.5, 5.5, 2.0, 8.0])
y_train = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
result = mean_absolute_scaled_error(y_true, y_pred, y_train)
print(result)

# 0.625

Source code in skforecast/metrics/metrics.py

def mean_absolute_scaled_error(
    y_true: np.ndarray | pd.Series,
    y_pred: np.ndarray | pd.Series,
    y_train: list[float] | np.ndarray | pd.Series,
) -> float:
    """
    Mean Absolute Scaled Error (MASE)

    MASE is a scale-independent error metric that measures the accuracy of
    a forecast. It is the mean absolute error of the forecast divided by the
    mean absolute error of a naive forecast in the training set. The naive
    forecast is the one obtained by shifting the time series by one period.
    If y_train is a list of numpy arrays or pandas Series, it is considered
    that each element is the true value of the target variable in the training
    set for each time series. In this case, the naive forecast is calculated
    for each time series separately.

    Parameters
    ----------
    y_true : pandas Series, numpy ndarray
        True values of the target variable.
    y_pred : pandas Series, numpy ndarray
        Predicted values of the target variable.
    y_train : list, pandas Series, numpy ndarray
        True values of the target variable in the training set. If `list`, it
        is consider that each element is the true value of the target variable
        in the training set for each time series.

    Returns
    -------
    mase : float
        MASE value.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import mean_absolute_scaled_error

    y_true = np.array([3.0, 5.0, 2.5, 7.0])
    y_pred = np.array([2.5, 5.5, 2.0, 8.0])
    y_train = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
    result = mean_absolute_scaled_error(y_true, y_pred, y_train)
    print(result)

    # 0.625
    ```

    """

    if not isinstance(y_true, (pd.Series, np.ndarray)):
        raise TypeError("`y_true` must be a pandas Series or numpy ndarray.")
    if not isinstance(y_pred, (pd.Series, np.ndarray)):
        raise TypeError("`y_pred` must be a pandas Series or numpy ndarray.")
    if not isinstance(y_train, (list, pd.Series, np.ndarray)):
        raise TypeError("`y_train` must be a list, pandas Series or numpy ndarray.")
    if isinstance(y_train, list):
        for x in y_train:
            if not isinstance(x, (pd.Series, np.ndarray)):
                raise TypeError(
                    "When `y_train` is a list, each element must be a pandas Series "
                    "or numpy ndarray."
                )
    if len(y_true) != len(y_pred):
        raise ValueError("`y_true` and `y_pred` must have the same length.")
    if len(y_true) == 0 or len(y_pred) == 0:
        raise ValueError("`y_true` and `y_pred` must have at least one element.")

    if isinstance(y_train, list):
        naive_forecast = np.concatenate([np.diff(x) for x in y_train])
    else:
        naive_forecast = np.diff(y_train)

    mase = np.mean(np.abs(y_true - y_pred)) / np.nanmean(np.abs(naive_forecast))

    return mase

skforecast.metrics.root_mean_squared_scaled_error

root_mean_squared_scaled_error(y_true, y_pred, y_train)

Root Mean Squared Scaled Error (RMSSE)

RMSSE is a scale-independent error metric that measures the accuracy of a forecast. It is the root mean squared error of the forecast divided by the root mean squared error of a naive forecast in the training set. The naive forecast is the one obtained by shifting the time series by one period. If y_train is a list of numpy arrays or pandas Series, it is considered that each element is the true value of the target variable in the training set for each time series. In this case, the naive forecast is calculated for each time series separately.

Parameters:

Name	Type	Description	Default
y_true	pandas Series, numpy ndarray	True values of the target variable.	required
y_pred	pandas Series, numpy ndarray	Predicted values of the target variable.	required
y_train	list, pandas Series, numpy ndarray	True values of the target variable in the training set. If list, it is consider that each element is the true value of the target variable in the training set for each time series.	required

Returns:

Name	Type	Description
rmsse	float	RMSSE value.

Examples:

import numpy as np
from skforecast.metrics import root_mean_squared_scaled_error

y_true = np.array([3.0, 5.0, 2.5, 7.0])
y_pred = np.array([2.5, 5.5, 2.0, 8.0])
y_train = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
result = root_mean_squared_scaled_error(y_true, y_pred, y_train)
print(result)

# 0.6614378277661477

Source code in skforecast/metrics/metrics.py

def root_mean_squared_scaled_error(
    y_true: np.ndarray | pd.Series,
    y_pred: np.ndarray | pd.Series,
    y_train: list[float] | np.ndarray | pd.Series,
) -> float:
    """
    Root Mean Squared Scaled Error (RMSSE)

    RMSSE is a scale-independent error metric that measures the accuracy of
    a forecast. It is the root mean squared error of the forecast divided by
    the root mean squared error of a naive forecast in the training set. The
    naive forecast is the one obtained by shifting the time series by one period.
    If y_train is a list of numpy arrays or pandas Series, it is considered
    that each element is the true value of the target variable in the training
    set for each time series. In this case, the naive forecast is calculated
    for each time series separately.

    Parameters
    ----------
    y_true : pandas Series, numpy ndarray
        True values of the target variable.
    y_pred : pandas Series, numpy ndarray
        Predicted values of the target variable.
    y_train : list, pandas Series, numpy ndarray
        True values of the target variable in the training set. If list, it
        is consider that each element is the true value of the target variable
        in the training set for each time series.

    Returns
    -------
    rmsse : float
        RMSSE value.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import root_mean_squared_scaled_error

    y_true = np.array([3.0, 5.0, 2.5, 7.0])
    y_pred = np.array([2.5, 5.5, 2.0, 8.0])
    y_train = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
    result = root_mean_squared_scaled_error(y_true, y_pred, y_train)
    print(result)

    # 0.6614378277661477
    ```

    """

    if not isinstance(y_true, (pd.Series, np.ndarray)):
        raise TypeError("`y_true` must be a pandas Series or numpy ndarray.")
    if not isinstance(y_pred, (pd.Series, np.ndarray)):
        raise TypeError("`y_pred` must be a pandas Series or numpy ndarray.")
    if not isinstance(y_train, (list, pd.Series, np.ndarray)):
        raise TypeError("`y_train` must be a list, pandas Series or numpy ndarray.")
    if isinstance(y_train, list):
        for x in y_train:
            if not isinstance(x, (pd.Series, np.ndarray)):
                raise TypeError(
                    "When `y_train` is a list, each element must be a pandas Series "
                    "or numpy ndarray."
                )
    if len(y_true) != len(y_pred):
        raise ValueError("`y_true` and `y_pred` must have the same length.")
    if len(y_true) == 0 or len(y_pred) == 0:
        raise ValueError("`y_true` and `y_pred` must have at least one element.")

    if isinstance(y_train, list):
        naive_forecast = np.concatenate([np.diff(x) for x in y_train])
    else:
        naive_forecast = np.diff(y_train)

    rmsse = np.sqrt(np.mean((y_true - y_pred) ** 2)) / np.sqrt(np.nanmean(naive_forecast ** 2))

    return rmsse

skforecast.metrics.symmetric_mean_absolute_percentage_error

symmetric_mean_absolute_percentage_error(y_true, y_pred)

Compute the Symmetric Mean Absolute Percentage Error (SMAPE).

SMAPE is a relative error metric used to measure the accuracy of forecasts. Unlike MAPE, it is symmetric and prevents division by zero by averaging the absolute values of actual and predicted values.

The result is expressed as a percentage and ranges from 0% (perfect prediction) to 200% (maximum error).

Parameters:

Name	Type	Description	Default
y_true	numpy ndarray, pandas Series	True values of the target variable.	required
y_pred	numpy ndarray, pandas Series	Predicted values of the target variable.	required

Returns:

Name	Type	Description
smape	float	SMAPE value as a percentage.

Notes

When both y_true and y_pred are zero, the corresponding term is treated as zero to avoid division by zero.

Examples:

import numpy as np
from skforecast.metrics import symmetric_mean_absolute_percentage_error

y_true = np.array([100, 200, 0])
y_pred = np.array([110, 180, 10])
result = symmetric_mean_absolute_percentage_error(y_true, y_pred)
print(f"SMAPE: {result:.2f}%")

# SMAPE: 73.35%

Source code in skforecast/metrics/metrics.py

def symmetric_mean_absolute_percentage_error(
    y_true: np.ndarray | pd.Series,
    y_pred: np.ndarray | pd.Series
) -> float:
    """
    Compute the Symmetric Mean Absolute Percentage Error (SMAPE).

    SMAPE is a relative error metric used to measure the accuracy 
    of forecasts. Unlike MAPE, it is symmetric and prevents division 
    by zero by averaging the absolute values of actual and predicted values.

    The result is expressed as a percentage and ranges from 0% 
    (perfect prediction) to 200% (maximum error).

    Parameters
    ----------
    y_true : numpy ndarray, pandas Series
        True values of the target variable.
    y_pred : numpy ndarray, pandas Series
        Predicted values of the target variable.

    Returns
    -------
    smape : float
        SMAPE value as a percentage.

    Notes
    -----
    When both `y_true` and `y_pred` are zero, the corresponding term is treated as zero
    to avoid division by zero.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import symmetric_mean_absolute_percentage_error

    y_true = np.array([100, 200, 0])
    y_pred = np.array([110, 180, 10])
    result = symmetric_mean_absolute_percentage_error(y_true, y_pred)
    print(f"SMAPE: {result:.2f}%")

    # SMAPE: 73.35%
    ```

    """

    if not isinstance(y_true, (pd.Series, np.ndarray)):
        raise TypeError("`y_true` must be a pandas Series or numpy ndarray.")
    if not isinstance(y_pred, (pd.Series, np.ndarray)):
        raise TypeError("`y_pred` must be a pandas Series or numpy ndarray.")
    if len(y_true) != len(y_pred):
        raise ValueError("`y_true` and `y_pred` must have the same length.")
    if len(y_true) == 0 or len(y_pred) == 0:
        raise ValueError("`y_true` and `y_pred` must have at least one element.")

    numerator = np.abs(y_true - y_pred)
    denominator = (np.abs(y_true) + np.abs(y_pred)) / 2

    # NOTE: Avoid division by zero
    mask = denominator != 0
    smape_values = np.zeros_like(denominator)
    smape_values[mask] = numerator[mask] / denominator[mask]

    smape = 100 * np.mean(smape_values)

    return smape

skforecast.metrics.calculate_coverage

calculate_coverage(y_true, lower_bound, upper_bound)

Calculate coverage of a given interval as the proportion of true values that fall within the interval.

Parameters:

Name	Type	Description	Default
y_true	numpy ndarray, pandas Series	True values of the target variable.	required
lower_bound	numpy ndarray, pandas Series	Lower bound of the interval.	required
upper_bound	numpy ndarray, pandas Series	Upper bound of the interval.	required

Returns:

Name	Type	Description
coverage	float	Coverage of the interval.

Examples:

import numpy as np
from skforecast.metrics import calculate_coverage

y_true = np.array([1.0, 2.0, 3.0, 4.0])
lower_bound = np.array([0.5, 1.5, 3.5, 3.0])
upper_bound = np.array([1.5, 2.5, 4.5, 5.0])
result = calculate_coverage(y_true, lower_bound, upper_bound)
print(result)

# 0.75

Source code in skforecast/metrics/metrics.py

def calculate_coverage(
    y_true: np.ndarray | pd.Series,
    lower_bound: np.ndarray | pd.Series,
    upper_bound: np.ndarray | pd.Series,
) -> float:
    """
    Calculate coverage of a given interval as the proportion of true values
    that fall within the interval.

    Parameters
    ----------
    y_true : numpy ndarray, pandas Series
        True values of the target variable.
    lower_bound : numpy ndarray, pandas Series
        Lower bound of the interval.
    upper_bound : numpy ndarray, pandas Series
        Upper bound of the interval.

    Returns
    -------
    coverage : float
        Coverage of the interval.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import calculate_coverage

    y_true = np.array([1.0, 2.0, 3.0, 4.0])
    lower_bound = np.array([0.5, 1.5, 3.5, 3.0])
    upper_bound = np.array([1.5, 2.5, 4.5, 5.0])
    result = calculate_coverage(y_true, lower_bound, upper_bound)
    print(result)

    # 0.75
    ```

    """
    if not isinstance(y_true, (np.ndarray, pd.Series)) or y_true.ndim != 1:
        raise TypeError("`y_true` must be a 1D numpy array or pandas Series.")

    if not isinstance(lower_bound, (np.ndarray, pd.Series)) or lower_bound.ndim != 1:
        raise TypeError("`lower_bound` must be a 1D numpy array or pandas Series.")

    if not isinstance(upper_bound, (np.ndarray, pd.Series)) or upper_bound.ndim != 1:
        raise TypeError("`upper_bound` must be a 1D numpy array or pandas Series.")

    y_true = np.asarray(y_true)
    lower_bound = np.asarray(lower_bound)
    upper_bound = np.asarray(upper_bound)

    if y_true.shape != lower_bound.shape or y_true.shape != upper_bound.shape:
        raise ValueError(
            "`y_true`, `lower_bound` and `upper_bound` must have the same shape."
        )

    coverage = np.mean(np.logical_and(y_true >= lower_bound, y_true <= upper_bound))

    return coverage

skforecast.metrics.crps_from_predictions

crps_from_predictions(y_true, y_pred)

Compute the Continuous Ranked Probability Score (CRPS) for a set of forecast realizations, for example from bootstrapping. The CRPS compares the empirical distribution of a set of forecasted values to a scalar observation. The smaller the CRPS, the better.

Parameters:

Name	Type	Description	Default
y_true	float	The true value of the random variable.	required
y_pred	ndarray	The predicted values of the random variable. These are the multiple forecasted values for a single observation.	required

Returns:

Name	Type	Description
crps	float	The CRPS score.

Examples:

import numpy as np
from skforecast.metrics import crps_from_predictions

y_true = 5.0
y_pred = np.array([4.5, 5.2, 5.0, 4.8, 5.5])
result = crps_from_predictions(y_true, y_pred)
print(result)

# 0.08800000000000005

Source code in skforecast/metrics/metrics.py

def crps_from_predictions(
    y_true: float, 
    y_pred: np.ndarray
) -> float:
    """
    Compute the Continuous Ranked Probability Score (CRPS) for a set of
    forecast realizations, for example from bootstrapping. The CRPS compares
    the empirical distribution of a set of forecasted values to a scalar
    observation. The smaller the CRPS, the better.

    Parameters
    ----------
    y_true : float
        The true value of the random variable.
    y_pred : np.ndarray
        The predicted values of the random variable. These are the multiple
        forecasted values for a single observation.

    Returns
    -------
    crps : float
        The CRPS score.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import crps_from_predictions

    y_true = 5.0
    y_pred = np.array([4.5, 5.2, 5.0, 4.8, 5.5])
    result = crps_from_predictions(y_true, y_pred)
    print(result)

    # 0.08800000000000005
    ```

    """
    if not isinstance(y_pred, np.ndarray) or y_pred.ndim != 1:
        raise TypeError("`y_pred` must be a 1D numpy array.")

    if not isinstance(y_true, (float, int)):
        raise TypeError("`y_true` must be a float or integer.")

    y_pred = np.sort(y_pred)
    # Define the grid for integration including the true value
    grid = np.concatenate(([y_true], y_pred))
    grid = np.sort(grid)
    cdf_values = np.searchsorted(y_pred, grid, side='right') / len(y_pred)
    indicator = grid >= y_true
    diffs = np.diff(grid)
    crps = np.sum(diffs * (cdf_values[:-1] - indicator[:-1])**2)

    return crps

skforecast.metrics.crps_from_quantiles

crps_from_quantiles(
    y_true, pred_quantiles, quantile_levels
)

Calculate the Continuous Ranked Probability Score (CRPS) for a given true value and predicted quantiles. The empirical cdf is approximated using linear interpolation between the predicted quantiles.

Parameters:

Name	Type	Description	Default
y_true	float	The true value of the random variable.	required
pred_quantiles	numpy ndarray	The predicted quantile values.	required
quantile_levels	numpy ndarray	The quantile levels corresponding to the predicted quantiles.	required

Returns:

Name	Type	Description
crps	float	The CRPS score.

Examples:

import numpy as np
from skforecast.metrics import crps_from_quantiles

y_true = 9.5
pred_quantiles = np.array([8.0, 10.0, 12.0])
quantile_levels = np.array([0.1, 0.5, 0.9])
result = crps_from_quantiles(y_true, pred_quantiles, quantile_levels)
print(result)

# 0.4634331094524621

Source code in skforecast/metrics/metrics.py

def crps_from_quantiles(
    y_true: float,
    pred_quantiles: np.ndarray,
    quantile_levels: np.ndarray,
) -> float:
    """
    Calculate the Continuous Ranked Probability Score (CRPS) for a given true value
    and predicted quantiles. The empirical cdf is approximated using linear interpolation
    between the predicted quantiles.

    Parameters
    ----------
    y_true : float
        The true value of the random variable.
    pred_quantiles : numpy ndarray
        The predicted quantile values.
    quantile_levels : numpy ndarray
        The quantile levels corresponding to the predicted quantiles.

    Returns
    -------
    crps : float
        The CRPS score.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import crps_from_quantiles

    y_true = 9.5
    pred_quantiles = np.array([8.0, 10.0, 12.0])
    quantile_levels = np.array([0.1, 0.5, 0.9])
    result = crps_from_quantiles(y_true, pred_quantiles, quantile_levels)
    print(result)

    # 0.4634331094524621
    ```

    """
    if not isinstance(y_true, (float, int)):
        raise TypeError("`y_true` must be a float or integer.")

    if not isinstance(pred_quantiles, np.ndarray) or pred_quantiles.ndim != 1:
        raise TypeError("`pred_quantiles` must be a 1D numpy array.")

    if not isinstance(quantile_levels, np.ndarray) or quantile_levels.ndim != 1:
        raise TypeError("`quantile_levels` must be a 1D numpy array.")

    if len(pred_quantiles) != len(quantile_levels):
        raise ValueError(
            "The number of predicted quantiles and quantile levels must be equal."
        )

    if np.any((quantile_levels < 0) | (quantile_levels > 1)):
        raise ValueError("All quantile levels must be between 0 and 1.")

    sorted_indices = np.argsort(pred_quantiles)
    pred_quantiles = pred_quantiles[sorted_indices]
    quantile_levels = quantile_levels[sorted_indices]

    # Define the empirical CDF function using interpolation
    def empirical_cdf(x):
        return np.interp(x, pred_quantiles, quantile_levels, left=0.0, right=1.0)

    # Define the CRPS integrand
    def crps_integrand(x):
        return (empirical_cdf(x) - (x >= y_true)) ** 2

    xmin = pred_quantiles[0]
    xmax = pred_quantiles[-1]

    # Outside the range of the predicted quantiles the integrand is exactly 0 or 1,
    # so its contribution is computed analytically instead of numerically. This also
    # keeps the score valid when `y_true` falls outside the predicted quantiles.
    tail_area = max(0.0, xmin - y_true) + max(0.0, y_true - xmax)

    if xmax > xmin:
        # Create a fine grid of x values for integration
        x_values = np.linspace(xmin, xmax, 1000)

        # Compute the integrand values and integrate using the trapezoidal rule
        integrand_values = crps_integrand(x_values)
        if hasattr(np, "trapezoid"):
            crps = np.trapezoid(integrand_values, x=x_values)
        else:
            crps = np.trapz(integrand_values, x_values)
        crps = crps + tail_area
    else:
        # All predicted quantiles are identical, the forecast is a point mass
        crps = tail_area

    return float(crps)

skforecast.metrics.winkler_score

winkler_score(y_true, lower_bound, upper_bound, alpha)

Winkler Score (Interval Score) for evaluating prediction intervals.

Penalises both wide intervals and observations outside the interval. A lower score indicates a better interval forecast.

score_i = (upper_i - lower_i)
          + (2/alpha) * max(0, lower_i - y_true_i)
          + (2/alpha) * max(0, y_true_i - upper_i)

Parameters:

Name	Type	Description	Default
y_true	numpy ndarray, pandas Series	True values of the target variable. Shape (n,).	required
lower_bound	numpy ndarray, pandas Series	Lower bound of the prediction interval. Shape (n,).	required
upper_bound	numpy ndarray, pandas Series	Upper bound of the prediction interval. Shape (n,).	required
alpha	float	Significance level (e.g. 0.05 for a 95% interval). Must be in (0, 1).	required

Returns:

Name	Type	Description
score	float	Mean Winkler Score across all observations. Lower is better.

Notes

Strictly proper scoring rule for interval forecasts. Standard metric in the M4 and M5 Forecasting Competitions. The penalty for observations that fall outside the interval is controlled by alpha, so the metric can be tuned to how costly a missed observation is relative to a wide interval [1]_.

References

.. [1] Winkler, R. L. (1972). A Decision-Theoretic Approach to Interval Estimation. Journal of the American Statistical Association, 67(337), 187-191.

.. [2] Gneiting, T., & Raftery, A. E. (2007). Strictly Proper Scoring Rules, Prediction, and Estimation. Journal of the American Statistical Association, 102(477), 359-378.

Examples:

import numpy as np
from skforecast.metrics import winkler_score

y_true = np.array([100., 200., 150., 80.])
lower_bound = np.array([90., 180., 140., 100.])
upper_bound = np.array([110., 220., 160., 120.])
result = winkler_score(y_true, lower_bound, upper_bound, alpha=0.05)
print(result)

# 225.0

Source code in skforecast/metrics/metrics.py

def winkler_score(
    y_true: np.ndarray | pd.Series,
    lower_bound: np.ndarray | pd.Series,
    upper_bound: np.ndarray | pd.Series,
    alpha: float,
) -> float:
    """
    Winkler Score (Interval Score) for evaluating prediction intervals.

    Penalises both wide intervals and observations outside the interval.
    A lower score indicates a better interval forecast.

        score_i = (upper_i - lower_i)
                  + (2/alpha) * max(0, lower_i - y_true_i)
                  + (2/alpha) * max(0, y_true_i - upper_i)

    Parameters
    ----------
    y_true : numpy ndarray, pandas Series
        True values of the target variable. Shape (n,).
    lower_bound : numpy ndarray, pandas Series
        Lower bound of the prediction interval. Shape (n,).
    upper_bound : numpy ndarray, pandas Series
        Upper bound of the prediction interval. Shape (n,).
    alpha : float
        Significance level (e.g. 0.05 for a 95% interval). Must be in (0, 1).

    Returns
    -------
    score : float
        Mean Winkler Score across all observations. Lower is better.

    Notes
    -----
    Strictly proper scoring rule for interval forecasts. Standard metric in
    the M4 and M5 Forecasting Competitions. The penalty for observations that
    fall outside the interval is controlled by `alpha`, so the metric can be
    tuned to how costly a missed observation is relative to a wide interval [1]_.

    References
    ----------
    .. [1] Winkler, R. L. (1972). A Decision-Theoretic Approach to Interval
           Estimation. Journal of the American Statistical Association,
           67(337), 187-191.

    .. [2] Gneiting, T., & Raftery, A. E. (2007). Strictly Proper Scoring Rules,
           Prediction, and Estimation. Journal of the American Statistical
           Association, 102(477), 359-378.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import winkler_score

    y_true = np.array([100., 200., 150., 80.])
    lower_bound = np.array([90., 180., 140., 100.])
    upper_bound = np.array([110., 220., 160., 120.])
    result = winkler_score(y_true, lower_bound, upper_bound, alpha=0.05)
    print(result)

    # 225.0
    ```

    """
    if not isinstance(y_true, (np.ndarray, pd.Series)) or np.asarray(y_true).ndim != 1:
        raise TypeError("`y_true` must be a 1D numpy array or pandas Series.")
    if (
        not isinstance(lower_bound, (np.ndarray, pd.Series))
        or np.asarray(lower_bound).ndim != 1
    ):
        raise TypeError("`lower_bound` must be a 1D numpy array or pandas Series.")
    if (
        not isinstance(upper_bound, (np.ndarray, pd.Series))
        or np.asarray(upper_bound).ndim != 1
    ):
        raise TypeError("`upper_bound` must be a 1D numpy array or pandas Series.")
    if not isinstance(alpha, (float, int)) or not (0 < float(alpha) < 1):
        raise ValueError("`alpha` must be a float strictly between 0 and 1.")

    y_true = np.asarray(y_true, dtype=float)
    lower_bound = np.asarray(lower_bound, dtype=float)
    upper_bound = np.asarray(upper_bound, dtype=float)

    if not (y_true.shape == lower_bound.shape == upper_bound.shape):
        raise ValueError(
            "`y_true`, `lower_bound`, and `upper_bound` must have the same shape."
        )
    if len(y_true) == 0:
        raise ValueError("`y_true` must have at least one element.")
    if np.any(upper_bound < lower_bound):
        raise ValueError(
            "All values in `upper_bound` must be >= corresponding `lower_bound`."
        )

    width = upper_bound - lower_bound
    penalty_lower = (2.0 / alpha) * np.maximum(0.0, lower_bound - y_true)
    penalty_upper = (2.0 / alpha) * np.maximum(0.0, y_true - upper_bound)

    return float(np.mean(width + penalty_lower + penalty_upper))

skforecast.metrics.weighted_interval_score

weighted_interval_score(
    y_true, y_pred, lower_bounds, upper_bounds, alphas
)

Weighted Interval Score (WIS) for evaluating probabilistic forecasts.

Generalises the Winkler Score to K prediction intervals plus a point (median) forecast. Converges to CRPS as K grows. Primary evaluation metric of the US CDC COVID-19 Forecast Hub.

WIS = (1 / (K + 0.5)) * (0.5*|y - m| + sum_k(alpha_k/2 * IS_k))

Parameters:

Name	Type	Description	Default
y_true	numpy ndarray, pandas Series	True values of the target variable. Shape (n,).	required
y_pred	numpy ndarray, pandas Series	Predicted median (0.5 quantile) forecast. Shape (n,).	required
lower_bounds	numpy ndarray, pandas DataFrame	Lower bounds of the K prediction intervals. Shape (n, K). Columns are matched to alphas by position, not by label.	required
upper_bounds	numpy ndarray, pandas DataFrame	Upper bounds of the K prediction intervals. Shape (n, K). Columns are matched to alphas by position, not by label.	required
alphas	list, numpy ndarray	Significance levels of the K intervals. Shape (K,). Values in (0, 1).	required

Returns:

Name	Type	Description
wis	float	Mean Weighted Interval Score across observations. Lower is better.

Notes

Proper scoring rule that decomposes into sharpness (interval width) and calibration (miscoverage penalty), providing a single score to compare full predictive distributions [1]_.

References

.. [1] Bracher, J., Ray, E. L., Gneiting, T., & Reich, N. G. (2021). Evaluating epidemic forecasts in an interval format. PLOS Computational Biology, 17(2), e1008618. https://doi.org/10.1371/journal.pcbi.1008618

Examples:

import numpy as np
from skforecast.metrics import weighted_interval_score

y_true = np.array([100., 200., 150.])
y_pred = np.array([98., 195., 155.])
lower_bounds = np.array([[88., 80.], [175., 165.], [138., 128.]])
upper_bounds = np.array([[108., 118.], [215., 225.], [168., 178.]])
alphas = np.array([0.20, 0.05])
result = weighted_interval_score(
    y_true, y_pred, lower_bounds, upper_bounds, alphas
)
print(result)

# 2.4933333333333336

Source code in skforecast/metrics/metrics.py

def weighted_interval_score(
    y_true: np.ndarray | pd.Series,
    y_pred: np.ndarray | pd.Series,
    lower_bounds: np.ndarray | pd.DataFrame,
    upper_bounds: np.ndarray | pd.DataFrame,
    alphas: list[float] | np.ndarray,
) -> float:
    """
    Weighted Interval Score (WIS) for evaluating probabilistic forecasts.

    Generalises the Winkler Score to K prediction intervals plus a point
    (median) forecast. Converges to CRPS as K grows. Primary evaluation
    metric of the US CDC COVID-19 Forecast Hub.

        WIS = (1 / (K + 0.5)) * (0.5*|y - m| + sum_k(alpha_k/2 * IS_k))

    Parameters
    ----------
    y_true : numpy ndarray, pandas Series
        True values of the target variable. Shape (n,).
    y_pred : numpy ndarray, pandas Series
        Predicted median (0.5 quantile) forecast. Shape (n,).
    lower_bounds : numpy ndarray, pandas DataFrame
        Lower bounds of the K prediction intervals. Shape (n, K). Columns are
        matched to `alphas` by position, not by label.
    upper_bounds : numpy ndarray, pandas DataFrame
        Upper bounds of the K prediction intervals. Shape (n, K). Columns are
        matched to `alphas` by position, not by label.
    alphas : list, numpy ndarray
        Significance levels of the K intervals. Shape (K,). Values in (0, 1).

    Returns
    -------
    wis : float
        Mean Weighted Interval Score across observations. Lower is better.

    Notes
    -----
    Proper scoring rule that decomposes into sharpness (interval width) and
    calibration (miscoverage penalty), providing a single score to compare
    full predictive distributions [1]_.

    References
    ----------
    .. [1] Bracher, J., Ray, E. L., Gneiting, T., & Reich, N. G. (2021).
           Evaluating epidemic forecasts in an interval format.
           PLOS Computational Biology, 17(2), e1008618.
           https://doi.org/10.1371/journal.pcbi.1008618

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import weighted_interval_score

    y_true = np.array([100., 200., 150.])
    y_pred = np.array([98., 195., 155.])
    lower_bounds = np.array([[88., 80.], [175., 165.], [138., 128.]])
    upper_bounds = np.array([[108., 118.], [215., 225.], [168., 178.]])
    alphas = np.array([0.20, 0.05])
    result = weighted_interval_score(
        y_true, y_pred, lower_bounds, upper_bounds, alphas
    )
    print(result)

    # 2.4933333333333336
    ```

    """
    if not isinstance(y_true, (np.ndarray, pd.Series)) or np.asarray(y_true).ndim != 1:
        raise TypeError("`y_true` must be a 1D numpy array or pandas Series.")
    if (
        not isinstance(y_pred, (np.ndarray, pd.Series))
        or np.asarray(y_pred).ndim != 1
    ):
        raise TypeError("`y_pred` must be a 1D numpy array or pandas Series.")

    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    lower_bounds = np.asarray(lower_bounds, dtype=float)
    upper_bounds = np.asarray(upper_bounds, dtype=float)
    alphas = np.asarray(alphas, dtype=float)

    if lower_bounds.ndim != 2:
        raise ValueError(
            "`lower_bounds` must be a 2D array with shape "
            "(n_observations, n_intervals)."
        )
    if upper_bounds.ndim != 2:
        raise ValueError(
            "`upper_bounds` must be a 2D array with shape "
            "(n_observations, n_intervals)."
        )
    if alphas.ndim != 1:
        raise ValueError("`alphas` must be a 1D array of significance levels.")
    if len(y_true) == 0:
        raise ValueError("`y_true` must have at least one element.")

    n_obs = len(y_true)
    n_intervals = len(alphas)

    if len(y_pred) != n_obs:
        raise ValueError("`y_true` and `y_pred` must have the same length.")
    if lower_bounds.shape != (n_obs, n_intervals):
        raise ValueError(
            f"`lower_bounds` must have shape ({n_obs}, {n_intervals}). "
            f"Got {lower_bounds.shape}."
        )
    if upper_bounds.shape != (n_obs, n_intervals):
        raise ValueError(
            f"`upper_bounds` must have shape ({n_obs}, {n_intervals}). "
            f"Got {upper_bounds.shape}."
        )
    if np.any((alphas <= 0) | (alphas >= 1)):
        raise ValueError("All values in `alphas` must be strictly between 0 and 1.")
    if np.any(upper_bounds < lower_bounds):
        raise ValueError(
            "All values in `upper_bounds` must be >= corresponding `lower_bounds`."
        )

    abs_error = np.abs(y_true - y_pred)
    width = upper_bounds - lower_bounds
    penalty_lower = (2.0 / alphas[np.newaxis, :]) * np.maximum(
        0.0, lower_bounds - y_true[:, np.newaxis]
    )
    penalty_upper = (2.0 / alphas[np.newaxis, :]) * np.maximum(
        0.0, y_true[:, np.newaxis] - upper_bounds
    )
    interval_scores = width + penalty_lower + penalty_upper
    weighted_sum = np.sum((alphas[np.newaxis, :] / 2.0) * interval_scores, axis=1)
    wis_per_obs = (1.0 / (n_intervals + 0.5)) * (0.5 * abs_error + weighted_sum)

    return float(np.mean(wis_per_obs))

skforecast.metrics.create_mean_pinball_loss

create_mean_pinball_loss(alpha)

Create pinball loss, also known as quantile loss, for a given quantile. Internally, it uses the mean_pinball_loss function from scikit-learn.

Parameters:

Name	Type	Description	Default
alpha	float	Quantile for which the Pinball loss is calculated. Must be between 0 and 1, inclusive.	required

Returns:

Name	Type	Description
mean_pinball_loss_q	callable	Mean Pinball loss for the given quantile. Its name includes the quantile, mean_pinball_loss_q_{alpha} (for example, mean_pinball_loss_q_0.9), so that the losses of several quantiles are different columns of the backtesting metrics.

Examples:

import numpy as np
from skforecast.metrics import create_mean_pinball_loss

pinball_loss_q90 = create_mean_pinball_loss(alpha=0.9)
y_true = np.array([3.0, 5.0, 2.5, 7.0])
y_pred = np.array([2.5, 5.5, 2.0, 8.0])
result = pinball_loss_q90(y_true, y_pred)
print(result)

# 0.26249999999999996

Source code in skforecast/metrics/metrics.py

def create_mean_pinball_loss(alpha: float) -> Callable:
    """
    Create pinball loss, also known as quantile loss, for a given quantile.
    Internally, it uses the `mean_pinball_loss` function from scikit-learn.

    Parameters
    ----------
    alpha: float
        Quantile for which the Pinball loss is calculated. Must be between 0 and 1, inclusive.

    Returns
    -------
    mean_pinball_loss_q: callable
        Mean Pinball loss for the given quantile. Its name includes the
        quantile, `mean_pinball_loss_q_{alpha}` (for example,
        `mean_pinball_loss_q_0.9`), so that the losses of several quantiles
        are different columns of the backtesting metrics.

    Examples
    --------
    ```python
    import numpy as np
    from skforecast.metrics import create_mean_pinball_loss

    pinball_loss_q90 = create_mean_pinball_loss(alpha=0.9)
    y_true = np.array([3.0, 5.0, 2.5, 7.0])
    y_pred = np.array([2.5, 5.5, 2.0, 8.0])
    result = pinball_loss_q90(y_true, y_pred)
    print(result)

    # 0.26249999999999996
    ```

    """
    if not (0 <= alpha <= 1):
        raise ValueError("alpha must be between 0 and 1, both inclusive.")

    def mean_pinball_loss_q(y_true, y_pred):
        return mean_pinball_loss(y_true, y_pred, alpha=alpha)

    mean_pinball_loss_q.__name__ = f"mean_pinball_loss_q_{alpha}"
    mean_pinball_loss_q.__qualname__ = f"mean_pinball_loss_q_{alpha}"

    return mean_pinball_loss_q

skforecast.metrics.add_y_train_argument

add_y_train_argument(func)

Add y_train argument to a function if it is not already present.

Parameters:

Name	Type	Description	Default
func	callable	Function to which the argument is added.	required

Returns:

Name	Type	Description
wrapper	callable	Function with y_train argument added.

Source code in skforecast/metrics/metrics.py

def add_y_train_argument(func: Callable) -> Callable:
    """
    Add `y_train` argument to a function if it is not already present.

    Parameters
    ----------
    func : callable
        Function to which the argument is added.

    Returns
    -------
    wrapper : callable
        Function with `y_train` argument added.

    """

    # NOTE: Functions already processed (for example, metrics converted by a
    # search before calling the backtesting) are returned as they are, so that
    # `_needs_y_train` is not overwritten.
    if hasattr(func, "_needs_y_train"):
        return func

    sig = inspect.signature(func)

    if "y_train" in sig.parameters:
        func._needs_y_train = True
        return func

    new_params = list(sig.parameters.values()) + [
        inspect.Parameter("y_train", inspect.Parameter.KEYWORD_ONLY, default=None)
    ]
    new_sig = sig.replace(parameters=new_params)

    @wraps(func)
    def wrapper(*args, y_train=None, **kwargs):
        return func(*args, **kwargs)

    wrapper.__signature__ = new_sig
    wrapper._needs_y_train = False

    return wrapper