added preliminary 17.1 stan model code

Author: PeterLi2016

ARM/Ch.17/radon_correlation.stan

@@ -0,0 +1,48 @@
+data {
+  int<lower=0> N;
+  int<lower=0> J;
+  vector[N] y;
+  int<lower=0,upper=1> x[N];
+  int county[N];
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> sigma_a;
+  real<lower=0> sigma_b;
+  real mu_a;
+  real mu_b;
+  real<lower=-1,upper=1> rho;
+  matrix[J,2] B;
+}
+model {
+  vector[N] y_hat;
+  vector[J] a;
+  vector[J] b;
+  matrix[J,2] B_hat;
+  matrix[2,2] Sigma_b;
+
+  mu_a ~ normal(0, 100);
+  mu_b ~ normal(0, 100);
+  rho ~ uniform(-1, 1);
+
+  Sigma_b[1,1] <- pow(sigma_a, 2);
+  Sigma_b[2,2] <- pow(sigma_b, 2);
+  Sigma_b[1,2] <- rho * sigma_a * sigma_b;
+  Sigma_b[2,1] <- Sigma_b[1,2];
+
+  for (j in 1:J) {
+    B_hat[j,1] <- mu_a;
+    B_hat[j,2] <- mu_b;
+    B[j,1:2] ~ multi_normal(B_hat[j,],Sigma_b);
+  }
+
+  for (j in 1:J) {
+    a[j] <- B[j,1];
+    b[j] <- B[j,2];
+  }
+
+  for (i in 1:N)
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+
+  y ~ normal(y_hat, sigma);
+}

ARM/Ch.17/radon_multi_varying_coef.stan

@@ -0,0 +1,49 @@
+data {
+  int<lower=0> N;
+  int<lower=0> J;
+  vector[N] y;
+  int<lower=0,upper=1> x[N];
+  int county[N];
+}
+parameters {
+  real<lower=0> sigma;
+  vector<lower=0>[K] sigma_b;
+  vector[K] mu_raw;
+  vector[K] xi;
+  matrix[J,K] B_raw;
+  matrix[K,K] W;
+  matrix[K,K] Tau_b_raw;
+}
+model {
+  vector[N] y_hat;
+  matrix[J,K] B;
+  matrix[K,K] rho_b;
+  matrix[K,K] Sigma_b_raw;
+  vector[K] mu_raw;
+
+  mu_raw ~ normal(0, 100);
+  xi ~ uniform(0, 100);
+
+  mu <- xi .* mu_raw;
+
+  Tau_b_raw ~ wishart(W, K+1);
+  Sigma_b_raw <- inverse(Tau_b_raw);
+
+  for (j in 1:J) {
+    for (k in 1:K)
+      B[j,k] <- xi[k] * B_raw[j,k];
+    B_raw[j,] ~ multi_normal(mu_raw, Sigma_b_raw);    
+  }
+
+  for (k in 1:K) {
+    for (k_prime in 1:K)
+      rho_b[k,k_prime] <- Sigma_b_raw[k,k_prime] 
+        / sqrt(Sigma_b_raw[k,k] * Sigma_b_raw[k_prime,k_prime]);
+    sigma_b[k] <- abs(xi[k]) * sqrt(Sigma_b_raw[k,k]);
+  }
+  
+  for (i in 1:N)
+    y_hat[i] <- dot_product(B[county[i],],X[i,]);
+
+  y ~ normal(y_hat, sigma);
+}

ARM/Ch.17/radon_multi_varying_coef2.stan

@@ -0,0 +1,52 @@
+data {
+  int<lower=0> N;
+  int<lower=0> J;
+  vector[N] y;
+  int<lower=0,upper=1> x[N];
+  int county[N];
+}
+parameters {
+  real<lower=0> sigma;
+  vector<lower=0>[K] sigma_b;
+  vector[K] mu_raw;
+  vector[K] xi;
+  matrix[J,K] B_raw;
+  matrix[K,K] W;
+  matrix[K,K] Tau_b_raw;
+  matrix[N,K] X_0;  
+  vector[K] b_0;      
+}
+model {
+  vector[N] y_hat;
+  matrix[J,K] B;
+  matrix[K,K] rho_b;
+  matrix[K,K] Sigma_b_raw;
+  vector[K] mu_raw;
+
+  mu_raw ~ normal(0, 100);
+  xi ~ uniform(0, 100);
+  b_0 ~ normal(0, 100);
+
+  mu <- xi .* mu_raw;
+
+  Tau_b_raw ~ wishart(W, K+1);
+  Sigma_b_raw <- inverse(Tau_b_raw);
+
+  for (j in 1:J) {
+    for (k in 1:K)
+      B[j,k] <- xi[k] * B_raw[j,k];
+    B_raw[j,] ~ multi_normal(mu_raw, Sigma_b_raw);    
+  }
+
+  for (k in 1:K) {
+    for (k_prime in 1:K)
+      rho_b[k,k_prime] <- Sigma_b_raw[k,k_prime] 
+        / sqrt(Sigma_b_raw[k,k] * Sigma_b_raw[k_prime,k_prime]);
+    sigma_b[k] <- abs(xi[k]) * sqrt(Sigma_b_raw[k,k]);
+  }
+  
+  for (i in 1:N)
+    y_hat[i] <- dot_product(b_0,X_0[i,]) + dot_product(B[county[i],],X[i,]);
+
+  y ~ normal(y_hat, sigma);
+}

ARM/Ch.17/radon_vary_inter_slope.stan

@@ -0,0 +1,30 @@
+data {
+  int<lower=0> N;
+  int<lower=0> J;
+  vector[N] y;
+  int<lower=0,upper=1> x[N];
+  int county[N];
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> sigma_a;
+  real<lower=0> sigma_b;
+  vector[J] a;
+  vector[J] b;
+  real mu_a;
+  real mu_b;
+}
+transformed parameters {
+  vector[N] y_hat;
+
+  for (i in 1:N)
+    y_hat[i] <- a[county[i]] + b[county[i]] * x[i]
+}
+model {
+  mu_a ~ normal(0, 100);
+  mu_b ~ normal(0, 100);
+
+  a ~ normal(mu_a, sigma_a);
+  b ~ normal(mu_b, sigma_b);
+  y ~ normal(y_hat, sigma);
+}

ARM/Ch.17/radon_wishart2.stan

@@ -0,0 +1,62 @@
+data {
+  int<lower=0> N;
+  int<lower=0> J;
+  vector[N] y;
+  int<lower=0,upper=1> x[N];
+  int county[N];
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> sigma_a;
+  real<lower=0> sigma_b;
+  real mu_a_raw;
+  real mu_b_raw;
+  real xi_a;
+  real xi_b;
+  real<lower=-1,upper=1> rho;
+  matrix[J,2] B_raw;
+  matrix[2,2] W;
+  matrix[2,2] Tau_b_raw;
+  vector[K] b_0;
+  matrix[N,K] X_0;
+}
+model {
+  vector[N] y_hat;
+  vector[J] a;
+  vector[J] b;
+  matrix[J,2] B_raw_hat;
+  matrix[J,2] B_raw;
+  matrix[2,2] Sigma_b_raw;
+  real mu_a;
+  real mu_b;
+
+  mu_a_raw ~ normal(0, 100);
+  mu_b_raw ~ normal(0, 100);
+  xi_a ~ uniform(0, 100);
+  xi_b ~ uniform(0, 100);
+  rho ~ uniform(-1, 1);
+  b_0 ~ normal(0, 100);
+
+  mu_a <- xi_a * mu_a_raw;
+  mu_b <- xi_b * mu_b_raw;
+
+  Tau_b_raw ~ wishart(W, 3);
+  Sigma_b_raw <- inverse(Tau_b_raw);
+
+  sigma_a <- xi_a * sqrt(Sigma_b_raw[1,1]);
+  sigma_b <- xi_b * sqrt(Sigma_b_raw[2,2]):
+  rho <- Sigma_b_raw[1,2] / sqrt(Sigma_b_raw[1,1] * Sigma_b_raw[2,2]);
+
+  for (j in 1:J) {
+    B_raw_hat[j,1] <- mu_a_raw;
+    B_raw_hat[j,2] <- mu_b_raw;
+    B_raw[j,1:2] ~ multi_normal(B_raw_hat[j,], Sigma_b_raw);
+    a[j] <- xi_a * B_raw[j,1];
+    b[j] <- xi_b * B_raw[j,2];
+  }
+
+  for (i in 1:N)
+    y_hat[i] <- dot_product(b_0,X_0[i,]) + a[county[i]] + b[county[i]] * x[i]
+
+  y ~ normal(y_hat, sigma);
+}