Change cross validation performance metrics to first aggregate over horizon before computing rolling average (#839, #863)

This commit is contained in:
Ben Letham 2019-05-08 10:47:26 -07:00
parent e9c017ba01
commit 269133c133
8 changed files with 393 additions and 205 deletions

View file

@ -207,14 +207,18 @@ prophet_copy <- function(m, cutoff = NULL) {
#' `metrics` argument.
#'
#' Metrics are calculated over a rolling window of cross validation
#' predictions, after sorting by horizon. 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 simulated
#' forecast point (i.e., 'mse' will actually be squared error with no mean).
#' 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.
#' predictions, after sorting by horizon. Averaging is first done within each
#' value of the 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.
@ -223,7 +227,7 @@ prophet_copy <- function(m, cutoff = NULL) {
#' @param metrics An array of performance metrics to compute. If not provided,
#' will use c('mse', 'rmse', 'mae', 'mape', 'coverage').
#' @param rolling_window Proportion of data to use in each rolling window for
#' computing the metrics. Should be in [0, 1].
#' computing the metrics. Should be in [0, 1] to average.
#'
#' @return A dataframe with a column for each metric, and column 'horizon'.
#'
@ -244,39 +248,89 @@ performance_metrics <- function(df, metrics = NULL, rolling_window = 0.1) {
df_m <- df
df_m$horizon <- df_m$ds - df_m$cutoff
df_m <- df_m[order(df_m$horizon),]
# Window size
if (('mape' %in% metrics) & (min(abs(df_m$y)) < 1e-8)) {
message('Skipping MAPE because y close to 0')
metrics <- metrics[metrics != 'mape']
}
if (length(metrics) == 0) {
return(NULL)
}
w <- as.integer(rolling_window * nrow(df_m))
if (w >= 0) {
w <- max(w, 1)
w <- min(w, nrow(df_m))
cols <- c('horizon')
for (metric in metrics) {
df_m[[metric]] <- get(metric)(df_m, w)
cols <- c(cols, metric)
}
df_m <- df_m[cols]
return(stats::na.omit(df_m))
# Compute all metrics
dfs = list()
for (metric in metrics) {
dfs[[metric]] <- get(metric)(df_m, w)
}
res <- dfs[[metrics[1]]]
for (i in 2:length(metrics)) {
res_m <- dfs[[metrics[i]]]
stopifnot(res$horizon == res_m$horizon)
res[[metrics[i]]] = res_m[[metrics[i]]]
}
return(res)
}
#' Compute a rolling mean of x
#' Compute a rolling mean of x, after first aggregating by h
#'
#' Right-aligned. Padded with NAs on the front so the output is the same
#' size as x.
#' Right-aligned. Computes a single mean for each unique value of h. Each mean
#' is over at least w samples.
#'
#' @param x Array.
#' @param h Array of horizon for each value in x.
#' @param w Integer window size (number of elements).
#' @param name String name for metric in result dataframe.
#'
#' @return Rolling mean of x with window size w.
#' @return Dataframe with columns horizon and name, the rolling mean of x.
#'
#' @importFrom dplyr "%>%"
#' @keywords internal
rolling_mean <- function(x, w) {
s <- cumsum(c(0, x))
prefix <- rep(NA, w - 1)
return(c(prefix, (s[(w + 1):length(s)] - s[1:(length(s) - w)]) / w))
rolling_mean_by_h <- function(x, h, w, name) {
# Aggregate over h
df <- data.frame(x=x, h=h)
df2 <- df %>%
dplyr::group_by(h) %>%
dplyr::summarise(mean = mean(x), n = dplyr::n())
xm <- df2$mean
ns <- df2$n
hs <- df2$h
res <- data.frame(horizon=c())
res[[name]] <- c()
# Start from the right and work backwards
i <- length(hs)
while (i > 0) {
# Construct a mean of at least w samples
n <- ns[i]
xbar <- xm[i]
j <- i - 1
while ((n < w) & (j > 0)) {
# Include points from the previous horizon. All of them if still less
# than w, otherwise just enough to get to w.
n2 <- min(w - n, ns[j])
xbar <- xbar * (n / (n + n2)) + xm[j] * (n2 / (n + n2))
n <- n + n2
j <- j - 1
}
if (n < w) {
# Ran out of horizons before enough points.
break
}
res.i <- data.frame(horizon=hs[i])
res.i[[name]] <- xbar
res <- rbind(res.i, res)
i <- i - 1
}
return(res)
}
# The functions below specify performance metrics for cross-validation results.
# Each takes as input the output of cross_validation, and returns the statistic
# as an array, given a window size for rolling aggregation.
# as a dataframe, given a window size for rolling aggregation.
#' Mean squared error
#'
@ -288,7 +342,10 @@ rolling_mean <- function(x, w) {
#' @keywords internal
mse <- function(df, w) {
se <- (df$y - df$yhat) ** 2
return(rolling_mean(se, w))
if (w < 0) {
return(data.frame(horizon = df$horizon, mse = se))
}
return(rolling_mean_by_h(x = se, h = df$horizon, w = w, name = 'mse'))
}
#' Root mean squared error
@ -300,7 +357,10 @@ mse <- function(df, w) {
#'
#' @keywords internal
rmse <- function(df, w) {
return(sqrt(mse(df, w)))
res <- mse(df, w)
res$mse <- sqrt(res$mse)
names(res)[names(res) == 'mse'] <- 'rmse'
return(res)
}
#' Mean absolute error
@ -313,7 +373,10 @@ rmse <- function(df, w) {
#' @keywords internal
mae <- function(df, w) {
ae <- abs(df$y - df$yhat)
return(rolling_mean(ae, w))
if (w < 0) {
return(data.frame(horizon = df$horizon, mae = ae))
}
return(rolling_mean_by_h(x = ae, h = df$horizon, w = w, name = 'mae'))
}
#' Mean absolute percent error
@ -326,7 +389,10 @@ mae <- function(df, w) {
#' @keywords internal
mape <- function(df, w) {
ape <- abs((df$y - df$yhat) / df$y)
return(rolling_mean(ape, w))
if (w < 0) {
return(data.frame(horizon = df$horizon, mape = ape))
}
return(rolling_mean_by_h(x = ape, h = df$horizon, w = w, name = 'mape'))
}
#' Coverage
@ -339,5 +405,10 @@ mape <- function(df, w) {
#' @keywords internal
coverage <- function(df, w) {
is_covered <- (df$y >= df$yhat_lower) & (df$y <= df$yhat_upper)
return(rolling_mean(is_covered, w))
if (w < 0) {
return(data.frame(horizon = df$horizon, coverage = is_covered))
}
return(
rolling_mean_by_h(x = is_covered, h = df$horizon, w = w, name = 'coverage')
)
}

View file

@ -489,7 +489,7 @@ dyplot.prophet <- function(x, fcst, uncertainty=TRUE,
#'
#' @export
plot_cross_validation_metric <- function(df_cv, metric, rolling_window=0.1) {
df_none <- performance_metrics(df_cv, metrics = metric, rolling_window = 0)
df_none <- performance_metrics(df_cv, metrics = metric, rolling_window = -1)
df_h <- performance_metrics(
df_cv, metrics = metric, rolling_window = rolling_window
)

View file

@ -91,16 +91,20 @@ test_that("performance_metrics", {
df_cv <- cross_validation(
m, horizon = 4, units = "days", period = 10, initial = 90)
# Aggregation level none
df_none <- performance_metrics(df_cv, rolling_window = 0)
df_none <- performance_metrics(df_cv, rolling_window = -1)
expect_true(all(
sort(colnames(df_none))
== sort(c('horizon', 'coverage', 'mae', 'mape', 'mse', 'rmse'))
))
expect_equal(nrow(df_none), 16)
# Aggregation level 0
df_0 <- performance_metrics(df_cv, rolling_window = 0)
expect_equal(nrow(df_0), 4)
expect_equal(length(unique(df_0$h)), 4)
# Aggregation level 0.2
df_horizon <- performance_metrics(df_cv, rolling_window = 0.2)
expect_equal(nrow(df_horizon), 4)
expect_equal(length(unique(df_horizon$horizon)), 4)
expect_equal(nrow(df_horizon), 14)
# Aggregation level all
df_all <- performance_metrics(df_cv, rolling_window = 1)
expect_equal(nrow(df_all), 1)
@ -112,6 +116,38 @@ test_that("performance_metrics", {
expect_true(all(
sort(colnames(df_horizon)) == sort(c('coverage', 'mse', 'horizon'))
))
# Skip MAPE
df_cv$y[1] <- 0.
df_horizon <- performance_metrics(df_cv, metrics = c('coverage', 'mape'))
expect_true(all(
sort(colnames(df_horizon)) == sort(c('coverage', 'horizon'))
))
df_horizon <- performance_metrics(df_cv, metrics = c('mape'))
expect_null(df_horizon)
})
test_that("rolling_mean", {
skip_if_not(Sys.getenv('R_ARCH') != '/i386')
x <- 0:9
h <- 0:9
df <- prophet:::rolling_mean_by_h(x=x, h=h, w=1, name='x')
expect_equal(x, df$x)
expect_equal(h, df$horizon)
df <- prophet:::rolling_mean_by_h(x=x, h=h, w=4, name='x')
expect_equal(x[4:10] - 1.5, df$x)
expect_equal(3:9, df$horizon)
h <- c(1., 2., 3., 4., 4., 4., 4., 4., 7., 7.)
x.true <- c(1., 5., 22/3)
h.true <- c(3., 4., 7.)
df <- prophet:::rolling_mean_by_h(x=x, h=h, w=3, name='x')
expect_equal(x.true, df$x)
expect_equal(h.true, df$horizon)
df <- prophet:::rolling_mean_by_h(x=x, h=h, w=10, name='x')
expect_equal(c(7.), df$horizon)
expect_equal(c(4.5), df$x)
})
test_that("copy", {

View file

@ -425,11 +425,9 @@ test_that("auto_weekly_seasonality", {
train.w <- DATA[1:N.w, ]
m <- prophet(train.w)
expect_false('weekly' %in% names(m$seasonalities))
expect_warning({
# prophet warning: non-zero return code in optimizing
m <- prophet(train.w, weekly.seasonality = TRUE)
expect_true('weekly' %in% names(m$seasonalities))
})
# Should be False due to weekly spacing
train.w <- DATA[seq(1, nrow(DATA), 7), ]
m <- prophet(train.w)

File diff suppressed because one or more lines are too long

View file

@ -205,14 +205,18 @@ def performance_metrics(df, metrics=None, rolling_window=0.1):
`metrics` argument.
Metrics are calculated over a rolling window of cross validation
predictions, after sorting by horizon. 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 simulated
forecast point (i.e., 'mse' will actually be squared error with no mean).
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.
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.
@ -223,7 +227,7 @@ def performance_metrics(df, metrics=None, rolling_window=0.1):
metrics: A list of performance metrics to compute. If not provided, will
use ['mse', 'rmse', 'mae', 'mape', 'coverage'].
rolling_window: Proportion of data to use in each rolling window for
computing the metrics. Should be in [0, 1].
computing the metrics. Should be in [0, 1] to average
Returns
-------
@ -241,42 +245,84 @@ def performance_metrics(df, metrics=None, rolling_window=0.1):
df_m = df.copy()
df_m['horizon'] = df_m['ds'] - df_m['cutoff']
df_m.sort_values('horizon', inplace=True)
# Window size
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])
cols = ['horizon']
# Compute all metrics
dfs = {}
for metric in metrics:
df_m[metric] = eval(metric)(df_m, w)
cols.append(metric)
df_m = df_m[cols]
return df_m.dropna()
dfs[metric] = eval(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(x, w):
"""Compute a rolling mean of x
def rolling_mean_by_h(x, h, w, name):
"""Compute a rolling mean of x, after first aggregating by h.
Right-aligned. Padded with NaNs on the front so the output is the same
size as x.
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
-------
Rolling mean of x with window size w.
Dataframe with columns horizon and name, the rolling mean of x.
"""
s = np.cumsum(np.insert(x, 0, 0))
prefix = np.empty(w - 1)
prefix.fill(np.nan)
return np.hstack((prefix, (s[w:] - s[:-w]) / float(w))) # right-aligned
# Aggregate over h
df = pd.DataFrame({'x': x, 'h': h})
df2 = (
df.groupby('h').agg(['mean', 'count']).reset_index().sort_values('h')
)
xm = df2['x']['mean'].values
ns = df2['x']['count'].values
hs = df2['h'].values
res_h = []
res_x = []
# Start from the right and work backwards
i = len(hs) - 1
while i >= 0:
# Construct a mean of at least w samples.
n = int(ns[i])
xbar = float(xm[i])
j = i - 1
while ((n < w) and j >= 0):
# Include points from the previous horizon. All of them if still
# less than w, otherwise just enough to get to w.
n2 = min(w - n, ns[j])
xbar = xbar * (n / (n + n2)) + xm[j] * (n2 / (n + n2))
n += n2
j -= 1
if n < w:
# Ran out of horizons before enough points.
break
res_h.append(hs[i])
res_x.append(xbar)
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 an array, given a window size for rolling aggregation.
# as a dataframe, given a window size for rolling aggregation.
def mse(df, w):
@ -289,10 +335,14 @@ def mse(df, w):
Returns
-------
Array of mean squared errors.
Dataframe with columns horizon and mse.
"""
se = (df['y'] - df['yhat']) ** 2
return rolling_mean(se.values, w)
if w < 0:
return pd.DataFrame({'horizon': df['horizon'], 'mse': se})
return rolling_mean_by_h(
x=se.values, h=df['horizon'].values, w=w, name='mse'
)
def rmse(df, w):
@ -305,9 +355,12 @@ def rmse(df, w):
Returns
-------
Array of root mean squared errors.
Dataframe with columns horizon and rmse.
"""
return np.sqrt(mse(df, w))
res = mse(df, w)
res['mse'] = np.sqrt(res['mse'])
res.rename({'mse': 'rmse'}, axis='columns', inplace=True)
return res
def mae(df, w):
@ -320,10 +373,14 @@ def mae(df, w):
Returns
-------
Array of mean absolute errors.
Dataframe with columns horizon and mae.
"""
ae = np.abs(df['y'] - df['yhat'])
return rolling_mean(ae.values, w)
if w < 0:
return pd.DataFrame({'horizon': df['horizon'], 'mae': ae})
return rolling_mean_by_h(
x=ae.values, h=df['horizon'].values, w=w, name='mae'
)
def mape(df, w):
@ -336,10 +393,14 @@ def mape(df, w):
Returns
-------
Array of mean absolute percent errors.
Dataframe with columns horizon and mape.
"""
ape = np.abs((df['y'] - df['yhat']) / df['y'])
return rolling_mean(ape.values, w)
if w < 0:
return pd.DataFrame({'horizon': df['horizon'], 'mape': ape})
return rolling_mean_by_h(
x=ape.values, h=df['horizon'].values, w=w, name='mape'
)
def smape(df, w):
@ -352,10 +413,14 @@ def smape(df, w):
Returns
-------
Array of symmetric mean absolute percent errors.
Dataframe with columns horizon and smape.
"""
sape = np.abs(df['yhat']-df['y']) / ((np.abs(df['y']) + np.abs(df['yhat'])) /2)
return rolling_mean(sape.values, w)
if w < 0:
return pd.DataFrame({'horizon': df['horizon'], 'smape': sape})
return rolling_mean_by_h(
x=sape.values, h=df['horizon'].values, w=w, name='smape'
)
def coverage(df, w):
@ -368,7 +433,11 @@ def coverage(df, w):
Returns
-------
Array of coverages.
Dataframe with columns horizon and coverage.
"""
is_covered = (df['y'] >= df['yhat_lower']) & (df['y'] <= df['yhat_upper'])
return rolling_mean(is_covered.values, w)
if w < 0:
return pd.DataFrame({'horizon': df['horizon'], 'coverage': is_covered})
return rolling_mean_by_h(
x=is_covered.values, h=df['horizon'].values, w=w, name='coverage'
)

View file

@ -465,7 +465,7 @@ def plot_cross_validation_metric(
else:
fig = ax.get_figure()
# Get the metric at the level of individual predictions, and with the rolling window.
df_none = performance_metrics(df_cv, metrics=[metric], rolling_window=0)
df_none = performance_metrics(df_cv, metrics=[metric], rolling_window=-1)
df_h = performance_metrics(df_cv, metrics=[metric], rolling_window=rolling_window)
# Some work because matplotlib does not handle timedelta

View file

@ -113,21 +113,25 @@ class TestDiagnostics(TestCase):
df_cv = diagnostics.cross_validation(
m, horizon='4 days', period='10 days', initial='90 days')
# Aggregation level none
df_none = diagnostics.performance_metrics(df_cv, rolling_window=0)
df_none = diagnostics.performance_metrics(df_cv, rolling_window=-1)
self.assertEqual(
set(df_none.columns),
{'horizon', 'coverage', 'mae', 'mape', 'mse', 'rmse'},
)
self.assertEqual(df_none.shape[0], 16)
# Aggregation level 0
df_0 = diagnostics.performance_metrics(df_cv, rolling_window=0)
self.assertEqual(len(df_0), 4)
self.assertEqual(len(df_0['horizon'].unique()), 4)
# Aggregation level 0.2
df_horizon = diagnostics.performance_metrics(df_cv, rolling_window=0.2)
self.assertEqual(len(df_horizon), 4)
self.assertEqual(len(df_horizon['horizon'].unique()), 4)
self.assertEqual(df_horizon.shape[0], 14)
# Aggregation level all
df_all = diagnostics.performance_metrics(df_cv, rolling_window=1)
self.assertEqual(df_all.shape[0], 1)
for metric in ['mse', 'mape', 'mae', 'coverage']:
self.assertEqual(df_all[metric].values[0], df_none[metric].mean())
self.assertAlmostEqual(df_all[metric].values[0], df_none[metric].mean())
# Custom list of metrics
df_horizon = diagnostics.performance_metrics(
df_cv, metrics=['coverage', 'mse'],
@ -136,6 +140,41 @@ class TestDiagnostics(TestCase):
set(df_horizon.columns),
{'coverage', 'mse', 'horizon'},
)
# Skip MAPE
df_cv.loc[0, 'y'] = 0.
df_horizon = diagnostics.performance_metrics(
df_cv, metrics=['coverage', 'mape'],
)
self.assertEqual(
set(df_horizon.columns),
{'coverage', 'horizon'},
)
df_horizon = diagnostics.performance_metrics(
df_cv, metrics=['mape'],
)
self.assertIsNone(df_horizon)
def test_rolling_mean(self):
x = np.arange(10)
h = np.arange(10)
df = diagnostics.rolling_mean_by_h(x=x, h=h, w=1, name='x')
self.assertTrue(np.array_equal(x, df['x'].values))
self.assertTrue(np.array_equal(h, df['horizon'].values))
df = diagnostics.rolling_mean_by_h(x, h, w=4, name='x')
self.assertTrue(np.allclose(x[3:] - 1.5, df['x'].values))
self.assertTrue(np.array_equal(np.arange(3, 10), df['horizon'].values))
h = np.array([1., 2., 3., 4., 4., 4., 4., 4., 7., 7.])
x_true = np.array([1.0, 5.0 , 22. / 3])
h_true = np.array([3., 4., 7.])
df = diagnostics.rolling_mean_by_h(x, h, w=3, name='x')
self.assertTrue(np.allclose(x_true, df['x'].values))
self.assertTrue(np.array_equal(h_true, df['horizon'].values))
df = diagnostics.rolling_mean_by_h(x, h, w=10, name='x')
self.assertTrue(np.allclose(np.array([7.]), df['horizon'].values))
self.assertTrue(np.allclose(np.array([4.5]), df['x'].values))
def test_copy(self):
df = DATA_all.copy()