stagnant models from bugs volume 2

Author: Bob Carpenter

bugs_examples/vol2/stagnant/stagnant.fit.R

@@ -0,0 +1,5 @@
+library('rstan')
+source('stagnant.data.R')
+fit <- stan('stagnant2.stan',
+            data=c("N","x","Y"),
+            chains=4, iter=2000);

bugs_examples/vol2/stagnant/stagnant.stan

@@ -0,0 +1,42 @@
+# Change point model with very poor parameterization from BUGS, vol 2
+
+// Bradley P. Carlin; Alan E. Gelfand; Adrian F. M. Smith.
+// Hierarchical Bayesian Analysis of Changepoint Problems
+// Applied Statistics, Vol. 41, No. 2. (1992), pp. 389-405.
+//
+// In these data, X represents the logarithm of the flow rate of
+// water down an inclined channel (grams per centimetre per second) 
+// and Y represents the logarithm of the height of the stagnant
+// surface layer (centimetres) for different surfactants.
+
+
+data {
+  int<lower=0> N;
+  real x[N];
+  real Y[N];
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> alpha;
+  real beta[2];
+  simplex[N] theta;
+}
+model {
+  // local variables
+  real log_probs[N];
+  real mu[N];
+
+  // priors
+  theta ~ dirichlet(rep_vector(0.01,N));
+  alpha  ~ normal(0,5);
+  beta ~ normal(0,5);
+  sigma ~ cauchy(0,5);
+
+  // mixture likelihood
+  for (k in 1:N) {
+    for (n in 1:N)
+      mu[n] <- alpha + if_else(n <= k, beta[1], beta[2]) * (x[n] - x[k]);
+    log_probs[k] <- log(theta[k]) + normal_log(Y, mu, sigma);
+  }
+  lp__ <- lp__ + log_sum_exp(log_probs);
+}

bugs_examples/vol2/stagnant/stagnant.stan.0

@@ -0,0 +1,9 @@
+"N" <- 29
+"x" <-
+c(-1.39, -1.39, -1.08, -1.08, -0.94, -0.8, -0.63, -0.63, -0.25, 
+-0.25, -0.12, -0.12, 0.01, 0.11, 0.11, 0.11, 0.25, 0.25, 0.34, 
+0.34, 0.44, 0.59, 0.7, 0.7, 0.85, 0.85, 0.99, 0.99, 1.19)
+"Y" <-
+c(1.12, 1.12, 0.99, 1.03, 0.92, 0.9, 0.81, 0.83, 0.65, 0.67, 
+0.6, 0.59, 0.51, 0.44, 0.43, 0.43, 0.33, 0.3, 0.25, 0.24, 0.13, 
+-0.01, -0.13, -0.14, -0.3, -0.33, -0.46, -0.43, -0.65)

bugs_examples/vol2/stagnant/stagnant2.data.R

@@ -0,0 +1,25 @@
+data {
+  int<lower=0> N;
+  real x[N];
+  real Y[N];
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> alpha;
+  real beta[2];
+  real<lower=min(x),upper=max(x)> x_change;
+}
+model {
+  real mu[N];
+
+  alpha  ~ normal(0,5);
+  beta ~ normal(0,5);
+  sigma ~ cauchy(0,5);
+
+  for (n in 1:N)
+    mu[n] <- alpha 
+      + if_else(x[n] < x_change, beta[1], beta[2]) * (x[n] - x_change);
+
+  Y ~ normal(mu,sigma);
+}
+