diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon.stan
new file mode 100644
index 000000000..3621522c4
--- /dev/null
+++ b/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));
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon_Stan.R
index 4381f0861..f8e32bb39 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChaSaSoon_Stan.R
+++ b/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,
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection.stan
new file mode 100644
index 000000000..0c17e1a0e
--- /dev/null
+++ b/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);
+  }
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection_Stan.R
index b5ac1710a..16ff0e604 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/ChangeDetection_Stan.R
+++ b/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,
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1.stan
new file mode 100644
index 000000000..2badd6f39
--- /dev/null
+++ b/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);
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1_Stan.R
index 595c177cf..97187d077 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_1_Stan.R
+++ b/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,
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2.stan
new file mode 100644
index 000000000..1c715a5e2
--- /dev/null
+++ b/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);
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2_Stan.R
index 125b7ba5c..25d947eb8 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Correlation_2_Stan.R
+++ b/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,
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa.stan
new file mode 100644
index 000000000..1e9b15062
--- /dev/null
+++ b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa.stan
@@ -0,0 +1,55 @@
+// #### Notes to Stan model #######################################################
+// ## 1) This is the first time we use simplex data type. Simplex is similar to 
+// ##    vector, but with a property that the sum of all its elements is equal to 1.
+// ## 2) Sampling statements for parameters alpha, beta and gamma could be removed
+// ##    leading to uniform prior on (0, 1) interval which is the same as beta(1, 1)
+// ## 3) Variable n was removed here. Stan doesn't need this information as
+// ##    an argument for multinomial distribution. Always make sure that you know 
+// ##    what arguments are required for a function / sampling statement. In many 
+// ##    cases these are different from BUGS. Very useful for this are last pages 
+// ##    of Stan manual
+// ################################################################################
+
+// Kappa Coefficient of Agreement
+data { 
+  array[4] int<lower=0> y;
+}
+
+parameters {
+  // Underlying Rates
+
+  // Rate Objective Method Decides 'one'
+  real<lower=0,upper=1> alpha;
+
+  // Rate Surrogate Method Decides 'one' When Objective Method Decides 'one'
+  real<lower=0,upper=1> beta;
+
+  // Rate Surrogate Method Decides 'zero' When Objective Method Decides 'zero'
+  real<lower=0,upper=1> gamma;
+} 
+
+transformed parameters {
+  simplex[4] pi;
+  real xi;
+  real psi;
+  real kappa;
+
+  // Probabilities For Each Count
+  pi[1] = alpha * beta;
+  pi[2] = alpha * (1 - beta);
+  pi[3] = (1 - alpha) * (1 - gamma);
+  pi[4] = (1 - alpha) * gamma;
+    
+  // Derived Measures   
+  // Observed Rate of Agreement
+  xi = alpha * beta + (1 - alpha) * gamma ;
+  // Rate of Chance Agreement
+  psi = (pi[1] + pi[2]) * (pi[1] + pi[3]) + (pi[2] + pi[4]) * (pi[3] + pi[4]);  
+  // Chance-Corrected Agreement
+  kappa = (xi - psi) / (1 - psi);
+}
+
+model {
+  // Count Data     
+  y ~ multinomial(pi);
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa_Stan.R
index 6612f80da..19fcf3ebc 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa_Stan.R
+++ b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Kappa_Stan.R
@@ -3,60 +3,6 @@ rm(list=ls())
 
 library(rstan)
 
-#### Notes to Stan model #######################################################
-## 1) This is the first time we use simplex data type. Simplex is similar to 
-##    vector, but with a property that sum of all it's elements is equal to 1.
-## 2) Sampling statements for parameters alpha, beta and gamma could be removed
-##    leading to uniform prior on (0, 1) interval which is the same as beta(1, 1)
-## 3) Variable n was removed here. Stan doesn't need this information as
-##    an argument for multinomial distribution. Always make sure that you know 
-##    what arguments are required for a function / sampling statement. In many 
-##    cases these are different from BUGS. Very useful for this are last pages 
-##    of Stan manual
-################################################################################
-model <- "
-// Kappa Coefficient of Agreement
-data { 
-  int<lower=0> y[4];
-}
-parameters {
-  // Underlying Rates
-  // Rate Objective Method Decides 'one'
-  real<lower=0,upper=1> alpha;
-  // Rate Surrogate Method Decides 'one' When Objective Method Decides 'one'
-  real<lower=0,upper=1> beta;
-  // Rate Surrogate Method Decides 'zero' When Objective Method Decides 'zero'
-  real<lower=0,upper=1> gamma;
-} 
-transformed parameters {
-  simplex[4] pi;
-  real xi;
-  real psi;
-  real kappa;
-
-  // Probabilities For Each Count
-  pi[1] <- alpha * beta;
-  pi[2] <- alpha * (1 - beta);
-  pi[3] <- (1 - alpha) * (1 - gamma);
-  pi[4] <- (1 - alpha) * gamma;
-    
-  // Derived Measures   
-  // Rate Surrogate Method Agrees With the Objective Method
-  xi <- alpha * beta + (1 - alpha) * gamma ;
-  // Rate of Chance Agreement
-  psi <- (pi[1] + pi[2]) * (pi[1] + pi[3]) + (pi[2] + pi[4]) * (pi[3] + pi[4]);  
-  // Chance-Corrected Agreement
-  kappa <- (xi - psi) / (1 - psi);
-}
-model {
-  alpha ~ beta(1, 1);  // could be removed
-  beta ~ beta(1, 1);  // could be removed
-  gamma ~ beta(1, 1);  // could be removed
-
-  // Count Data     
-  y ~ multinomial(pi);
-}"
-
 # CHOOSE a data set:
 # Influenza 
 y <- c(14, 4, 5, 210)
@@ -74,7 +20,7 @@ parameters <- c("kappa", "xi", "psi", "alpha", "beta", "gamma", "pi")
 
 # The following command calls Stan with specific options.
 # For a detailed description type "?rstan".
-samples <- stan(model_code=model,   
+samples <- stan(file="Kappa.stan",   
                 data=data, 
                 init=myinits,  # If not specified, gives random inits
                 pars=parameters,
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes.stan b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes.stan
new file mode 100644
index 000000000..febe295dd
--- /dev/null
+++ b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes.stan
@@ -0,0 +1,42 @@
+// #### Notes to Stan model #######################################################
+// ## Stan has hypergeometric distribution implemented, so in one way the code
+// ## is more intuitive. On the other hand, Stan can't sample discrete parameters,
+// ## therefore we have to increment log probability manually (as we did in 
+// ## Survey example).  
+// ################################################################################
+
+// Planes
+data { 
+  int<lower=1> x;
+  int<lower=1> n;
+  int<lower=0, upper=n> k;
+  int<lower=x> tmax;
+}
+
+transformed data {
+  int<lower=x> tmin;
+    
+  tmin = x + n - k;
+}
+
+transformed parameters {
+  vector[tmax] lp_parts;
+
+  for (t in 1:tmax) 
+    if (t < tmin)
+      lp_parts[t] = log(1.0 / tmax) + negative_infinity();  // Zero probability
+    else
+      lp_parts[t] = log(1.0 / tmax) + hypergeometric_lpmf(k | n, x, t - x);
+}
+
+model {
+  target += log_sum_exp(lp_parts);
+}
+
+generated quantities {
+  int<lower=tmin, upper=tmax> t;
+  simplex[tmax] tp;
+
+  tp = softmax(lp_parts);
+  t = categorical_rng(tp);
+}
diff --git a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes_Stan.R b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes_Stan.R
index 03477c72b..78c8a9f80 100644
--- a/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes_Stan.R
+++ b/Bayesian_Cognitive_Modeling/ParameterEstimation/DataAnalysis/Planes_Stan.R
@@ -3,45 +3,6 @@ rm(list=ls())
 
 library(rstan)
 
-#### Notes to Stan model #######################################################
-## Stan has hypergeometric distribution implemented, so in one way the code
-## is more intuitive. On the other hand, Stan can't sample discrete parameters,
-## therefore we have to increment log probability manually (as we did in 
-## Survey example).  
-################################################################################
-model <- "
-// Planes
-data { 
-  int<lower=1> x;
-  int<lower=1> n;
-  int<lower=0,upper=n> k;
-  int<lower=x> tmax;
-}
-transformed data {
-  int<lower=x> tmin;
-    
-  tmin <- x + n - k;
-}
-transformed parameters {
-  vector[tmax] lp_parts;
-
-  for (t in 1:tmax) 
-    if (t < tmin)
-      lp_parts[t] <- log(1.0 / tmax) + negative_infinity();  // Zero probability
-    else
-      lp_parts[t] <- log(1.0 / tmax) + hypergeometric_log(k, n, x, t - x);
-}
-model {
-  increment_log_prob(log_sum_exp(lp_parts));
-}
-generated quantities {
-  int<lower=tmin,upper=tmax> t;
-  simplex[tmax] tp;
-
-  tp <- softmax(lp_parts);
-  t <- categorical_rng(tp);
-}"
-
 x <- 10  # number of captures
 k <- 4  # number of recaptures from n
 n <- 5  # size of second sample