manual updates: separate out time series, add moving average, add standardization for linear regression, with examples

Author: Bob Carpenter

misc/linear-regression/regression.stan

@@ -0,0 +1,17 @@
+data {
+  int<lower=0> N;
+  vector[N] y;
+  vector[N] x;
+}
+parameters {
+  real alpha;
+  real beta;
+  real<lower=0> sigma;
+}
+model {
+  alpha ~ normal(0,10);    
+  beta ~ normal(0,10);
+  sigma ~ cauchy(0,5);
+  for (n in 1:N)
+    y[n] ~ normal(alpha + beta * x[n], sigma);
+}

misc/linear-regression/regression_std.stan

@@ -0,0 +1,30 @@
+data {
+  int<lower=0> N;
+  vector[N] y;
+  vector[N] x;
+}
+transformed data {
+  vector[N] x_std;
+  vector[N] y_std;
+  x_std <- (x - mean(x)) / sd(x);
+  y_std <- (y - mean(y)) / sd(y);
+}
+parameters {
+  real alpha_std;
+  real beta_std;
+  real<lower=0> sigma_std;
+}
+model {
+  alpha_std ~ normal(0,10);    
+  beta_std ~ normal(0,10);
+  sigma_std ~ cauchy(0,5);
+  y_std ~ normal(alpha_std + beta_std * x_std, sigma_std);
+}
+generated quantities {
+  real alpha;
+  real beta;
+  real<lower=0> sigma;
+  alpha <- sd(y) * (alpha_std + beta_std * mean(x) / sd(x)) + mean(y);
+  beta <- beta_std * sd(y) / sd(x);
+  sigma <- sd(y) * sigma_std;
+}

misc/moving-avg/ma2-sim.R

@@ -0,0 +1,24 @@
+mu <- -1.25;
+sigma <- 0.75;
+theta <- c(0.7,0.3);
+T <- 1000;
+T <- 1000;
+y <- rep(0,T);
+epsilon <- rep(0,T);
+
+predictor <- mu;
+y[1] <- rnorm(1,predictor,sigma);
+epsilon[1] <- y[1] - predictor;
+
+predictor <- mu + theta[1] * epsilon[1];
+y[2] <- rnorm(1,predictor,sigma);
+epsilon[2] <- y[2] - predictor;
+
+for (t in 3:T) {
+  predictor <- mu + theta[1] * epsilon[t - 1] + theta[2] * epsilon[t - 2];
+  y[t] <- rnorm(1, predictor, sigma);
+  epsilon[t] <- y[t] - predictor;
+}
+
+library('rstan')
+# fit <- stan('ma2.stan', data=list(T=T, y=y), iter=500, chains=2, init=0);

misc/moving-avg/ma2.stan

@@ -0,0 +1,28 @@
+data {
+  int<lower=3> T;  // number of observations
+  vector[T] y;     // observation at time T
+}
+parameters {
+  real mu;              // mean
+  real<lower=0> sigma;  // error scale
+  vector[2] theta;      // lag coefficients
+}
+transformed parameters {
+  vector[T] epsilon;    // error terms
+  epsilon[1] <- y[1] - mu;
+  epsilon[2] <- y[2] - mu - theta[1] * epsilon[1];
+  for (t in 3:T)
+    epsilon[t] <- ( y[t] - mu
+                    - theta[1] * epsilon[t - 1]
+                    - theta[2] * epsilon[t - 2] );
+}
+model {
+  mu ~ cauchy(0,2.5);
+  theta ~ cauchy(0,2.5);
+  sigma ~ cauchy(0,2.5);
+  for (t in 3:T)
+    y[t] ~ normal(mu 
+                  + theta[1] * epsilon[t - 1]
+                  + theta[2] * epsilon[t - 2],
+                  sigma);
+}

misc/moving-avg/maQ.stan

@@ -0,0 +1,30 @@
+data {
+  int<lower=0> Q;  // num previous noise terms
+  int<lower=3> T;  // num observations
+  vector[T] y;     // observation at time t
+}
+parameters {
+  real mu;              // mean
+  real<lower=0> sigma;  // error scale
+  vector[2] theta;      // error coeff, lag -t
+}
+transformed parameters {
+  vector[T] epsilon;    // error term at time t
+  for (t in 1:T) {
+    epsilon[t] <- y[t] - mu;
+    for (q in 1:min(t-1,Q))
+      epsilon[t] <- epsilon[t] - theta[q] * epsilon[t - q];
+  }
+}
+model {
+  vector[T] eta;
+  mu ~ cauchy(0,2.5);
+  theta ~ cauchy(0,2.5);
+  sigma ~ cauchy(0,2.5);
+  for (t in 1:T) {
+    eta[t] <- mu;
+    for (q in 1:min(t-1,Q))
+      eta[t] <- eta[t] + theta[q] * epsilon[t - q];
+  }
+  y ~ normal(eta,sigma);
+}