prophet_logistic_growth.stan 1.7 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253
  1. data {
  2. int T; // Sample size
  3. int<lower=1> K; // Number of seasonal vectors
  4. vector[T] t; // Day
  5. vector[T] cap; // Capacities
  6. vector[T] y; // Time-series
  7. int S; // Number of changepoints
  8. matrix[T, S] A; // Split indicators
  9. real t_change[S]; // Index of changepoints
  10. matrix[T,K] X; // season vectors
  11. vector[K] sigmas; // scale on seasonality prior
  12. real<lower=0> tau; // scale on changepoints prior
  13. }
  14. parameters {
  15. real k; // Base growth rate
  16. real m; // offset
  17. vector[S] delta; // Rate adjustments
  18. real<lower=0> sigma_obs; // Observation noise (incl. seasonal variation)
  19. vector[K] beta; // seasonal vector
  20. }
  21. transformed parameters {
  22. vector[S] gamma; // adjusted offsets, for piecewise continuity
  23. vector[S + 1] k_s; // actual rate in each segment
  24. real m_pr;
  25. // Compute the rate in each segment
  26. k_s[1] = k;
  27. for (i in 1:S) {
  28. k_s[i + 1] = k_s[i] + delta[i];
  29. }
  30. // Piecewise offsets
  31. m_pr = m; // The offset in the previous segment
  32. for (i in 1:S) {
  33. gamma[i] = (t_change[i] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
  34. m_pr = m_pr + gamma[i]; // update for the next segment
  35. }
  36. }
  37. model {
  38. //priors
  39. k ~ normal(0, 5);
  40. m ~ normal(0, 5);
  41. delta ~ double_exponential(0, tau);
  42. sigma_obs ~ normal(0, 0.1);
  43. beta ~ normal(0, sigmas);
  44. // Likelihood
  45. y ~ normal(cap ./ (1 + exp(-(k + A * delta) .* (t - (m + A * gamma)))) + X * beta, sigma_obs);
  46. }