prophet.metrics

  1# Copyright (c) Facebook, Inc. and its affiliates.
  2
  3# This source code is licensed under the MIT license found in the
  4# LICENSE file in the root directory of this source tree.
  5
  6from __future__ import annotations
  7
  8import logging
  9from typing import TYPE_CHECKING, Any, Callable, Final, Sequence, cast
 10
 11import numpy as np
 12import pandas as pd
 13
 14if TYPE_CHECKING:
 15    from typing import TypeAlias, TypeVar
 16
 17    _PerformanceMetric: TypeAlias = Callable[[pd.DataFrame, int], pd.DataFrame]
 18    _PerformanceMetricT = TypeVar("_PerformanceMetricT", bound=_PerformanceMetric)
 19
 20
 21logger: logging.Logger = logging.getLogger('prophet')
 22
 23PERFORMANCE_METRICS: Final[dict[str, Callable[..., Any]]] = {}
 24
 25
 26def register_performance_metric(func: Callable[..., Any]) -> Callable[..., Any]:
 27    """Register custom performance metric
 28
 29    Parameters that your metric should contain
 30    ----------
 31    df: Cross-validation results dataframe.
 32    w: Aggregation window size.
 33
 34    Registered metric should return following
 35    -------
 36    Dataframe with columns horizon and metric.
 37    """
 38    PERFORMANCE_METRICS[func.__name__] = func
 39    return func
 40
 41
 42def performance_metrics(
 43    df: pd.DataFrame,
 44    metrics: Sequence[str] | None = None,
 45    rolling_window: float = 0.1,
 46    monthly: bool = False,
 47) -> pd.DataFrame | None:
 48    """Compute performance metrics from cross-validation results.
 49
 50    Computes a suite of performance metrics on the output of cross-validation.
 51    By default, the following metrics are included:
 52    'mse': mean squared error
 53    'rmse': root mean squared error
 54    'mae': mean absolute error
 55    'mape': mean absolute percent error
 56    'mdape': median absolute percent error
 57    'smape': symmetric mean absolute percentage error
 58    'coverage': coverage of the upper and lower intervals
 59
 60    A subset of these can be specified by passing a list of names as the
 61    `metrics` argument.
 62
 63    Metrics are calculated over a rolling window of cross validation
 64    predictions, after sorting by horizon. Averaging is first done within each
 65    value of horizon, and then across horizons as needed to reach the window
 66    size. The size of that window (number of simulated forecast points) is
 67    determined by the rolling_window argument, which specifies a proportion of
 68    simulated forecast points to include in each window. rolling_window=0 will
 69    compute it separately for each horizon. The default of rolling_window=0.1
 70    will use 10% of the rows in df in each window. rolling_window=1 will
 71    compute the metric across all simulated forecast points. The results are
 72    set to the right edge of the window.
 73
 74    If rolling_window < 0, then metrics are computed at each datapoint with no
 75    averaging (i.e., 'mse' will actually be squared error with no mean).
 76
 77    The output is a dataframe containing column 'horizon' along with columns
 78    for each of the metrics computed.
 79
 80    Parameters
 81    ----------
 82    df: The dataframe returned by cross_validation.
 83    metrics: Performance metrics to compute. A sequence of names; if not
 84        provided, uses ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape',
 85        'coverage']. The caller's sequence is not modified.
 86    rolling_window: Proportion of data to use in each rolling window for
 87        computing the metrics. Should be in [0, 1] to average.
 88    monthly: monthly=True will compute horizons as numbers of calendar months
 89        from the cutoff date, starting from 0 for the cutoff month.
 90
 91    Returns
 92    -------
 93    Dataframe with a column for each metric, and column 'horizon'
 94    """
 95    valid_metrics = ['mse', 'rmse', 'mae', 'mape', 'mdape', 'smape', 'coverage']
 96    if metrics is None:
 97        metrics = valid_metrics
 98    # Copy so that removing skipped metrics below does not modify the caller's list
 99    metrics = list(metrics)
100    if ('yhat_lower' not in df or 'yhat_upper' not in df) and ('coverage' in metrics):
101        metrics.remove('coverage')
102    if len(set(metrics)) != len(metrics):
103        raise ValueError('Input metrics must be a list of unique values')
104    if not set(metrics).issubset(set(PERFORMANCE_METRICS)):
105        raise ValueError(
106            'Valid values for metrics are: {}'.format(valid_metrics)
107        )
108    df_m = df.copy()
109    if monthly:
110        df_m['horizon'] = df_m['ds'].dt.to_period('M').astype(int) - df_m['cutoff'].dt.to_period('M').astype(int)
111    else:
112        df_m['horizon'] = df_m['ds'] - df_m['cutoff']
113    df_m.sort_values('horizon', inplace=True)
114    if 'mape' in metrics and df_m['y'].abs().min() < 1e-8:
115        logger.info('Skipping MAPE because y close to 0')
116        metrics.remove('mape')
117    if len(metrics) == 0:
118        return None
119    w = int(rolling_window * df_m.shape[0])
120    if w >= 0:
121        w = max(w, 1)
122        w = min(w, df_m.shape[0])
123    # Compute all metrics
124    dfs = {}
125    for metric in metrics:
126        dfs[metric] = PERFORMANCE_METRICS[metric](df_m, w)
127    res = dfs[metrics[0]]
128    for i in range(1, len(metrics)):
129        res_m = dfs[metrics[i]]
130        assert np.array_equal(res['horizon'].values, res_m['horizon'].values)
131        res[metrics[i]] = res_m[metrics[i]]
132    return res
133
134
135def rolling_mean_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame:
136    """Compute a rolling mean of x, after first aggregating by h.
137
138    Right-aligned. Computes a single mean for each unique value of h. Each
139    mean is over at least w samples.
140
141    Parameters
142    ----------
143    x: Array.
144    h: Array of horizon for each value in x.
145    w: Integer window size (number of elements).
146    name: Name for metric in result dataframe
147
148    Returns
149    -------
150    Dataframe with columns horizon and name, the rolling mean of x.
151    """
152    # Aggregate over h
153    df = pd.DataFrame({'x': x, 'h': h})
154    df2 = (
155        df.groupby('h').agg(['sum', 'count']).reset_index().sort_values('h')
156    )
157    xs = df2['x']['sum'].values
158    ns = df2['x']['count'].values
159    hs = df2.h.values
160
161    trailing_i = len(df2) - 1
162    x_sum = 0
163    n_sum = 0
164    # We don't know output size but it is bounded by len(df2)
165    res_x = np.empty(len(df2))
166
167    # Start from the right and work backwards
168    for i in range(len(df2) - 1, -1, -1):
169        x_sum += xs[i]
170        n_sum += ns[i]
171        while n_sum >= w:
172            # Include points from the previous horizon. All of them if still
173            # less than w, otherwise weight the mean by the difference
174            excess_n = n_sum - w
175            excess_x = excess_n * xs[i] / ns[i]
176            res_x[trailing_i] = (x_sum - excess_x)/ w
177            x_sum -= xs[trailing_i]
178            n_sum -= ns[trailing_i]
179            trailing_i -= 1
180
181    res_h = hs[(trailing_i + 1):]
182    res_x = res_x[(trailing_i + 1):]
183
184    return pd.DataFrame({'horizon': res_h, name: res_x})
185
186
187
188def rolling_median_by_h(x: np.ndarray, h: np.ndarray, w: int, name: str) -> pd.DataFrame:
189    """Compute a rolling median of x, after first aggregating by h.
190
191    Right-aligned. Computes a single median for each unique value of h. Each
192    median is over at least w samples.
193
194    For each h where there are fewer than w samples, we take samples from the previous h,
195    moving backwards. (In other words, we ~ assume that the x's are shuffled within each h.)
196
197    Parameters
198    ----------
199    x: Array.
200    h: Array of horizon for each value in x.
201    w: Integer window size (number of elements).
202    name: Name for metric in result dataframe
203
204    Returns
205    -------
206    Dataframe with columns horizon and name, the rolling median of x.
207    """
208    # Aggregate over h
209    df = pd.DataFrame({'x': x, 'h': h})
210    grouped = df.groupby('h')
211    df2 = grouped.size().reset_index().sort_values('h')
212    hs = df2['h']
213
214    res_h = []
215    res_x = []
216    # Start from the right and work backwards
217    i = len(hs) - 1
218    while i >= 0:
219        h_i = hs[i]
220        xs = grouped.get_group(h_i).x.tolist()
221
222        # wrap in array so this works if h is pandas Series with custom index or numpy array
223        next_idx_to_add = np.array(h == h_i).argmax() - 1
224        while (len(xs) < w) and (next_idx_to_add >= 0):
225            # Include points from the previous horizon. All of them if still
226            # less than w, otherwise just enough to get to w.
227            xs.append(x[next_idx_to_add])
228            next_idx_to_add -= 1
229        if len(xs) < w:
230            # Ran out of points before getting enough.
231            break
232        res_h.append(hs[i])
233        res_x.append(np.median(xs))
234        i -= 1
235    res_h.reverse()
236    res_x.reverse()
237    return pd.DataFrame({'horizon': res_h, name: res_x})
238
239
240# The functions below specify performance metrics for cross-validation results.
241# Each takes as input the output of cross_validation, and returns the statistic
242# as a dataframe, given a window size for rolling aggregation.
243
244
245@register_performance_metric
246def mse(df: pd.DataFrame, w: int) -> pd.DataFrame:
247    """Mean squared error
248
249    Parameters
250    ----------
251    df: Cross-validation results dataframe.
252    w: Aggregation window size.
253
254    Returns
255    -------
256    Dataframe with columns horizon and mse.
257    """
258    se = (df['y'] - df['yhat']) ** 2
259    if w < 0:
260        return pd.DataFrame({'horizon': df['horizon'], 'mse': se})
261    return rolling_mean_by_h(
262        x=cast(np.ndarray, se.values),
263        h=cast(np.ndarray, df['horizon'].values),
264        w=w,
265        name='mse',
266    )
267
268
269@register_performance_metric
270def rmse(df: pd.DataFrame, w: int) -> pd.DataFrame:
271    """Root mean squared error
272
273    Parameters
274    ----------
275    df: Cross-validation results dataframe.
276    w: Aggregation window size.
277
278    Returns
279    -------
280    Dataframe with columns horizon and rmse.
281    """
282    res = mse(df, w)
283    res['mse'] = np.sqrt(res['mse'])
284    res.rename({'mse': 'rmse'}, axis='columns', inplace=True)
285    return res
286
287
288@register_performance_metric
289def mae(df: pd.DataFrame, w: int) -> pd.DataFrame:
290    """Mean absolute error
291
292    Parameters
293    ----------
294    df: Cross-validation results dataframe.
295    w: Aggregation window size.
296
297    Returns
298    -------
299    Dataframe with columns horizon and mae.
300    """
301    ae = (df['y'] - df['yhat']).abs()
302    if w < 0:
303        return pd.DataFrame({'horizon': df['horizon'], 'mae': ae})
304    return rolling_mean_by_h(
305        x=cast(np.ndarray, ae.values),
306        h=cast(np.ndarray, df['horizon'].values),
307        w=w,
308        name='mae',
309    )
310
311
312@register_performance_metric
313def mape(df: pd.DataFrame, w: int) -> pd.DataFrame:
314    """Mean absolute percent error
315
316    Parameters
317    ----------
318    df: Cross-validation results dataframe.
319    w: Aggregation window size.
320
321    Returns
322    -------
323    Dataframe with columns horizon and mape.
324    """
325    ape = ((df['y'] - df['yhat']) / df['y']).abs()
326    if w < 0:
327        return pd.DataFrame({'horizon': df['horizon'], 'mape': ape})
328    return rolling_mean_by_h(
329        x=cast(np.ndarray, ape.values),
330        h=cast(np.ndarray, df['horizon'].values),
331        w=w,
332        name='mape',
333    )
334
335
336@register_performance_metric
337def mdape(df: pd.DataFrame, w: int) -> pd.DataFrame:
338    """Median absolute percent error
339
340    Parameters
341    ----------
342    df: Cross-validation results dataframe.
343    w: Aggregation window size.
344
345    Returns
346    -------
347    Dataframe with columns horizon and mdape.
348    """
349    ape = ((df['y'] - df['yhat']) / df['y']).abs()
350    if w < 0:
351        return pd.DataFrame({'horizon': df['horizon'], 'mdape': ape})
352    return rolling_median_by_h(
353        x=cast(np.ndarray, ape.values),
354        h=cast(np.ndarray, df['horizon'].values),
355        w=w,
356        name='mdape',
357    )
358
359
360@register_performance_metric
361def smape(df: pd.DataFrame, w: int) -> pd.DataFrame:
362    """Symmetric mean absolute percentage error
363    based on Chen and Yang (2004) formula
364
365    Parameters
366    ----------
367    df: Cross-validation results dataframe.
368    w: Aggregation window size.
369
370    Returns
371    -------
372    Dataframe with columns horizon and smape.
373    """
374    sape = (df['y'] - df['yhat']).abs() / ((np.abs(df['y']) + np.abs(df['yhat'])) / 2)
375    sape = sape.fillna(0)
376    if w < 0:
377        return pd.DataFrame({'horizon': df['horizon'], 'smape': sape})
378    return rolling_mean_by_h(
379        x=cast(np.ndarray, sape.values),
380        h=cast(np.ndarray, df['horizon'].values),
381        w=w,
382        name='smape',
383    )
384
385
386@register_performance_metric
387def coverage(df: pd.DataFrame, w: int) -> pd.DataFrame:
388    """Coverage
389
390    Parameters
391    ----------
392    df: Cross-validation results dataframe.
393    w: Aggregation window size.
394
395    Returns
396    -------
397    Dataframe with columns horizon and coverage.
398    """
399    is_covered = (df['y'] >= df['yhat_lower']) & (df['y'] <= df['yhat_upper'])
400    if w < 0:
401        return pd.DataFrame({'horizon': df['horizon'], 'coverage': is_covered})
402    return rolling_mean_by_h(
403        x=cast(np.ndarray, is_covered.values),
404        h=cast(np.ndarray, df['horizon'].values),
405        w=w,
406        name='coverage',
407    )
408
409
410__all__ = [
411    "PERFORMANCE_METRICS",
412    "register_performance_metric",
413    "performance_metrics",
414    "rolling_mean_by_h",
415    "rolling_median_by_h",
416    "mse",
417    "rmse",
418    "mae",
419    "mape",
420    "mdape",
421    "smape",
422    "coverage",
423]