From 19e95311c27bb89d430b6a236e12fda149c746f4 Mon Sep 17 00:00:00 2001
From: Simon Kim <sandrokim@gmail.com>
Date: Tue, 4 Jul 2017 17:27:55 -0700
Subject: [PATCH] Add a posterior analysis function to R (#238)

* Add a posterior analysis function to R

* Add a predictive_samples function to R
---
 R/NAMESPACE                 |  3 ++-
 R/R/prophet.R               | 46 +++++++++++++++++++++++++++++--------
 R/man/plot.prophet.Rd       |  4 ++--
 R/man/predictive_samples.Rd | 21 +++++++++++++++++
 4 files changed, 62 insertions(+), 12 deletions(-)
 create mode 100644 R/man/predictive_samples.Rd

diff --git a/R/NAMESPACE b/R/NAMESPACE
index 82c16db..5e558a4 100644
--- a/R/NAMESPACE
+++ b/R/NAMESPACE
@@ -1,9 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method(plot,prophet)
 S3method(predict,prophet)
 export(fit.prophet)
 export(make_future_dataframe)
+export(plot.prophet)
+export(predictive_samples)
 export(prophet)
 export(prophet_plot_components)
 import(Rcpp)
diff --git a/R/R/prophet.R b/R/R/prophet.R
index d705169..f81fff8 100644
--- a/R/R/prophet.R
+++ b/R/R/prophet.R
@@ -768,24 +768,24 @@ predict_seasonal_components <- function(m, df) {
   return(component.predictions)
 }
 
-#' Prophet uncertainty intervals.
+#' Prophet posterior predictive samples.
 #'
 #' @param m Prophet object.
 #' @param df Prediction dataframe.
 #'
-#' @return Dataframe with uncertainty intervals.
+#' @return List with posterior predictive samples for each component.
 #'
-predict_uncertainty <- function(m, df) {
+sample_posterior_predictive <- function(m, df) {
   # Sample trend, seasonality, and yhat from the extrapolation model.
   n.iterations <- length(m$params$k)
   samp.per.iter <- max(1, ceiling(m$uncertainty.samples / n.iterations))
   nsamp <- n.iterations * samp.per.iter  # The actual number of samples
-
+  
   seasonal.features <- make_all_seasonality_features(m, df)
   sim.values <- list("trend" = matrix(, nrow = nrow(df), ncol = nsamp),
                      "seasonal" = matrix(, nrow = nrow(df), ncol = nsamp),
                      "yhat" = matrix(, nrow = nrow(df), ncol = nsamp))
-
+  
   for (i in 1:n.iterations) {
     # For each set of parameters from MCMC (or just 1 set for MAP),
     for (j in 1:samp.per.iter) {
@@ -797,11 +797,22 @@ predict_uncertainty <- function(m, df) {
       }
     }
   }
-
+  return(sim.values)
+}  
+  
+#' Prophet uncertainty intervals.
+#'
+#' @param m Prophet object.
+#' @param df Prediction dataframe.
+#'
+#' @return Dataframe with uncertainty intervals.
+#'
+predict_uncertainty <- function(m, df) {
+  sim.values <- sample_posterior_predictive(m, df)
   # Add uncertainty estimates
   lower.p <- (1 - m$interval.width)/2
   upper.p <- (1 + m$interval.width)/2
-
+  
   intervals <- cbind(
     t(apply(t(sim.values$yhat), 2, stats::quantile, c(lower.p, upper.p),
             na.rm = TRUE)),
@@ -810,12 +821,13 @@ predict_uncertainty <- function(m, df) {
     t(apply(t(sim.values$seasonal), 2, stats::quantile, c(lower.p, upper.p),
             na.rm = TRUE))
   ) %>% dplyr::as_data_frame()
-
+  
   colnames(intervals) <- paste(rep(c('yhat', 'trend', 'seasonal'), each=2),
                                c('lower', 'upper'), sep = "_")
   return(intervals)
 }
 
+
 #' Simulate observations from the extrapolated generative model.
 #'
 #' @param m Prophet object.
@@ -1163,4 +1175,20 @@ plot_yearly <- function(m, uncertainty = TRUE, yearly_start = 0) {
   return(gg.yearly)
 }
 
-# fb-block 3
+#' Sample from the posterior predictive distribution.
+#'
+#' @param m Prophet object.
+#' @param df Dataframe with dates for predictions (column ds), and capacity
+#'  (column cap) if logistic growth.
+#'
+#' @return A list with items "trend", "seasonal", and "yhat" containing 
+#'  posterior predictive samples for that component.
+#'
+#' @export
+predictive_samples <- function(m, df) {
+    df <- setup_dataframe(m, df)$df
+    sim.values <- sample_posterior_predictive(m, df)
+    return(sim.values)
+}
+
+# fb-block 3
\ No newline at end of file
diff --git a/R/man/plot.prophet.Rd b/R/man/plot.prophet.Rd
index 25c520e..4595b44 100644
--- a/R/man/plot.prophet.Rd
+++ b/R/man/plot.prophet.Rd
@@ -4,8 +4,8 @@
 \alias{plot.prophet}
 \title{Plot the prophet forecast.}
 \usage{
-\method{plot}{prophet}(x, fcst, uncertainty = TRUE, plot_cap = TRUE,
-  xlabel = "ds", ylabel = "y", ...)
+plot.prophet(x, fcst, uncertainty = TRUE, plot_cap = TRUE, xlabel = "ds",
+  ylabel = "y", ...)
 }
 \arguments{
 \item{x}{Prophet object.}
diff --git a/R/man/predictive_samples.Rd b/R/man/predictive_samples.Rd
new file mode 100644
index 0000000..3ade7e5
--- /dev/null
+++ b/R/man/predictive_samples.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prophet.R
+\name{predictive_samples}
+\alias{predictive_samples}
+\title{Sample from the posterior predictive distribution.}
+\usage{
+predictive_samples(m, df)
+}
+\arguments{
+\item{m}{Prophet object.}
+
+\item{df}{Dataframe with dates for predictions (column ds), and capacity
+(column cap) if logistic growth.}
+}
+\value{
+A list with items "trend", "seasonal", and "yhat" containing 
+ posterior predictive samples for that component.
+}
+\description{
+Sample from the posterior predictive distribution.
+}