Stan code is now in separate files; updated some deprecated Stan functions

Author: gilleslijnzaad

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon.stan

@@ -0,0 +1,25 @@
+// #### Notes to Stan model #######################################################
+// ## Implementation of this model can be difficult to understand for beginners. 
+// ## Therefore I suggest either not trying to understand it and look on WinBUGS
+// ## version or go deep into Stan manual.
+// ################################################################################
+ 
+// ChaSaSoon Censored Data
+data { 
+  int<lower=0> n_fails;
+  int<lower=0> n;
+  int<lower=0> z_observed;
+}
+
+parameters { 
+  real<lower=0.25, upper=1> theta;  // Uniform Prior on Rate Theta
+}
+
+model { 
+  // Observed Data
+  z_observed ~ binomial(n, theta); 
+  
+  // Unobserved Data
+  target += n_fails * log(binomial_cdf(25 | n, theta) 
+                        - binomial_cdf(14 | n, theta));
+}

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon_Stan.R

@@ -2,30 +2,6 @@
 rm(list=ls()) 
 
 library(rstan)
-
-#### Notes to Stan model #######################################################
-## Implementation of this model can be difficult to understand for beginners. 
-## Therefore I suggest either not trying to understand it and look on WinBUGS
-## version or go deep into Stan manual.
-################################################################################
-model <- "
-# ChaSaSoon Censored Data
-data { 
-  int<lower=0> nfails;
-  int<lower=0> n;
-  int<lower=0> z_observed;
-} 
-parameters { 
-  real<lower=.25,upper=1> theta;  // Uniform Prior on Rate Theta
-} 
-model { 
-  // Observed Data
-  z_observed ~ binomial(n, theta); 
-  
-  // Unobserved Data
-  increment_log_prob(nfails * log(binomial_cdf(25, n, theta) 
-                                  - binomial_cdf(14, n, theta)));
-}"
 
 nfails <- 949  
 n <- 50  # Number of questions  
@@ -40,7 +16,7 @@ parameters <- c("theta")
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?rstan".
-samples <- stan(model_code=model,   
+samples <- stan(file="ChaSaSoon.stan",   
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection.stan

@@ -0,0 +1,42 @@
+// #### Notes to Stan model #######################################################
+// ## The model is not very effective.  
+// ## 1) Don't change seed or lower iterations. This model converges slowly so if
+// ##    you change the values, you'll need to increment iterations significantly  
+// ## 2) Code is quite dissimilar to original WinBUGS model - using conditionals 
+// ##    instead of step function. This will happen in further models more often.
+// ##    There is a difference in what functions are efficient in BUGS and Stan.
+// ################################################################################
+
+// Change Detection
+data { 
+  int n;
+  array[n] int t;
+  array[n] real c;
+}
+
+parameters {
+  array[2] real mu;
+  real<lower=0> lambda;
+  real<lower=0, upper=n> tau;
+}
+
+transformed parameters {
+  real<lower=0> sigma;
+
+  sigma = inv_sqrt(lambda);
+}
+model { 
+  // Group Means
+  mu ~ normal(0, sqrt(1000));
+  // Common Precision
+  lambda ~ gamma(0.001, 0.001);
+    
+  // Which Side is Time of Change Point?
+  // Data Come From A Gaussian
+  for (i in 1:n) {
+    if ((t[i] - tau) < 0.0)
+      c[i] ~ normal(mu[1], sigma);
+    else 
+      c[i] ~ normal(mu[2], sigma);
+  }
+}

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection_Stan.R

@@ -3,47 +3,6 @@ rm(list=ls())
 
 library(rstan)
 
-#### Notes to Stan model #######################################################
-## The model is not very effective.  
-## 1) Don't change seed or lower iterations. This model converges slowly so if
-##    you change the values, you'll need to increment iterations significantly  
-## 2) Code is quite dissimilar to original WinBUGS model - using conditionals 
-##    instead of step function. This will happen in further models more often.
-##    There is a difference in what functions are efficient in BUGS and Stan.
-################################################################################
-model <- "
-// Change Detection
-data { 
-  int n;
-  vector[n] t;
-  vector[n] c;
-}
-parameters {
-  vector[2] mu;
-  real<lower=0> lambda;
-  real<lower=0,upper=n> tau;
-} 
-transformed parameters {
-  real<lower=0> sigma;
-
-  sigma <- inv_sqrt(lambda);
-}
-model { 
-  // Group Means
-  mu ~ normal(0, inv_sqrt(.001));
-  // Common Precision
-  lambda ~ gamma(.001, .001);
-    
-  // Which Side is Time of Change Point?
-  // Data Come From A Gaussian
-  for (i in 1:n) {
-    if ((t[i] - tau) < 0.0)
-      c[i] ~ normal(mu[1], sigma);
-    else 
-      c[i] ~ normal(mu[2], sigma);
-  }
-}"
-
 c <- scan("changepointdata.txt")
 n <- length(c)
 t <- 1:n
@@ -57,7 +16,7 @@ parameters <- c("mu", "sigma", "tau")
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?rstan".
-samples <- stan(model_code=model,   
+samples <- stan(file="ChangeDetection.stan",   
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1.stan

@@ -0,0 +1,41 @@
+// #### Notes to Stan model #######################################################
+// ## 1) Multivariate normal distribution in Stan uses covariance matrix instead of 
+// ##    precision matrix.
+// ## 2) Multivariate normal distribution can be (and is) also vectorized.
+// ## 3) Warnings may occur during sampling, ignore them.
+// ################################################################################
+
+// Pearson Correlation
+data { 
+  int<lower=0> n;
+  array[n] vector[2] x;
+}
+
+parameters {
+  vector[2] mu;
+  vector<lower=0>[2] lambda;
+  real<lower=-1, upper=1> r;
+}
+
+transformed parameters {
+  vector<lower=0>[2] sigma;
+  cov_matrix[2] T;
+
+  // Reparameterization
+  sigma[1] = inv_sqrt(lambda[1]);
+  sigma[2] = inv_sqrt(lambda[2]);
+  
+  T[1,1] = square(sigma[1]);
+  T[1,2] = r * sigma[1] * sigma[2];
+  T[2,1] = r * sigma[1] * sigma[2];
+  T[2,2] = square(sigma[2]);
+}
+
+model {
+  // Priors
+  mu ~ normal(0, sqrt(1000));
+  lambda ~ gamma(0.001, 0.001);
+  
+  // Data
+  x ~ multi_normal(mu, T);
+}

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1_Stan.R

@@ -3,44 +3,6 @@ rm(list=ls())
 
 library(rstan)
 
-#### Notes to Stan model #######################################################
-## 1) Multivariate normal distribution in Stan uses covariance matrix instead of 
-##    precision matrix.
-## 2) Multivariate normal distribution can be (and is) also vectorized.
-## 3) Warnings may occur during sampling, ignore them.
-################################################################################
-model <- "
-// Pearson Correlation
-data { 
-  int<lower=0> n;
-  vector[2] x[n];
-}
-parameters {
-  vector[2] mu;
-  vector<lower=0>[2] lambda;
-  real<lower=-1,upper=1> r;
-} 
-transformed parameters {
-  vector<lower=0>[2] sigma;
-  cov_matrix[2] T;
-
-  // Reparameterization
-  sigma[1] <- inv_sqrt(lambda[1]);
-  sigma[2] <- inv_sqrt(lambda[2]);
-  T[1,1] <- square(sigma[1]);
-  T[1,2] <- r * sigma[1] * sigma[2];
-  T[2,1] <- r * sigma[1] * sigma[2];
-  T[2,2] <- square(sigma[2]);
-}
-model {
-  // Priors
-  mu ~ normal(0, inv_sqrt(.001));
-  lambda ~ gamma(.001, .001);
-  
-  // Data
-  x ~ multi_normal(mu, T);
-}"
-
 # Choose a dataset:
 dataset <- 1
 
@@ -94,7 +56,7 @@ parameters <- c("r", "mu", "sigma")
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?rstan".
-samples <- stan(model_code=model,   
+samples <- stan(file="Correlation_1.stan",   
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2.stan

@@ -0,0 +1,44 @@
+// #### Notes to Stan model #######################################################
+// ## 1) All notes from previous model Correlation_1 also apply to this model.
+// ## 2) If you change sigmaerror to c(0.03, 10) as suggested in excercise 5.2.2
+// ##    warning statements will be more frequent and posterior less smooth.
+// ################################################################################
+
+// Pearson Correlation With Uncertainty in Measurement
+data { 
+  int<lower=0> n;
+  array[n] vector[2] x;
+  vector[2] sigmaerror;
+}
+
+parameters {
+  vector[2] mu;
+  vector<lower=0>[2] lambda;
+  real<lower=-1, upper=1> r;
+  array[n] vector[2] y;
+}
+
+transformed parameters {
+  vector<lower=0>[2] sigma;
+  cov_matrix[2] T;
+
+  // Reparameterization
+  sigma[1] = inv_sqrt(lambda[1]);
+  sigma[2] = inv_sqrt(lambda[2]);
+  
+  T[1,1] = square(sigma[1]);
+  T[1,2] = r * sigma[1] * sigma[2];
+  T[2,1] = r * sigma[1] * sigma[2];
+  T[2,2] = square(sigma[2]);
+}
+
+model {
+  // Priors
+  mu ~ normal(0, sqrt(1000));
+  lambda ~ gamma(0.001, 0.001);
+
+  // Data
+  y ~ multi_normal(mu, T);
+  for (i in 1:n)
+    x[i] ~ normal(y[i], sigmaerror);
+}

Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2_Stan.R

@@ -3,48 +3,6 @@ rm(list=ls())
 
 library(rstan)
 
-#### Notes to Stan model #######################################################
-## 1) All notes from previous model Correlation_1 also apply to this model.
-## 2) If you change sigmaerror to c(0.03, 10) as suggested in excercise 5.2.2
-##    warning statements will be more frequent and posterior less smooth.
-################################################################################
-model <- "
-// Pearson Correlation With Uncertainty in Measurement
-data { 
-  int<lower=0> n;
-  vector[2] x[n];
-  vector[2] sigmaerror;
-}
-parameters {
-  vector[2] mu;
-  vector<lower=0>[2] lambda;
-  real<lower=-1,upper=1> r;
-  vector[2] y[n];
-} 
-transformed parameters {
-  vector<lower=0>[2] sigma;
-  cov_matrix[2] T;
-
-  // Reparameterization
-  sigma[1] <- inv_sqrt(lambda[1]);
-  sigma[2] <- inv_sqrt(lambda[2]);
-  
-  T[1,1] <- square(sigma[1]);
-  T[1,2] <- r * sigma[1] * sigma[2];
-  T[2,1] <- r * sigma[1] * sigma[2];
-  T[2,2] <- square(sigma[2]);
-}
-model {
-  // Priors
-  mu ~ normal(0, inv_sqrt(.001));
-  lambda ~ gamma(.001, .001);
-
-  // Data
-  y ~ multi_normal(mu, T);
-  for (i in 1:n)
-    x[i] ~ normal(y[i], sigmaerror);
-}"
-
 x <- matrix(c( .8, 102, 
               1.0,  98, 
                .5, 100,
@@ -72,7 +30,7 @@ parameters <- c("r", "mu", "sigma")
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?rstan".
-samples <- stan(model_code=model,   
+samples <- stan(file="Correlation_2.stan",   
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,