Bearloga has uploaded a new change for review. ( https://gerrit.wikimedia.org/r/344677 )
Change subject: Adds Prophet forecasting ...................................................................... Adds Prophet forecasting Adds forecasting via Facebook's recently released Prophet procedure based on a Bayesian additive model with non-linear trends (https://facebookincubator.github.io/prophet/) Bug: T112170 Change-Id: I3797f8363056d749fd0f4dc8737468531920796d --- M README.md M modules/forecasts/forecast.R M modules/forecasts/models.R M modules/forecasts/search/api_cirrus_arima M modules/forecasts/search/api_cirrus_bsts A modules/forecasts/search/api_cirrus_prophet M modules/forecasts/search/config.yaml M modules/forecasts/search/zrr_overall_arima M modules/forecasts/search/zrr_overall_bsts A modules/forecasts/search/zrr_overall_prophet M modules/forecasts/wdqs/config.yaml M modules/forecasts/wdqs/homepage_traffic_arima M modules/forecasts/wdqs/homepage_traffic_bsts A modules/forecasts/wdqs/homepage_traffic_prophet M modules/forecasts/wdqs/sparql_usage_arima M modules/forecasts/wdqs/sparql_usage_bsts A modules/forecasts/wdqs/sparql_usage_prophet M test.R 18 files changed, 129 insertions(+), 52 deletions(-) git pull ssh://gerrit.wikimedia.org:29418/wikimedia/discovery/golden refs/changes/77/344677/1 diff --git a/README.md b/README.md index f49c6e9..c9a847f 100644 --- a/README.md +++ b/README.md @@ -221,17 +221,21 @@ - [x] Cirrus API usage - [x] [ARIMA-modelled forecasts](modules/forecasts/search/api_cirrus_arima) - [x] [BSTS-modelled forecasts](modules/forecasts/search/api_cirrus_bsts) + - [x] [Prophet-modelled forecasts](modules/forecasts/search/api_cirrus_prophet) - [x] Overall zero results rate - [x] [ARIMA-modelled forecasts](modules/forecasts/search/zrr_overall_arima) - [x] [BSTS-modelled forecasts](modules/forecasts/search/zrr_overall_bsts) + - [x] [Prophet-modelled forecasts](modules/forecasts/search/zrr_overall_prophet) - Wikipedia.org Portal (planned) - [x] WDQS ([configuration](modules/forecasts/wdqs/config.yaml)) - [x] Homepage traffic - [x] [ARIMA-modelled forecasts](modules/forecasts/wdqs/homepage_traffic_arima) - [x] [BSTS-modelled forecasts](modules/forecasts/wdqs/homepage_traffic_bsts) + - [x] [Prophet-modelled forecasts](modules/forecasts/wdqs/homepage_traffic_prophet) - [x] SPARQL endpoint usage - [x] [ARIMA-modelled forecasts](modules/forecasts/wdqs/sparql_usage_arima) - [x] [BSTS-modelled forecasts](modules/forecasts/wdqs/sparql_usage_bsts) + - [x] [Prophet-modelled forecasts](modules/forecasts/wdqs/sparql_usage_prophet) - Maps (planned) - External Traffic (planned) @@ -384,7 +388,7 @@ Forecasting modules assume that all the data is current (hence why they are scheduled to run after the metrics modules in **main.sh**) and the forecast is made for the next day. For example, if backfilling a forecast for 2016-12-01, the model is fit using all available data up to and including 2016-11-30. -There are two model wrappers in [modules/forecasts/models.R](modules/forecasts/models.R): +There are three model wrappers in [modules/forecasts/models.R](modules/forecasts/models.R): - `forecast_arima()` which models the time series via [ARIMA](https://en.wikipedia.org/wiki/Autoregressive_integrated_moving_average) and accepts the following inputs: - `x`: a 1-column `xts` object - `arima_params`: a list w/ order & seasonal components @@ -397,13 +401,17 @@ - `burn_in`: number of MCMC iterations to throw away as burn-in, - `transformation`: a transformation to apply to the data ("none", "log", "logit", or "in millions"); the function back-transforms the predictions to the original scale depending on the transformation chosen - `ar_lags`: number of lags ("p") in the AR(p) process, omitted by default so an AR(p) state component is *NOT* added to the state specification +- `forecast_prophet()` which models the time series via [Facebook's Core Data Science team](https://research.fb.com/category/data-science/)'s open source forecasting procedure [Prophet](https://facebookincubator.github.io/prophet/) and accepts the following inputs: + - `x`: a 1-column `xts` object + - `n_iter`: number of MCMC samples (default 500). If greater than 0, will perform a full Bayesian inference using [Stan](http://mc-stan.org/) with 4 chains. If 0, will do perform a fast [maximum a posteriori probability (MAP) estimatation](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation). + - `transformation`: a transformation to apply to the data ("none", "log", "logit", or "in millions"); the function back-transforms the predictions to the original scale depending on the transformation chosen When adding a new forecasting module, add a script-type report to the respective **config.yaml** and use the following template for the script: ```bash #!/bin/bash -Rscript modules/forecasts/forecast.R --date=$2 --metric=[your forecasted metric] --model=[ARIMA [--bootstrap_ci]|BSTS] +Rscript modules/forecasts/forecast.R --date=$2 --metric=[your forecasted metric] --model=[ARIMA [--bootstrap_ci]|BSTS|Prophet] ``` Change the `--metric` and `--model` arguments accordingly. The actual data-reading and metric-forecasting calls are in a switch statement in [modules/forecasts/forecast.R](modules/forecasts/forecast.R). Don't forget to add the forecasted metric to the `--metric` option's help text at the top of **forecast.R** and don't forget to subset the data after reading it in (e.g. `dplyr::filter(data, date < as.Date(opt$date))`) diff --git a/modules/forecasts/forecast.R b/modules/forecasts/forecast.R old mode 100644 new mode 100755 index bcbc850..04e157e --- a/modules/forecasts/forecast.R +++ b/modules/forecasts/forecast.R @@ -13,7 +13,9 @@ make_option("--model", default = NA, action = "store", type = "character", help = "Available: ARIMA, BSTS"), make_option("--iters", default = 10000, action = "store", type = "numeric", - help = "Number of MCMC iterations to keep in BSTS models [default %default]"), + help = "Number of MCMC iterations to keep in BSTS models [default %default] + If using a Prophet model, then setting this to 0 will perform MAP estimation, + otherwise a full Bayesian inference is performed using Stan."), make_option("--burnin", default = 1000, action = "store", type = "numeric", help = "Number of iterations to use as burn-in in BSTS models [default %default]") ) @@ -76,8 +78,10 @@ try( forecast_arima(api_usage[, "cirrus"], arima_params = list(order = c(0, 1, 2), seasonal = list(order = c(2, 1, 1), period = 7))) ) - } else { # BSTS + } else if (opt$model == "BSTS") { forecast_bsts(api_usage[, "cirrus"], transformation = "log", ar_lags = 1, n_iter = opt$iters, burn_in = opt$burnin) + } else if (opt$model == "Prophet") { + forecast_prophet(api_usage[, "cirrus"], transformation = "log", n_iter = opt$iters) } }, @@ -93,8 +97,10 @@ try( forecast_arima(zrr_overall[, "rate"], arima_params = list(order = c(2, 1, 2), seasonal = list(order = c(1, 0, 0), period = 7))) ) - } else { # BSTS + } else if (opt$model == "BSTS") { forecast_bsts(zrr_overall[, "rate"], transformation = "logit", ar_lags = 1, n_iter = opt$iters, burn_in = opt$burnin) + } else if (opt$model == "Prophet") { + forecast_prophet(zrr_overall[, "rate"], transformation = "logit", n_iter = opt$iters) } }, @@ -103,8 +109,10 @@ try( forecast_arima(wdqs_usage[, "homepage"], transformation = "log", arima_params = list(order = c(1, 1, 1), seasonal = list(order = c(1, 0, 0), period = 7))) ) - } else { # BSTS + } else if (opt$model == "BSTS") { forecast_bsts(wdqs_usage[, "homepage"], transformation = "log", ar_lags = 1, n_iter = opt$iters, burn_in = opt$burnin) + } else if (opt$model == "Prophet") { + forecast_prophet(wdqs_usage[, "homepage"], transformation = "log", n_iter = opt$iters) } }, @@ -113,8 +121,10 @@ try( forecast_arima(wdqs_usage[, "sparql"], transformation = "log", arima_params = list(order = c(1, 1, 2), seasonal = list(order = c(1, 0, 0), period = 7))) ) - } else { # BSTS + } else if (opt$model == "BSTS") { forecast_bsts(wdqs_usage[, "sparql"], transformation = "log", ar_lags = 1, n_iter = opt$iters, burn_in = opt$burnin) + } else if (opt$model == "Prophet") { + forecast_prophet(wdqs_usage[, "sparql"], transformation = "log", n_iter = opt$iters) } } diff --git a/modules/forecasts/models.R b/modules/forecasts/models.R old mode 100644 new mode 100755 index bb2ebe8..e4bded1 --- a/modules/forecasts/models.R +++ b/modules/forecasts/models.R @@ -1,11 +1,29 @@ -# .libPaths("/a/discovery/r-library") - suppressPackageStartupMessages(suppressWarnings(suppressMessages({ library(magrittr) # install.packages("tidyverse") library(xts) # install.packages("xts") library(bsts) # install.packages("bsts") library(forecast) # install.packages("forecast") + library(prophet) # install.packages("prophet") }))) + +transforms <- list( + "none" = list( + to = identity, + from = identity + ), + "log" = list( + to = log10, + from = function(x) { return(10^x) } + ), + "logit" = list( + to = function(p) { return(log(p/(1-p))) }, + from = function(x) { return(exp(x)/(exp(x)+1)) } + ), + "in millions" = list( + to = function(x) { return(x/1e6) }, + from = function(x) { return(x*1e6) } + ) +) forecast_arima <- function( x, # a 1-column xts object @@ -15,30 +33,14 @@ transformation = c("none", "log", "logit", "in millions") ) { if (is.null(arima_params)) { - arima_params <- list(order = c(0L, 0L, 0L), - seasonal = list(order = c(0L, 0L, 0L), period = NA)) + arima_params <- list( + order = c(0L, 0L, 0L), + seasonal = list(order = c(0L, 0L, 0L), period = NA) + ) } if (!(transformation[1] %in% c("none", "log", "logit", "in millions"))) { stop("transformation must be one of: 'none', 'log', 'logit', 'in millions'") } - transforms <- list( - "none" = list( - to = identity, - from = identity - ), - "log" = list( - to = log10, - from = function(x) { return(10^x) } - ), - "logit" = list( - to = function(p) { return(log(p/(1-p))) }, - from = function(x) { return(exp(x)/(exp(x)+1)) } - ), - "in millions" = list( - to = function(x) { return(x/1e6) }, - from = function(x) { return(x*1e6) } - ) - ) transform <- transforms[[transformation[1]]] # Fit: fit <- arima(transform$to(x), order = arima_params$order, seasonal = arima_params$seasonal) @@ -68,24 +70,6 @@ if (!(transformation[1] %in% c("none", "log", "logit", "in millions"))) { stop("transformation must be one of: 'none', 'log', 'logit', 'in millions'") } - transforms <- list( - "none" = list( - to = identity, - from = identity - ), - "log" = list( - to = log10, - from = function(x) { return(10^x) } - ), - "logit" = list( - to = function(p) { return(log(p/(1-p))) }, - from = function(x) { return(exp(x)/(exp(x)+1)) } - ), - "in millions" = list( - to = function(x) { return(x/1e6) }, - from = function(x) { return(x*1e6) } - ) - ) transform <- transforms[[transformation[1]]] # Pre-processing because Some days may be missing, so this ensures that we # have something for every day, even if that something is a NA. @@ -103,10 +87,12 @@ ss <- AddAr(ss, y, lags = ar_lags) } # Fit: - model <- bsts(y, family = "gaussian", - state.specification = ss, - niter = burn_in + n_iter, seed = 0, - ping = 0) + model <- bsts( + y, family = "gaussian", + state.specification = ss, + niter = burn_in + n_iter, + seed = 0, ping = 0 + ) # Forecast: predicted <- predict(model, horizon = 1, burn = burn_in) # Post-processing: @@ -125,3 +111,44 @@ # Return: return(output) } + +forecast_prophet <- function( + x, # a 1-column xts object + n_iter = 500, + transformation = c("none", "log", "logit", "in millions") +) { + if (!(transformation[1] %in% c("none", "log", "logit", "in millions"))) { + stop("transformation must be one of: 'none', 'log', 'logit', 'in millions'") + } + transform <- transforms[[transformation[1]]] + df <- data.frame(ds = zoo::index(x), y = transform$to(as.numeric(x))) + # Fit: + model <- prophet( + df, + weekly.seasonality = TRUE, + yearly.seasonality = TRUE, + mcmc.samples = n_iter + ) + predicted_95 <- predict( + model, + df = data.frame(ds = tail(df$ds, 1) + 1), + interval_width = 0.95 + ) %>% + .[, c("yhat_lower", "yhat_upper")] %>% + transform$from() %>% + set_colnames(c("lower_95", "upper_95")) + predicted_80 <- predict( + model, + df = data.frame(ds = tail(df$ds, 1) + 1), + interval_width = 0.80 + ) %>% + .[, c("yhat", "yhat_lower", "yhat_upper")] %>% + transform$from() %>% + set_colnames(c("point_est", "lower_80", "upper_80")) + predicted <- cbind(predicted_80, predicted_95) + if (!is.data.frame(output)) { + stop("Output is not a data frame for some reason?!?") + } + # Return: + return(output) +} diff --git a/modules/forecasts/search/api_cirrus_arima b/modules/forecasts/search/api_cirrus_arima old mode 100644 new mode 100755 diff --git a/modules/forecasts/search/api_cirrus_bsts b/modules/forecasts/search/api_cirrus_bsts old mode 100644 new mode 100755 diff --git a/modules/forecasts/search/api_cirrus_prophet b/modules/forecasts/search/api_cirrus_prophet new file mode 100755 index 0000000..36ecde7 --- /dev/null +++ b/modules/forecasts/search/api_cirrus_prophet @@ -0,0 +1,3 @@ +#!/bin/bash + +Rscript modules/forecasts/forecast.R --date=$2 --metric=search_api_cirrus --model=Prophet --iters=0 diff --git a/modules/forecasts/search/config.yaml b/modules/forecasts/search/config.yaml index 3f5439a..a46deaa 100644 --- a/modules/forecasts/search/config.yaml +++ b/modules/forecasts/search/config.yaml @@ -9,6 +9,11 @@ granularity: days starts: 2017-02-01 type: script + api_cirrus_prophet: + description: Prophet-modelled forecasts of Cirrus API usage by non-automata users + granularity: days + starts: 2017-02-01 + type: script zrr_overall_arima: description: ARIMA-modelled forecasts of zero results rate, excluding known bots/tools granularity: days @@ -19,3 +24,8 @@ granularity: days starts: 2017-02-01 type: script + zrr_overall_prophet: + description: Prophet-modelled forecasts of zero results rate, excluding known bots/tools + granularity: days + starts: 2017-02-01 + type: script diff --git a/modules/forecasts/search/zrr_overall_arima b/modules/forecasts/search/zrr_overall_arima old mode 100644 new mode 100755 diff --git a/modules/forecasts/search/zrr_overall_bsts b/modules/forecasts/search/zrr_overall_bsts old mode 100644 new mode 100755 diff --git a/modules/forecasts/search/zrr_overall_prophet b/modules/forecasts/search/zrr_overall_prophet new file mode 100755 index 0000000..2fef5ca --- /dev/null +++ b/modules/forecasts/search/zrr_overall_prophet @@ -0,0 +1,3 @@ +#!/bin/bash + +Rscript modules/forecasts/forecast.R --date=$2 --metric=search_zrr_overall --model=Prophet --iters=0 diff --git a/modules/forecasts/wdqs/config.yaml b/modules/forecasts/wdqs/config.yaml index ade052d..f4e4e3a 100644 --- a/modules/forecasts/wdqs/config.yaml +++ b/modules/forecasts/wdqs/config.yaml @@ -9,6 +9,11 @@ granularity: days starts: 2017-02-01 type: script + homepage_traffic_prophet: + description: Prophet-modelled forecasts of WDQS homepage traffic by non-automata users + granularity: days + starts: 2017-02-01 + type: script sparql_usage_arima: description: ARIMA-modelled forecasts of WDQS SPARQL endpoint usage by non-automata granularity: days @@ -19,3 +24,8 @@ granularity: days starts: 2017-02-01 type: script + sparql_usage_prophet: + description: Prophet-modelled forecasts of WDQS SPARQL endpoint usage by non-automata + granularity: days + starts: 2017-02-01 + type: script diff --git a/modules/forecasts/wdqs/homepage_traffic_arima b/modules/forecasts/wdqs/homepage_traffic_arima old mode 100644 new mode 100755 diff --git a/modules/forecasts/wdqs/homepage_traffic_bsts b/modules/forecasts/wdqs/homepage_traffic_bsts old mode 100644 new mode 100755 diff --git a/modules/forecasts/wdqs/homepage_traffic_prophet b/modules/forecasts/wdqs/homepage_traffic_prophet new file mode 100755 index 0000000..00615d3 --- /dev/null +++ b/modules/forecasts/wdqs/homepage_traffic_prophet @@ -0,0 +1,3 @@ +#!/bin/bash + +Rscript modules/forecasts/forecast.R --date=$2 --metric=wdqs_homepage --model=Prophet --iters=0 diff --git a/modules/forecasts/wdqs/sparql_usage_arima b/modules/forecasts/wdqs/sparql_usage_arima old mode 100644 new mode 100755 diff --git a/modules/forecasts/wdqs/sparql_usage_bsts b/modules/forecasts/wdqs/sparql_usage_bsts old mode 100644 new mode 100755 diff --git a/modules/forecasts/wdqs/sparql_usage_prophet b/modules/forecasts/wdqs/sparql_usage_prophet new file mode 100755 index 0000000..070878c --- /dev/null +++ b/modules/forecasts/wdqs/sparql_usage_prophet @@ -0,0 +1,3 @@ +#!/bin/bash + +Rscript modules/forecasts/forecast.R --date=$2 --metric=wdqs_sparql --model=Prophet --iters=0 diff --git a/test.R b/test.R index 5cbd1f6..390ceba 100644 --- a/test.R +++ b/test.R @@ -11,7 +11,7 @@ "knitr", # For forecasting modules: - "bsts", "forecast", + "bsts", "forecast", "prophet", # For querying, etc.: "ISOcodes", "uaparser", "ortiz", "wmf", "polloi" -- To view, visit https://gerrit.wikimedia.org/r/344677 To unsubscribe, visit https://gerrit.wikimedia.org/r/settings Gerrit-MessageType: newchange Gerrit-Change-Id: I3797f8363056d749fd0f4dc8737468531920796d Gerrit-PatchSet: 1 Gerrit-Project: wikimedia/discovery/golden Gerrit-Branch: master Gerrit-Owner: Bearloga <[email protected]> _______________________________________________ MediaWiki-commits mailing list [email protected] https://lists.wikimedia.org/mailman/listinfo/mediawiki-commits
