From a3fc16846431aaad86d47b3d9af0c9d9464ccee1 Mon Sep 17 00:00:00 2001
From: Kim Seonghyun <shyeon.kim@scipi.net>
Date: Sun, 26 Feb 2017 04:43:02 +0100
Subject: [PATCH] Add inst/stan

Signed-off-by: Kim Seonghyun <shyeon.kim@scipi.net>
---
 R/inst/stan/prophet_linear_growth.stan   | 40 +++++++++++++++
 R/inst/stan/prophet_logistic_growth.stan | 63 ++++++++++++++++++++++++
 2 files changed, 103 insertions(+)
 create mode 100644 R/inst/stan/prophet_linear_growth.stan
 create mode 100644 R/inst/stan/prophet_logistic_growth.stan

diff --git a/R/inst/stan/prophet_linear_growth.stan b/R/inst/stan/prophet_linear_growth.stan
new file mode 100644
index 0000000..27c6458
--- /dev/null
+++ b/R/inst/stan/prophet_linear_growth.stan
@@ -0,0 +1,40 @@
+data {
+  int T;                                // Sample size
+  int<lower=1> K;                       // Number of seasonal vectors
+  vector[T] t;                            // Day
+  vector[T] y;                            // Time-series
+  int S;                                // Number of split points
+  matrix[T, S] A;                   // Split indicators
+  int s_indx[S];                 // Index of split points
+  matrix[T,K] X;                // season vectors
+  real<lower=0> sigma;              // scale on seasonality prior
+  real<lower=0> tau;                  // scale on changepoints prior
+}
+
+parameters {
+  real k;                            // Base growth rate
+  real m;                            // offset
+  vector[S] delta;                       // Rate adjustments
+  real<lower=0> sigma_obs;               // Observation noise (incl. seasonal variation)
+  vector[K] beta;                    // seasonal vector
+}
+
+transformed parameters {
+  vector[S] gamma;                  // adjusted offsets, for piecewise continuity
+
+  for (i in 1:S) {
+    gamma[i] = -t[s_indx[i]] * delta[i];
+  }
+}
+
+model {
+  //priors
+  k ~ normal(0, 5);
+  m ~ normal(0, 5);
+  delta ~ double_exponential(0, tau);
+  sigma_obs ~ normal(0, 0.5);
+  beta ~ normal(0, sigma);
+
+  // Likelihood
+  y ~ normal((k + A * delta) .* t + (m + A * gamma) + X * beta, sigma_obs);
+}
diff --git a/R/inst/stan/prophet_logistic_growth.stan b/R/inst/stan/prophet_logistic_growth.stan
new file mode 100644
index 0000000..bea299c
--- /dev/null
+++ b/R/inst/stan/prophet_logistic_growth.stan
@@ -0,0 +1,63 @@
+data {
+  int T;                                // Sample size
+  int<lower=1> K;                       // Number of seasonal vectors
+  vector[T] t;                            // Day
+  vector[T] cap;                          // Capacities
+  vector[T] y;                            // Time-series
+  int S;                                // Number of split points
+  matrix[T, S] A;                   // Split indicators
+  int s_indx[S];                 // Index of split points
+  matrix[T,K] X;                    // season vectors
+  real<lower=0> sigma;              // scale on seasonality prior
+  real<lower=0> tau;                  // scale on changepoints prior
+}
+
+
+transformed data {
+  int s_ext[S + 1];  // Segment endpoints
+  for (j in 1:S) {
+    s_ext[j] = s_indx[j];
+  }
+  s_ext[S + 1] = T + 1;  // Used for the m_adj loop below.
+}
+
+
+parameters {
+  real k;                            // Base growth rate
+  real m;                            // offset
+  vector[S] delta;                       // Rate adjustments
+  real<lower=0> sigma_obs;               // Observation noise (incl. seasonal variation)
+  vector[K] beta;                    // seasonal vector
+}
+
+
+transformed parameters {
+  vector[S] gamma;                  // adjusted offsets, for piecewise continuity
+  vector[S + 1] k_s;                 // actual rate in each segment
+  real m_pr;
+
+  // Compute the rate in each segment
+  k_s[1] = k;
+  for (i in 1:S) {
+    k_s[i + 1] = k_s[i] + delta[i];
+  }
+
+  // Piecewise offsets
+  m_pr = m; // The offset in the previous segment
+  for (i in 1:S) {
+    gamma[i] = (t[s_indx[i]] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
+    m_pr = m_pr + gamma[i];  // update for the next segment
+  }
+}
+
+model {
+  //priors
+  k ~ normal(0, 5);
+  m ~ normal(0, 5);
+  delta ~ double_exponential(0, tau);
+  sigma_obs ~ normal(0, 0.1);
+  beta ~ normal(0, sigma);
+
+  // Likelihood
+  y ~ normal(cap ./ (1 + exp(-(k + A * delta) .* (t - (m + A * gamma)))) + X * beta, sigma_obs);
+}