diff --git a/python/prophet/diagnostics.py b/python/prophet/diagnostics.py index 677922d3d..5aa5c2aa7 100644 --- a/python/prophet/diagnostics.py +++ b/python/prophet/diagnostics.py @@ -11,18 +11,32 @@ import concurrent.futures import multiprocessing import sys -from typing import TYPE_CHECKING, Final, Literal, cast +from typing import TYPE_CHECKING, Literal, cast import numpy as np import pandas as pd +from prophet.metrics import ( # noqa: F401 – re-exported for backward compatibility + PERFORMANCE_METRICS, + coverage, + mae, + mape, + mdape, + mse, + performance_metrics, + register_performance_metric, + rmse, + rolling_mean_by_h, + rolling_median_by_h, + smape, +) + if TYPE_CHECKING: from typing import ( Any, Callable, Iterable, Protocol, - TypeAlias, TypeVar, type_check_only, ) @@ -33,8 +47,7 @@ _ModelT = TypeVar('_ModelT', bound=Prophet) - _PerformanceMetric: TypeAlias = Callable[[pd.DataFrame, int], pd.DataFrame] - _PerformanceMetricT = TypeVar("_PerformanceMetricT", bound=_PerformanceMetric) + _PerformanceMetricT = TypeVar("_PerformanceMetricT", bound=Callable[..., Any]) @type_check_only class _SupportsMap(Protocol): @@ -388,385 +401,4 @@ def prophet_copy(m: _ModelT, cutoff: pd.Timestamp | None = None) -> _ModelT: return m2 -PERFORMANCE_METRICS: Final[dict[str, Callable[..., Any]]] = {} -def register_performance_metric(func: _PerformanceMetricT) -> _PerformanceMetricT: - """Register custom performance metric - - Parameters that your metric should contain - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Registered metric should return following - ------- - Dataframe with columns horizon and metric. - """ - PERFORMANCE_METRICS[func.__name__] = func - return func - - -def performance_metrics( - df: pd.DataFrame, - metrics: list[str] | None = None, - rolling_window: float = 0.1, - monthly: bool = False, -) -> pd.DataFrame | None: - """Compute performance metrics from cross-validation results. - - Computes a suite of performance metrics on the output of cross-validation. - By default, the following metrics are included: - 'mse': mean squared error - 'rmse': root mean squared error - 'mae': mean absolute error - 'mape': mean absolute percent error - 'mdape': median absolute percent error - 'smape': symmetric mean absolute percentage error - 'coverage': coverage of the upper and lower intervals - - A subset of these can be specified by passing a list of names as the - `metrics` argument. - - Metrics are calculated over a rolling window of cross validation - predictions, after sorting by horizon. Averaging is first done within each - value of horizon, and then across horizons as needed to reach the window - size. The size of that window (number of simulated forecast points) is - determined by the rolling_window argument, which specifies a proportion of - simulated forecast points to include in each window. rolling_window=0 will - compute it separately for each horizon. The default of rolling_window=0.1 - will use 10% of the rows in df in each window. rolling_window=1 will - compute the metric across all simulated forecast points. The results are - set to the right edge of the window. - - If rolling_window < 0, then metrics are computed at each datapoint with no - averaging (i.e., 'mse' will actually be squared error with no mean). - - The output is a dataframe containing column 'horizon' along with columns - for each of the metrics computed. - - Parameters - ---------- - df: The dataframe returned by cross_validation. - metrics: A list of performance metrics to compute. If not provided, will - use ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', 'coverage']. - rolling_window: Proportion of data to use in each rolling window for - computing the metrics. Should be in [0, 1] to average. - monthly: monthly=True will compute horizons as numbers of calendar months - from the cutoff date, starting from 0 for the cutoff month. - - Returns - ------- - Dataframe with a column for each metric, and column 'horizon' - """ - valid_metrics = ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', 'coverage'] - if metrics is None: - metrics = valid_metrics - # Copy so that removing skipped metrics below does not modify the caller's list - metrics = list(metrics) - if ('yhat_lower' not in df or 'yhat_upper' not in df) and ('coverage' in metrics): - metrics.remove('coverage') - if len(set(metrics)) != len(metrics): - raise ValueError('Input metrics must be a list of unique values') - if not set(metrics).issubset(set(PERFORMANCE_METRICS)): - raise ValueError( - 'Valid values for metrics are: {}'.format(valid_metrics) - ) - df_m = df.copy() - if monthly: - df_m['horizon'] = df_m['ds'].dt.to_period('M').astype(int) - df_m['cutoff'].dt.to_period('M').astype(int) - else: - df_m['horizon'] = df_m['ds'] - df_m['cutoff'] - df_m.sort_values('horizon', inplace=True) - if 'mape' in metrics and df_m['y'].abs().min() < 1e-8: - logger.info('Skipping MAPE because y close to 0') - metrics.remove('mape') - if len(metrics) == 0: - return None - w = int(rolling_window * df_m.shape[0]) - if w >= 0: - w = max(w, 1) - w = min(w, df_m.shape[0]) - # Compute all metrics - dfs = {} - for metric in metrics: - dfs[metric] = PERFORMANCE_METRICS[metric](df_m, w) - res = dfs[metrics[0]] - for i in range(1, len(metrics)): - res_m = dfs[metrics[i]] - assert np.array_equal(res['horizon'].values, res_m['horizon'].values) - res[metrics[i]] = res_m[metrics[i]] - return res - - -def rolling_mean_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame: - """Compute a rolling mean of x, after first aggregating by h. - - Right-aligned. Computes a single mean for each unique value of h. Each - mean is over at least w samples. - - Parameters - ---------- - x: Array. - h: Array of horizon for each value in x. - w: Integer window size (number of elements). - name: Name for metric in result dataframe - - Returns - ------- - Dataframe with columns horizon and name, the rolling mean of x. - """ - # Aggregate over h - df = pd.DataFrame({'x': x, 'h': h}) - df2 = ( - df.groupby('h').agg(['sum', 'count']).reset_index().sort_values('h') - ) - xs = df2['x']['sum'].values - ns = df2['x']['count'].values - hs = df2.h.values - - trailing_i = len(df2) - 1 - x_sum = 0 - n_sum = 0 - # We don't know output size but it is bounded by len(df2) - res_x = np.empty(len(df2)) - # Start from the right and work backwards - for i in range(len(df2) - 1, -1, -1): - x_sum += xs[i] - n_sum += ns[i] - while n_sum >= w: - # Include points from the previous horizon. All of them if still - # less than w, otherwise weight the mean by the difference - excess_n = n_sum - w - excess_x = excess_n * xs[i] / ns[i] - res_x[trailing_i] = (x_sum - excess_x)/ w - x_sum -= xs[trailing_i] - n_sum -= ns[trailing_i] - trailing_i -= 1 - - res_h = hs[(trailing_i + 1):] - res_x = res_x[(trailing_i + 1):] - - return pd.DataFrame({'horizon': res_h, name: res_x}) - - - -def rolling_median_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame: - """Compute a rolling median of x, after first aggregating by h. - - Right-aligned. Computes a single median for each unique value of h. Each - median is over at least w samples. - - For each h where there are fewer than w samples, we take samples from the previous h, - moving backwards. (In other words, we ~ assume that the x's are shuffled within each h.) - - Parameters - ---------- - x: Array. - h: Array of horizon for each value in x. - w: Integer window size (number of elements). - name: Name for metric in result dataframe - - Returns - ------- - Dataframe with columns horizon and name, the rolling median of x. - """ - # Aggregate over h - df = pd.DataFrame({'x': x, 'h': h}) - grouped = df.groupby('h') - df2 = grouped.size().reset_index().sort_values('h') - hs = df2['h'] - - res_h = [] - res_x = [] - # Start from the right and work backwards - i = len(hs) - 1 - while i >= 0: - h_i = hs[i] - xs = grouped.get_group(h_i).x.tolist() - - # wrap in array so this works if h is pandas Series with custom index or numpy array - next_idx_to_add = np.array(h == h_i).argmax() - 1 - while (len(xs) < w) and (next_idx_to_add >= 0): - # Include points from the previous horizon. All of them if still - # less than w, otherwise just enough to get to w. - xs.append(x[next_idx_to_add]) - next_idx_to_add -= 1 - if len(xs) < w: - # Ran out of points before getting enough. - break - res_h.append(hs[i]) - res_x.append(np.median(xs)) - i -= 1 - res_h.reverse() - res_x.reverse() - return pd.DataFrame({'horizon': res_h, name: res_x}) - - -# The functions below specify performance metrics for cross-validation results. -# Each takes as input the output of cross_validation, and returns the statistic -# as a dataframe, given a window size for rolling aggregation. - - -@register_performance_metric -def mse(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Mean squared error - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and mse. - """ - se = (df['y'] - df['yhat']) ** 2 - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'mse': se}) - return rolling_mean_by_h( - x=cast(np.ndarray, se.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='mse', - ) - - -@register_performance_metric -def rmse(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Root mean squared error - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and rmse. - """ - res = mse(df, w) - res['mse'] = np.sqrt(res['mse']) - res.rename({'mse': 'rmse'}, axis='columns', inplace=True) - return res - - -@register_performance_metric -def mae(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Mean absolute error - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and mae. - """ - ae = (df['y'] - df['yhat']).abs() - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'mae': ae}) - return rolling_mean_by_h( - x=cast(np.ndarray, ae.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='mae', - ) - - -@register_performance_metric -def mape(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Mean absolute percent error - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and mape. - """ - ape = ((df['y'] - df['yhat']) / df['y']).abs() - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'mape': ape}) - return rolling_mean_by_h( - x=cast(np.ndarray, ape.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='mape', - ) - - -@register_performance_metric -def mdape(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Median absolute percent error - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and mdape. - """ - ape = ((df['y'] - df['yhat']) / df['y']).abs() - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'mdape': ape}) - return rolling_median_by_h( - x=cast(np.ndarray, ape.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='mdape', - ) - - -@register_performance_metric -def smape(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Symmetric mean absolute percentage error - based on Chen and Yang (2004) formula - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and smape. - """ - sape = (df['y'] - df['yhat']).abs() / ((np.abs(df['y']) + np.abs(df['yhat'])) / 2) - sape = sape.fillna(0) - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'smape': sape}) - return rolling_mean_by_h( - x=cast(np.ndarray, sape.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='smape', - ) - - -@register_performance_metric -def coverage(df: pd.DataFrame, w: int) -> pd.DataFrame: - """Coverage - - Parameters - ---------- - df: Cross-validation results dataframe. - w: Aggregation window size. - - Returns - ------- - Dataframe with columns horizon and coverage. - """ - is_covered = (df['y'] >= df['yhat_lower']) & (df['y'] <= df['yhat_upper']) - if w < 0: - return pd.DataFrame({'horizon': df['horizon'], 'coverage': is_covered}) - return rolling_mean_by_h( - x=cast(np.ndarray, is_covered.values), - h=cast(np.ndarray, df['horizon'].values), - w=w, - name='coverage', - ) diff --git a/python/prophet/metrics.py b/python/prophet/metrics.py new file mode 100644 index 000000000..388b0e137 --- /dev/null +++ b/python/prophet/metrics.py @@ -0,0 +1,423 @@ +# Copyright (c) Facebook, Inc. and its affiliates. + +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING, Any, Callable, Final, Sequence, cast + +import numpy as np +import pandas as pd + +if TYPE_CHECKING: + from typing import TypeAlias, TypeVar + + _PerformanceMetric: TypeAlias = Callable[[pd.DataFrame, int], pd.DataFrame] + _PerformanceMetricT = TypeVar("_PerformanceMetricT", bound=_PerformanceMetric) + + +logger: logging.Logger = logging.getLogger('prophet') + +PERFORMANCE_METRICS: Final[dict[str, Callable[..., Any]]] = {} + + +def register_performance_metric(func: Callable[..., Any]) -> Callable[..., Any]: + """Register custom performance metric + + Parameters that your metric should contain + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Registered metric should return following + ------- + Dataframe with columns horizon and metric. + """ + PERFORMANCE_METRICS[func.__name__] = func + return func + + +def performance_metrics( + df: pd.DataFrame, + metrics: Sequence[str] | None = None, + rolling_window: float = 0.1, + monthly: bool = False, +) -> pd.DataFrame | None: + """Compute performance metrics from cross-validation results. + + Computes a suite of performance metrics on the output of cross-validation. + By default, the following metrics are included: + 'mse': mean squared error + 'rmse': root mean squared error + 'mae': mean absolute error + 'mape': mean absolute percent error + 'mdape': median absolute percent error + 'smape': symmetric mean absolute percentage error + 'coverage': coverage of the upper and lower intervals + + A subset of these can be specified by passing a list of names as the + `metrics` argument. + + Metrics are calculated over a rolling window of cross validation + predictions, after sorting by horizon. Averaging is first done within each + value of horizon, and then across horizons as needed to reach the window + size. The size of that window (number of simulated forecast points) is + determined by the rolling_window argument, which specifies a proportion of + simulated forecast points to include in each window. rolling_window=0 will + compute it separately for each horizon. The default of rolling_window=0.1 + will use 10% of the rows in df in each window. rolling_window=1 will + compute the metric across all simulated forecast points. The results are + set to the right edge of the window. + + If rolling_window < 0, then metrics are computed at each datapoint with no + averaging (i.e., 'mse' will actually be squared error with no mean). + + The output is a dataframe containing column 'horizon' along with columns + for each of the metrics computed. + + Parameters + ---------- + df: The dataframe returned by cross_validation. + metrics: Performance metrics to compute. A sequence of names; if not + provided, uses ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', + 'coverage']. The caller's sequence is not modified. + rolling_window: Proportion of data to use in each rolling window for + computing the metrics. Should be in [0, 1] to average. + monthly: monthly=True will compute horizons as numbers of calendar months + from the cutoff date, starting from 0 for the cutoff month. + + Returns + ------- + Dataframe with a column for each metric, and column 'horizon' + """ + valid_metrics = ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', 'coverage'] + if metrics is None: + metrics = valid_metrics + # Copy so that removing skipped metrics below does not modify the caller's list + metrics = list(metrics) + if ('yhat_lower' not in df or 'yhat_upper' not in df) and ('coverage' in metrics): + metrics.remove('coverage') + if len(set(metrics)) != len(metrics): + raise ValueError('Input metrics must be a list of unique values') + if not set(metrics).issubset(set(PERFORMANCE_METRICS)): + raise ValueError( + 'Valid values for metrics are: {}'.format(valid_metrics) + ) + df_m = df.copy() + if monthly: + df_m['horizon'] = df_m['ds'].dt.to_period('M').astype(int) - df_m['cutoff'].dt.to_period('M').astype(int) + else: + df_m['horizon'] = df_m['ds'] - df_m['cutoff'] + df_m.sort_values('horizon', inplace=True) + if 'mape' in metrics and df_m['y'].abs().min() < 1e-8: + logger.info('Skipping MAPE because y close to 0') + metrics.remove('mape') + if len(metrics) == 0: + return None + w = int(rolling_window * df_m.shape[0]) + if w >= 0: + w = max(w, 1) + w = min(w, df_m.shape[0]) + # Compute all metrics + dfs = {} + for metric in metrics: + dfs[metric] = PERFORMANCE_METRICS[metric](df_m, w) + res = dfs[metrics[0]] + for i in range(1, len(metrics)): + res_m = dfs[metrics[i]] + assert np.array_equal(res['horizon'].values, res_m['horizon'].values) + res[metrics[i]] = res_m[metrics[i]] + return res + + +def rolling_mean_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame: + """Compute a rolling mean of x, after first aggregating by h. + + Right-aligned. Computes a single mean for each unique value of h. Each + mean is over at least w samples. + + Parameters + ---------- + x: Array. + h: Array of horizon for each value in x. + w: Integer window size (number of elements). + name: Name for metric in result dataframe + + Returns + ------- + Dataframe with columns horizon and name, the rolling mean of x. + """ + # Aggregate over h + df = pd.DataFrame({'x': x, 'h': h}) + df2 = ( + df.groupby('h').agg(['sum', 'count']).reset_index().sort_values('h') + ) + xs = df2['x']['sum'].values + ns = df2['x']['count'].values + hs = df2.h.values + + trailing_i = len(df2) - 1 + x_sum = 0 + n_sum = 0 + # We don't know output size but it is bounded by len(df2) + res_x = np.empty(len(df2)) + + # Start from the right and work backwards + for i in range(len(df2) - 1, -1, -1): + x_sum += xs[i] + n_sum += ns[i] + while n_sum >= w: + # Include points from the previous horizon. All of them if still + # less than w, otherwise weight the mean by the difference + excess_n = n_sum - w + excess_x = excess_n * xs[i] / ns[i] + res_x[trailing_i] = (x_sum - excess_x)/ w + x_sum -= xs[trailing_i] + n_sum -= ns[trailing_i] + trailing_i -= 1 + + res_h = hs[(trailing_i + 1):] + res_x = res_x[(trailing_i + 1):] + + return pd.DataFrame({'horizon': res_h, name: res_x}) + + + +def rolling_median_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame: + """Compute a rolling median of x, after first aggregating by h. + + Right-aligned. Computes a single median for each unique value of h. Each + median is over at least w samples. + + For each h where there are fewer than w samples, we take samples from the previous h, + moving backwards. (In other words, we ~ assume that the x's are shuffled within each h.) + + Parameters + ---------- + x: Array. + h: Array of horizon for each value in x. + w: Integer window size (number of elements). + name: Name for metric in result dataframe + + Returns + ------- + Dataframe with columns horizon and name, the rolling median of x. + """ + # Aggregate over h + df = pd.DataFrame({'x': x, 'h': h}) + grouped = df.groupby('h') + df2 = grouped.size().reset_index().sort_values('h') + hs = df2['h'] + + res_h = [] + res_x = [] + # Start from the right and work backwards + i = len(hs) - 1 + while i >= 0: + h_i = hs[i] + xs = grouped.get_group(h_i).x.tolist() + + # wrap in array so this works if h is pandas Series with custom index or numpy array + next_idx_to_add = np.array(h == h_i).argmax() - 1 + while (len(xs) < w) and (next_idx_to_add >= 0): + # Include points from the previous horizon. All of them if still + # less than w, otherwise just enough to get to w. + xs.append(x[next_idx_to_add]) + next_idx_to_add -= 1 + if len(xs) < w: + # Ran out of points before getting enough. + break + res_h.append(hs[i]) + res_x.append(np.median(xs)) + i -= 1 + res_h.reverse() + res_x.reverse() + return pd.DataFrame({'horizon': res_h, name: res_x}) + + +# The functions below specify performance metrics for cross-validation results. +# Each takes as input the output of cross_validation, and returns the statistic +# as a dataframe, given a window size for rolling aggregation. + + +@register_performance_metric +def mse(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Mean squared error + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and mse. + """ + se = (df['y'] - df['yhat']) ** 2 + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'mse': se}) + return rolling_mean_by_h( + x=cast(np.ndarray, se.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='mse', + ) + + +@register_performance_metric +def rmse(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Root mean squared error + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and rmse. + """ + res = mse(df, w) + res['mse'] = np.sqrt(res['mse']) + res.rename({'mse': 'rmse'}, axis='columns', inplace=True) + return res + + +@register_performance_metric +def mae(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Mean absolute error + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and mae. + """ + ae = (df['y'] - df['yhat']).abs() + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'mae': ae}) + return rolling_mean_by_h( + x=cast(np.ndarray, ae.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='mae', + ) + + +@register_performance_metric +def mape(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Mean absolute percent error + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and mape. + """ + ape = ((df['y'] - df['yhat']) / df['y']).abs() + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'mape': ape}) + return rolling_mean_by_h( + x=cast(np.ndarray, ape.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='mape', + ) + + +@register_performance_metric +def mdape(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Median absolute percent error + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and mdape. + """ + ape = ((df['y'] - df['yhat']) / df['y']).abs() + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'mdape': ape}) + return rolling_median_by_h( + x=cast(np.ndarray, ape.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='mdape', + ) + + +@register_performance_metric +def smape(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Symmetric mean absolute percentage error + based on Chen and Yang (2004) formula + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and smape. + """ + sape = (df['y'] - df['yhat']).abs() / ((np.abs(df['y']) + np.abs(df['yhat'])) / 2) + sape = sape.fillna(0) + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'smape': sape}) + return rolling_mean_by_h( + x=cast(np.ndarray, sape.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='smape', + ) + + +@register_performance_metric +def coverage(df: pd.DataFrame, w: int) -> pd.DataFrame: + """Coverage + + Parameters + ---------- + df: Cross-validation results dataframe. + w: Aggregation window size. + + Returns + ------- + Dataframe with columns horizon and coverage. + """ + is_covered = (df['y'] >= df['yhat_lower']) & (df['y'] <= df['yhat_upper']) + if w < 0: + return pd.DataFrame({'horizon': df['horizon'], 'coverage': is_covered}) + return rolling_mean_by_h( + x=cast(np.ndarray, is_covered.values), + h=cast(np.ndarray, df['horizon'].values), + w=w, + name='coverage', + ) + + +__all__ = [ + "PERFORMANCE_METRICS", + "register_performance_metric", + "performance_metrics", + "rolling_mean_by_h", + "rolling_median_by_h", + "mse", + "rmse", + "mae", + "mape", + "mdape", + "smape", + "coverage", +] diff --git a/python/prophet/plot.py b/python/prophet/plot.py index e450ecb79..96b3be11a 100644 --- a/python/prophet/plot.py +++ b/python/prophet/plot.py @@ -11,8 +11,7 @@ import numpy as np import pandas as pd -# TODO: separate performance_metrics into a different module. there is an implicit circular import between forecaster.py and diagnostics.py -from prophet.diagnostics import performance_metrics +from prophet.metrics import performance_metrics if TYPE_CHECKING: from typing import Literal, Sequence, TypeVar, type_check_only diff --git a/python/prophet/tests/test_metrics.py b/python/prophet/tests/test_metrics.py new file mode 100644 index 000000000..53563814d --- /dev/null +++ b/python/prophet/tests/test_metrics.py @@ -0,0 +1,221 @@ +# Copyright (c) Facebook, Inc. and its affiliates. + +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + +import numpy as np +import pandas as pd +import pytest + +from prophet import metrics + + +class TestRollingMeanByH: + def test_basic(self): + x = np.arange(10) + h = np.arange(10) + df = metrics.rolling_mean_by_h(x=x, h=h, w=1, name="x") + assert np.array_equal(x, df["x"].values) + assert np.array_equal(h, df["horizon"].values) + + def test_window_4(self): + x = np.arange(10) + h = np.arange(10) + df = metrics.rolling_mean_by_h(x, h, w=4, name="x") + assert np.allclose(x[3:] - 1.5, df["x"].values) + assert np.array_equal(np.arange(3, 10), df["horizon"].values) + + def test_grouped_horizons(self): + x = np.arange(10) + h = np.array([1.0, 2.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 7.0, 7.0]) + x_true = np.array([1.0, 5.0, 22.0 / 3]) + h_true = np.array([3.0, 4.0, 7.0]) + df = metrics.rolling_mean_by_h(x, h, w=3, name="x") + assert np.allclose(x_true, df["x"].values) + assert np.array_equal(h_true, df["horizon"].values) + + def test_large_window(self): + x = np.arange(10) + h = np.array([1.0, 2.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 7.0, 7.0]) + df = metrics.rolling_mean_by_h(x, h, w=10, name="x") + assert np.allclose(np.array([7.0]), df["horizon"].values) + assert np.allclose(np.array([4.5]), df["x"].values) + + +class TestRollingMedianByH: + def test_basic(self): + x = np.arange(10) + h = np.arange(10) + df = metrics.rolling_median_by_h(x=x, h=h, w=1, name="x") + assert np.array_equal(x, df["x"].values) + assert np.array_equal(h, df["horizon"].values) + + def test_window_4(self): + x = np.arange(10) + h = np.arange(10) + df = metrics.rolling_median_by_h(x, h, w=4, name="x") + x_true = x[3:] - 1.5 + assert np.allclose(x_true, df["x"].values) + assert np.array_equal(np.arange(3, 10), df["horizon"].values) + + def test_grouped_horizons(self): + x = np.arange(10) + h = np.array([1.0, 2.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 7.0, 7.0]) + x_true = np.array([1.0, 5.0, 8.0]) + h_true = np.array([3.0, 4.0, 7.0]) + df = metrics.rolling_median_by_h(x, h, w=3, name="x") + assert np.allclose(x_true, df["x"].values) + assert np.array_equal(h_true, df["horizon"].values) + + def test_large_window(self): + x = np.arange(10) + h = np.array([1.0, 2.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 7.0, 7.0]) + df = metrics.rolling_median_by_h(x, h, w=10, name="x") + assert np.allclose(np.array([7.0]), df["horizon"].values) + assert np.allclose(np.array([4.5]), df["x"].values) + + +class TestMetricsRegistered: + def test_all_registered(self): + expected = {'mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', 'coverage'} + assert expected.issubset(set(metrics.PERFORMANCE_METRICS.keys())) + for metric_name in expected: + assert metrics.PERFORMANCE_METRICS[metric_name] is getattr(metrics, metric_name) + + def test_register_custom(self): + @metrics.register_performance_metric + def custom_metric(df, w): + return pd.DataFrame({'horizon': df['horizon'], 'custom': 0.0}) + + assert 'custom_metric' in metrics.PERFORMANCE_METRICS + assert metrics.PERFORMANCE_METRICS['custom_metric'] is custom_metric + del metrics.PERFORMANCE_METRICS['custom_metric'] + + +class TestPerformanceMetrics: + def _make_cv_df(self): + np.random.seed(42) + ds = pd.date_range("2020-01-01", periods=20, freq="D") + y = np.random.randn(20).cumsum() + 10 + return pd.DataFrame({ + 'ds': ds, + 'y': y, + 'yhat': y + np.random.randn(20) * 0.5, + 'yhat_lower': y - 1, + 'yhat_upper': y + 1, + 'cutoff': [ds[0]] * 20, + }) + + def test_basic(self): + df = self._make_cv_df() + res = metrics.performance_metrics(df, rolling_window=0) + assert res is not None + assert 'horizon' in res.columns + + def test_all_metrics(self): + df = self._make_cv_df() + res = metrics.performance_metrics(df, rolling_window=-1) + assert res is not None + assert set(res.columns) == { + 'horizon', 'coverage', 'mae', 'mape', 'mdape', 'mse', 'rmse', 'smape' + } + + def test_no_coverage_without_uncertainty(self): + df = self._make_cv_df() + df = df.drop(columns=['yhat_lower', 'yhat_upper']) + res = metrics.performance_metrics(df) + assert res is not None + assert 'coverage' not in res.columns + + def test_skip_mape_near_zero(self): + df = self._make_cv_df() + df['y'] = 0.0 + df['yhat'] = 0.0 + res = metrics.performance_metrics(df, metrics=['mape']) + assert res is None + + def test_does_not_modify_metrics_when_coverage_dropped(self): + df = self._make_cv_df().drop(columns=['yhat_lower', 'yhat_upper']) + caller_metrics = ['mae', 'coverage'] + res = metrics.performance_metrics(df, metrics=caller_metrics) + assert caller_metrics == ['mae', 'coverage'] + assert res is not None + assert 'coverage' not in res.columns + assert 'mae' in res.columns + + def test_does_not_modify_metrics_when_mape_dropped(self): + df = self._make_cv_df() + df['y'] = 0.0 + caller_metrics = ['mae', 'mape'] + res = metrics.performance_metrics(df, metrics=caller_metrics) + assert caller_metrics == ['mae', 'mape'] + assert res is not None + assert 'mape' not in res.columns + assert 'mae' in res.columns + + def test_does_not_modify_metrics_when_all_skipped(self): + df = self._make_cv_df() + df['y'] = 0.0 + caller_metrics = ['mape'] + res = metrics.performance_metrics(df, metrics=caller_metrics) + assert res is None + assert caller_metrics == ['mape'] + + def test_accepts_immutable_metrics_sequence(self): + df = self._make_cv_df().drop(columns=['yhat_lower', 'yhat_upper']) + res = metrics.performance_metrics(df, metrics=('mae', 'coverage')) + assert res is not None + assert 'mae' in res.columns + assert 'coverage' not in res.columns + + def test_invalid_metric(self): + df = self._make_cv_df() + with pytest.raises(ValueError): + metrics.performance_metrics(df, metrics=['mse', 'invalid']) + + def test_duplicate_metrics(self): + df = self._make_cv_df() + with pytest.raises(ValueError): + metrics.performance_metrics(df, metrics=['mse', 'mse']) + + def test_monthly_horizon(self): + ds = pd.date_range("2020-01-01", periods=12, freq="MS") + y = np.arange(12, dtype=float) + cutoff = pd.Timestamp("2020-01-01") + df = pd.DataFrame({ + 'ds': ds, + 'y': y, + 'yhat': y + 0.5, + 'yhat_lower': y - 1, + 'yhat_upper': y + 1, + 'cutoff': [cutoff] * len(ds), + }) + res = metrics.performance_metrics(df, monthly=True, rolling_window=0) + assert res is not None + assert res['horizon'].dtype in (np.int64, int) + + +class TestBackwardCompatibility: + def test_import_from_diagnostics(self): + from prophet.diagnostics import ( + PERFORMANCE_METRICS, + coverage, + mae, + mape, + mdape, + mse, + performance_metrics, + register_performance_metric, + rmse, + rolling_mean_by_h, + rolling_median_by_h, + smape, + ) + assert callable(mse) + assert callable(performance_metrics) + assert isinstance(PERFORMANCE_METRICS, dict) + + def test_import_from_metrics(self): + assert callable(metrics.mse) + assert callable(metrics.performance_metrics) + assert isinstance(metrics.PERFORMANCE_METRICS, dict)