From 5329b967937fa05197caf95cb68caa2892b8ba3e Mon Sep 17 00:00:00 2001
From: jrnold <jeffrey.arnold@gmail.com>
Date: Thu, 24 Dec 2015 15:55:14 -0800
Subject: [PATCH] add Stan samples

---
 samples/Stan/congress.stan | 14 ++++++++++++++
 samples/Stan/dogs.stan     | 31 +++++++++++++++++++++++++++++++
 samples/Stan/schools.stan  | 26 ++++++++++++++++++++++++++
 3 files changed, 71 insertions(+)
 create mode 100644 samples/Stan/congress.stan
 create mode 100644 samples/Stan/dogs.stan
 create mode 100644 samples/Stan/schools.stan

diff --git a/samples/Stan/congress.stan b/samples/Stan/congress.stan
new file mode 100644
index 00000000..6e7408e5
--- /dev/null
+++ b/samples/Stan/congress.stan
@@ -0,0 +1,14 @@
+data {
+  int<lower=0> N;
+  vector[N] incumbency_88;
+  vector[N] vote_86;
+  vector[N] vote_88;
+}
+parameters {
+  vector[3] beta;
+  real<lower=0> sigma;
+}
+model {
+    vote_88 ~ normal(beta[1] + beta[2] * vote_86
+                     + beta[3] * incumbency_88,sigma);
+}
diff --git a/samples/Stan/dogs.stan b/samples/Stan/dogs.stan
new file mode 100644
index 00000000..ebca16d7
--- /dev/null
+++ b/samples/Stan/dogs.stan
@@ -0,0 +1,31 @@
+data {
+  int<lower=0> n_dogs;
+  int<lower=0> n_trials;
+  int<lower=0,upper=1> y[n_dogs,n_trials];
+}
+parameters {
+  vector[3] beta;
+}
+transformed parameters {
+  matrix[n_dogs,n_trials] n_avoid;
+  matrix[n_dogs,n_trials] n_shock;
+  matrix[n_dogs,n_trials] p;
+
+  for (j in 1:n_dogs) {
+    n_avoid[j,1] <- 0;
+    n_shock[j,1] <- 0;
+    for (t in 2:n_trials) {
+      n_avoid[j,t] <- n_avoid[j,t-1] + 1 - y[j,t-1];
+      n_shock[j,t] <- n_shock[j,t-1] + y[j,t-1];
+    }
+    for (t in 1:n_trials)
+      p[j,t] <- beta[1] + beta[2] * n_avoid[j,t] + beta[3] * n_shock[j,t];
+  }
+}
+model {
+  beta ~ normal(0, 100);
+  for (i in 1:n_dogs) {
+    for (j in 1:n_trials)
+      y[i,j] ~ bernoulli_logit(p[i,j]);
+  }
+}
diff --git a/samples/Stan/schools.stan b/samples/Stan/schools.stan
new file mode 100644
index 00000000..171864a1
--- /dev/null
+++ b/samples/Stan/schools.stan
@@ -0,0 +1,26 @@
+data {
+  int<lower=0> N;
+  vector[N] y;
+  vector[N] sigma_y;
+}
+parameters {
+  vector[N] eta;
+  real mu_theta;
+  real<lower=0,upper=100> sigma_eta;
+  real xi;
+}
+transformed parameters {
+  real<lower=0> sigma_theta;
+  vector[N] theta;
+
+  theta <- mu_theta + xi * eta;
+  sigma_theta <- fabs(xi) / sigma_eta;
+}
+model {
+  mu_theta ~ normal(0, 100);
+  sigma_eta ~ inv_gamma(1, 1); //prior distribution can be changed to uniform
+
+  eta ~ normal(0, sigma_eta);
+  xi ~ normal(0, 5);
+  y ~ normal(theta,sigma_y);
+}