# Chapter 10 Introduction to the probabilistic programming language Stan

Stan is a probabilistic programming language for statistical inference written in C++ that can be accessed through several interfaces (e.g., R, Python, Matlab, etc.). Stan uses an advanced dynamic Hamiltonian Monte Carlo algorithm (Betancourt 2016) based on a variant of the No-U-Turn sampler (known as NUTS: Hoffman and Gelman 2014), which is, in general, more efficient than the traditional Gibbs sampler used in other probabilistic languages such as (Win)BUGS (Lunn et al. 2000) and JAGS (Plummer 2016). In this part of the book, we will focus on the package `rstan`

(Guo, Gabry, and Goodrich 2019) that integrates Stan (Carpenter et al. 2017) with R (R Core Team 2019).

In order to understand how to fit a model in Stan and the difficulties we might face, a minimal understanding of the Stan sampling algorithm is needed.
Stan takes advantage of the fact that the *shape* of the posterior distribution is completely determined by the priors and the likelihood we have defined; the phrase completely determined means that we know the unnormalized posterior distribution, which is the upper part of the Bayes rule (abbreviated below as \(up\)). This is because the denominator, or marginal likelihood, “only” constitutes a normalizing constant:

\[\begin{equation} p(\Theta|y) = \cfrac{ p(y|\Theta) \cdot p(\Theta) }{p(y)} \tag{10.1} \end{equation}\]

\[\begin{equation} up(\Theta|y) = p(y|\Theta) \cdot p(\Theta) \tag{10.2} \end{equation}\]

Thus the unnormalized posterior is proportional to the posterior distribution:

\[\begin{equation} up(\Theta|y) \propto p(\Theta|y) \tag{10.3} \end{equation}\]

(The notation \(up(\cdot)\) that we are using here is not standard in statistics; we are using it only for pedagogical convenience.)

The Stan sampler uses Hamiltonian dynamics and treats the vector of parameters, \(\Theta\) (that could range from a vector containing a couple of parameters, e.g., \(<\mu,\sigma>\), to a vector of hundreds of parameters in hierarchical models), as the position of a frictionless particle that glides on the *negative logarithm of the unnormalized posterior*.
That means that high probability places are valleys and low probability places are peaks in this space.^{30} However, Stan doesn’t just let the particle glide until the bottom of this space. If we let that happen, we would find the mode of the posterior distribution, rather than samples. Stan uses a complex algorithm to determine the weight of the particle and the momentum that we apply to it, as well as when to stop the particle trajectory to take a sample. Because we need to know the speed of this particle, Stan needs to be able to calculate the derivative of the log unnormalized posterior with respect to the parameters (recall that speed is the first derivative of position). This means that if the parameter space is differentiable (is relatively smooth, and does not have any break or angle) and if the parameters of the Stan algorithm are well adjusted–as should happen in the warm-up period–these samples are going to represent samples of the true posterior distribution. Bear in mind that the geometry of the posterior has a big influence on whether the algorithm will converge (fast) or not: If the space is very flat, because there isn’t much data and the priors are not informative, then the particle may need to glide for a long time before it gets to a high probability area; if there are several valleys (multimodality) the particle may never leave the vicinity of one of them; and if the space is funnel shaped, the particle may never explore the funnel. One of the reasons for the difficulties in exploring complicated spaces is that the continuous path of the “particle” is discretized and divided into steps, and the *step size* is optimized for the entire posterior space. In spaces that are too complex, such as a funnel, a step size might be too small to explore the wide part of the funnel, but too large to explore the narrow part; we will deal with this problem in section 11.1.2.

Although our following example assumes a vector of two parameters and thus a simple geometry, real world examples can easily have hundreds of parameters defining an unnormalized posterior space with hundreds of dimensions.

One question that might arise here is the following: Given that we already know the shape of the posterior, why do we need samples? After all, the posterior is just the unnormalized posterior multiplied by some number, the normalizing constant.

To make this discussion concrete, let’s say that we have a subject that participates in a memory test, and in each trial we get a noisy score from their true working memory score. We assume that at each trial, the score is elicited with normally distributed noise. If we want to estimate the score and how much the noise makes it vary from trial to trial, we are assuming a normal likelihood and we want to estimate its mean and standard deviation.

We will use simulated data produced by a normal distribution with a true mean of 3 and a true standard deviation of 10:

`## [1] 11.26 -18.74 -11.88 -8.62 -12.89 7.20`

As always, given our prior knowledge, we decide on priors. In this case, we use a log-normal prior for the standard deviation, \(\sigma\), since it can only be positive, but except for that, the prior distributions are quite arbitrary in this example.

\[\begin{equation} \begin{aligned} \mu &\sim \mathit{Normal}(0, 20)\\ \sigma &\sim \mathit{LogNormal}(3, 1) \end{aligned} \tag{10.4} \end{equation}\]

The unnormalized posterior will be the product of the likelihood of each data point times the prior for each parameter:

\[\begin{equation} up(\mu, \sigma |y) = \prod_n^{100} \mathit{Normal}(y_n|\mu, \sigma) \cdot \mathit{Normal}(\mu | 0, 20) \cdot \mathit{LogNormal}(\sigma | 3, 1) \tag{10.5} \end{equation}\]

where \(y = {11.258, -18.741, \ldots}\)

We can also define the unnormalized posterior, \(up(\cdot)\), as a function in R:

```
up <- function(y, mu, sigma) {
dnorm(x = mu, mean = 0, sd = 20) *
dlnorm(x = sigma, mean = 3, sd = 1) *
prod(dnorm(x = y, mean = mu, sd = sigma))
}
```

For example, if we want to know the unnormalized posterior density for the vector of parameters \(<\mu,\sigma> = <0, 5>\), we do the following:

`## [1] 3e-194`

The shape of the unnormalized posterior density is completely defined and it will look like Figure 10.1.

Why is the shape of the unnormalized posterior density not enough? The main reason is that unless we already know which probability distribution we are dealing with (e.g., normal, Bernoulli, etc.) or we can easily integrate it (which can only be done in simpler cases) we cannot do much with the analytical form of the unnormalized posterior: We cannot calculate credible intervals, or know how likely it is that the true score is above or below zero, and even the mean of the posterior is impossible to calculate. This is because the unnormalized posterior distribution represents the general shape of the posterior distribution. With just the shape of an unknown distribution, we can only answer the following question: What is the most (or least) likely value of the vector of parameters? We can answer this question by searching for the highest (or lowest) place in that shape. This leads us to the the maximum a posteriori (MAP) estimate, which is the Bayesian counterpart of the maximum likelihood estimate (MLE). (However, if we can recognize the shape as a distribution, we are in a different situation. In that case, we might know already the formulas for the expectation, variance, etc. This is what we did in chapter 2, but this is an unusual situation in realistic analyses.) As we mentioned before, if we want to get posterior density values, we need the denominator of the Bayes rule (or marginal likelihood), \(p(y)\), which requires integrating the unnormalized posterior. Even this is not too useful if we want to communicate findings, almost every summary statistic requires us to solve more integrals, and except for a handful of cases, these integrals might not have an analytical solution.

As we mentioned before, the only summary that we can get with an unnormalized posterior shape (that we don’t recognize as a familiar distribution) is its mode (or the MAP: the highest point in Figure 10.1, or the lowest point of the negative log unnormalized log posterior). However, even calculating the mode is not always trivial. In simple cases as this one, one can calculate it analytically; but in more complex cases relatively complicated algorithms are needed.

If we want to be able to calculate summary statistics of the posterior distribution (mean, quantiles, etc.), we are going to need samples from this distribution. This is because with enough samples of a probability distribution, we can achieve very good approximations of summary statistics. Stan will take care of returning samples from the posterior distribution, if the log unnormalized posterior distribution is differentiable and can be expressed as follows:^{31}

\[\begin{equation} \log(up(\Theta|y)) = \sum_n \log(p(y_n|\Theta)) + \sum_q \log(p(\Theta_q)) \tag{10.6} \end{equation}\]

where \(n\) indicates each data point and \(q\) each parameter. In our case, this corresponds to the following:

\[\begin{equation} \begin{aligned} \log(up(\mu, \sigma |y)) =& \sum_n^{100} \log(Normal(y_n|\mu, \sigma)) + \log(Normal(\mu | 0, 20)) \\ &+ \log(LogNormal(\sigma | 3, 1)) \end{aligned} \tag{10.7} \end{equation}\]

In the following sections, we’ll see how we can implement this model and many others in Stan.

## 10.1 Stan syntax

A Stan program is usually saved as a `.stan`

file and accessed through R (or other interfaces) and it is organized into a sequence of optional and obligatory blocks, which must be written in order. The Stan language is different from R and it is loosely based on C++; one important aspect to pay attention to is that every statement ends in a semi-colon, `;`

. Blocks (`{}`

) do not end in semi-colons. Some functions in Stan are written in the same way as in R (e.g., `mean`

, `sum`

, `max`

, `min`

). But some are different; when in doubt, Stan documentation can be extremely helpful. In addition, the package `rstan`

provides the function `lookup()`

to look up for translations of functions. For example, in 4.3, we saw that the R function `plogis()`

is needed to convert from log-odds to probability space. If we need it in a Stan program, we can look for it in the following way:

```
## StanFunction Arguments ReturnType
## 227 inv_logit (T x) R
## 260 log_inv_logit (T x) R
## 261 logistic_cdf (reals y, reals mu, reals sigma) real
## 262 logistic_lccdf (reals y , reals mu, reals sigma) real
## 263 logistic_lcdf (reals y , reals mu, reals sigma) real
```

There are three columns in the output of this call. The first one indicates Stan function names, the second one their arguments with their type, and the third one the type they return. Unlike R, Stan is strict with the type of the variables.^{32} In order to decide on which function to use, it is necessary to look at the Stan documentation and find the function that matches our specific needs (for `plogis`

, the corresponding function would be `inv_logit()`

).

Another important difference with R is that every variable needs to be declared at the beginning of a block with its type (real, integer, vector, matrix, etc.). The next two sections exemplify these details through basic Stan programs.

## 10.2 A first simple example with Stan: Normal likelihood

Let’s fit a Stan model to estimate the simple example given at the introduction of this chapter, where we simulate data in R from a normal distribution with a true mean of `3`

and a true standard deviation of `10`

:

As mentioned earlier, Stan code is organized in blocks. The first block indicates what constitutes “data” for the model:

```
data {
int<lower = 1> N; // Total number of trials
vector[N] y; // Score in each trial
}
```

The variable of type `int`

(integer) represents the number of trials. In addition to the type, some constraints can be indicated with `lower`

and `upper`

. In this case, `N`

can’t be smaller than `1`

. These constraints serve as a sanity check; if they are not satisfied, we get an error and the model won’t run. The data are stored in a vector of length `N`

, unlike R, vectors (and matrices and arrays) need to be defined with their dimensions. Comments are indicated with `//`

rather than `#`

.

The next block indicates the parameters of the model:

```
parameters {
real mu;
real<lower = 0> sigma;
}
```

The two parameters are real numbers, and `sigma`

is constrained to be positive.

Finally, we indicate the prior distributions and likelihood functions in the model block:

```
model {
// Priors:
target += normal_lpdf(mu | 0, 20);
target += lognormal_lpdf(sigma | 3, 1);
// Likelihood:
for(i in 1:N)
target += normal_lpdf(y[i] | mu, sigma);
}
```

The variable `target`

is a reserved word in Stan; every statement with `target +=`

adds terms to the unnormalized *log* posterior density. We do this because adding to the unnormalized log posterior amounts to multiplying a term in the numerator of the unnormalized posterior. As explained earlier, Stan uses the shape of the unnormalized posterior to sample from the actual posterior distribution. See Box 10.1 for a more detailed explanation, and see Box 10.2 for alternative notations.

**Box 10.1 ****What does target do?**

We can exemplify how `target`

works with one hypothetical iteration of the sampler.

In every iteration where the sampler explores the posterior space, `mu`

and `sigma`

acquire different values (this is where the Stan algorithm stops the movement of the particle in the Hamiltonian space). Say that in an iteration, `mu = 1.77`

and `sigma = 10.703`

. Then the following happens in the model block:

- At the beginning of the iteration,
`target`

is zero. - The transformations that the sampler
*automatically*does are taken into account. In our case, although`sigma`

is constrained to be positive in our model, inside Stan’s sampler it is transformed to an “unconstrained” space amenable to Hamiltonian Monte Carlo. That is, Stan samples from an auxiliary parameter that ranges from minus infinity to infinity, which is equivalent to`log(sigma)`

. This auxiliary parameter is then exponentiated, when it is incorporated into our model. Because of the mismatch between the constrained parameter space that we defined and the unconstrained space that it is converted to by Stan, an adjustment to the unnormalized posterior is required and added*automatically*. The reasons for this requirement are somewhat complex and will be discussed in section 12. In this particular case, this adjustment (which is the log absolute value of the Jacobian determinant), is equivalent to adding`log(sigma) = 2.371`

to`target`

. - After
`target += normal_lpdf(mu | 0, 20);`

the log of the density of \(\mathit{Normal}(0,20)\) is evaluated at a given sample of mu (specifically 1.77) and this is added to`target`

. In R, this would be`dnorm(x = 1.77, mean = 0, sd = 20, log = TRUE)`

, which is equal to`-3.919`

. Thus,`target`

should be`-3.919 + 2.371 = -1.548`

. - After
`target += lognormal_lpdf(sigma | 3, 1)`

, we add the log of the density of \(\mathit{LogNormal}(3, 1)\) evaluated at`10.703`

to the previous value of the target. In R, this would be`dlnorm(x = 10.703 , mean = 3, sd = 1, log = TRUE)`

, which is equal to`-3.488`

. Thus,`target`

should be updated to`-1.548 + -3.488 = -5.036`

. - After each iteration of the for-loop in the model block, we add to the target the log density of \(\mathit{Normal}( 1.77, 10.703)\) evaluated at each of the values of Y. In R, this would be to add
`sum(dnorm(Y, 1.77, 10.703, log = TRUE))`

(which is equal to`-375.269`

) to the current value of`target`

`-5.036 + -375.269 = -380.305`

.

This means that for the coordinates <mu = 1.77, sigma = 10.703>, the height of the unnormalized posterior would be the value `exp(target) =`

\(\exp( -380.305 ) = 6.852\times 10^{-166}\). Incidentally, the value of `target`

is returned as `lp__`

(log probability) in an object storing a fit model with Stan.

It is possible to expose the value of `target`

, by printing `target()`

inside a Stan model. The value of `target`

after each iteration is named `lp__`

in the Stan object. This can be useful for troubleshooting a problematic model.

**Box 10.2 ****Explicitly incrementing the log probability function ( target) vs. using the sampling (~) notation.**

In this book we specify priors and likelihoods by explicitly incrementing the log-probability function using the following syntax:

`target += pdf_name_lpdf(parameter | ...)`

However, Stan also allows for specifying priors and likelihood with the so-called sampling notation with the following code.

`parameter ~ pdf_name(..)`

Confusingly enough a sampling statement does not perform any actual sampling, and it is meant to be a notational convenience.

The following two lines of code lead to the same behavior in Stan with respect to parameter estimation. There is, nonetheless, an important difference between them.

```
y ~ normal(mu, sigma);
target += normal_lpdf(y | mu, sigma);
```

The important difference is that the sampling notation (the notation with the \(\sim\)) will *drop normalizing constants*. Consider the following formula that corresponds to the log-transformed PDF of a normal distribution:

\[\begin{equation} -log(\sigma) - \frac{log(2 \pi)}{2} - \frac{(y-\mu)^2}{2 \sigma^2} \end{equation}\]

If one uses the sampling notation, Stan will ignore the terms that don’t contain parameters, such as \(- \frac{log(2 \pi)}{2}\), and depending on whether \(\sigma\), \(\mu\), and \(y\) are data or parameters, Stan will either ignore them or not. For example, if \(\sigma\) and \(\mu\) are data and \(y\) is a parameter, it means that \(-log(\sigma)\) is a constant term that can be ignored, but not \(- \frac{(y-\mu)^2}{2 \sigma^2}\) because it contains the parameter \(y\). Dropping constant terms does not affect parameter estimation because it only affects the unnormalized likelihood in the same way in all the parameter space. To make this more concrete, the whole plot in Figure 10.1 will move up or down by some constant amount, and this won’t affect the Hamiltionian dynamics that determine how we sample from the posterior.

The advantage of the sampling notation is that it can be faster (when many terms are ignored), but the disadvantage is that (i) it is not compatible with the calculation of Bayes factor with bridge sampling (see section 15.4 in chapter 15), or the calculation of the log-likelihood for cross validation (see chapter 16), (ii) it misleads us into thinking that Stan is actually sampling the left term in the sampling statement, e.g., drawing \(y\) from a normal distribution in the previous example, when in fact at each step the log-probability (`target`

) is incremented based on the parameter values determined by Hamiltonian dynamics (as explained before), and (iii) it makes it less straightforward to transition to more complex models where the sampling notation cannot be used (as in, for example, mixture models in chapter 19).

If one is not going to use Bayes factor with bridge sampling or cross-validation, the same speed advantage of the sampling notation can also be achieved by incrementing the log-probability with log-unnormalized probability density or mass functions (functions ending with `_lupdf`

or `_lupmf`

). The previous example would be translated into the following:

`target += normal_lupdf(y | mu, sigma);`

We didn’t use curly brackets with the for-loop; this is a common practice if the for-loop has only one line, but brackets can be added and are obligatory if the for-loop spans several lines.

It’s also possible to avoid the for-loop since many functions are vectorized in Stan:

```
model {
// Priors:
target += normal_lpdf(mu | 0, 20);
target += lognormal_lpdf(sigma | 3, 1);
// Likelihood:
target += normal_lpdf(y | mu, sigma);
}
```

The for-loop and vectorized versions give us the same output: The for-loop version evaluated the log-likelihood at each value of `y`

and added it to `target`

. The vectorized version does not create a vector of log-likelihoods; instead, it sums up the log-likelihood evaluated at each element of `y`

and then it adds that to `target`

.

The complete model looks like this:

```
data {
int<lower = 1> N; // Total number of trials
vector[N] y; // Score in each trial
}
parameters {
real mu;
real<lower = 0> sigma;
}
model {
// Priors:
target += normal_lpdf(mu | 0, 20);
target += lognormal_lpdf(sigma | 3, 1);
// Likelihood:
target += normal_lpdf(y | mu, sigma);
}
```

You can save the above code as `normal.stan`

. Alternatively, you can use the version stored in the package `bcogsci`

. (Typing `?stan-normal`

in the R console provides some documentation for the model.) You can access the code of the models of this book by using `system.file("stan_models", "name_of_the_model.stan", package = "bcogsci")`

.

This command just points to a text file that the package `bcogsci`

stores on your computer. You can open it to read the code (with any text editor, or `readLines()`

in R). You’ll need to compile this code and run it with `stan()`

.

Stan requires the data to be in a list object in R. Below, we fit the model with the default number of chains and iterations.

```
Y <- rnorm(n = 100, mean = 3, sd = 10)
lst_score_data <- list(y = Y, N = length(Y))
# Fit the model with the default values of number of
# chains and iterations: chains = 4, iter = 2000
fit_score <- stan(normal, data = lst_score_data)
# alternatively:
# stan("normal.stan", data = lst_score_data)
```

Inspect how well the chains mixed in Figure 10.2. The chains for each parameter should look like a “fat hairy caterpillar” (Lunn et al. 2012); see section 3.2.1.2 for a brief discussion about convergence.

We can see a summary of the posterior by either printing out the model fit, or by plotting it. The summary displayed by the function `print`

includes means, standard deviations (`sd`

), quantiles, Monte Carlo standard errors for the mean of the posterior (`se_mean`

), split Rhats, and effective sample sizes (`n_eff`

). The summaries are computed after removing the warmup and merging together all chains. The `se_mean`

is unrelated to the `se`

of an estimate in the parallel frequentist model. Similarly to a large effective sample size, small Monte Carlo standard errors indicate an “efficient” sampling procedure: with a large value of `n_eff`

and a small value for `se_mean`

we can be relatively sure of the reliability of the mean of the posterior. However, what constitutes a large or small `se_mean`

is harder to define (see Vehtari, Gelman, Simpson, Carpenter, and Bürkner 2019 for a more extensive discussion).^{33}

```
## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 2.76 0.02 1.04 0.72 2.08 2.75 3.43 4.83 2787 1
## sigma 10.43 0.01 0.77 9.07 9.91 10.38 10.91 12.05 2864 1
##
## Samples were drawn using NUTS(diag_e) at Thu Feb 16 17:28:01 2023.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
```

After transforming the `stanfit`

object into a data frame, it’s possible to provide summary plots as the one shown in 10.3. The package `bayesplot`

(Gabry and Mahr 2019) is a wrapper around `ggplot2`

(Wickham, Chang, et al. 2019) and has several convenient functions to plot the samples. Bayesplot functions for posterior summaries start with `mcmc_`

:

There are also several ways to get the samples for other summaries or customized plots, depending on whether we want a list, a data frame, or an array.

```
# extract from rstan is sometimes overwritten by
# a tidyverse version, we make sure that it's right one:
rstan::extract(fit_score) %>%
str()
```

```
## List of 3
## $ mu : num [1:4000(1d)] 1.77 3.71 2.43 1.06 2.09 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ sigma: num [1:4000(1d)] 10.7 11.81 9.93 11.06 11.07 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
## $ lp__ : num [1:4000(1d)] -380 -382 -380 -381 -380 ...
## ..- attr(*, "dimnames")=List of 1
## .. ..$ iterations: NULL
```

```
## 'data.frame': 4000 obs. of 3 variables:
## $ mu : num 1.85 3.36 1.88 3.21 4.29 ...
## $ sigma: num 9.59 8.82 11.51 9.31 11.47 ...
## $ lp__ : num -381 -383 -381 -381 -382 ...
```

```
## num [1:1000, 1:4, 1:3] 1.85 3.36 1.88 3.21 4.29 ...
## - attr(*, "dimnames")=List of 3
## ..$ iterations: NULL
## ..$ chains : chr [1:4] "chain:1" "chain:2" "chain:3" "chain:4"
## ..$ parameters: chr [1:3] "mu" "sigma" "lp__"
```

**Box 10.3 ****An alternative R interface to Stan: cmdstanr**

At the time of writing this, there are two major nuisances with `rstan`

, (i) the R code interfaces directly with C++ creating installation problems in many systems, (ii) `rstan`

releases lag behind Stan language releases considerably preventing the user from taking advantage of the latest features of Stan. The package `cmdstanr`

(https://mc-stan.org/cmdstanr/) is a lightweight interface to Stan for R that solves these problems.
The downside (at the moment of writing this) is that being lightweight some functionality of `rstan`

is lost, such as looking up functions with `lookup()`

, exposing functions with `expose_stan_function()`

, as well as using the fitted model with the `bridgesampling`

package to generate Bayes factors. Furthermore, the package `cmdstanr`

is currently under development and the application programming interface (API) might still change. However, the user interested in an easy (and painless) installation and the latest features of Stan might find it useful.

Once `cmdstanr`

is installed, we can use it as follows:

First create a new `CmdStanModel`

object from a file containing a Stan program using `cmdstan_model()`

```
normal <- system.file("stan_models",
"normal.stan",
package = "bcogsci")
normal_mod <- cmdstan_model(normal)
```

The object `normal_mod`

is an R6 reference object (https://r6.r-lib.org/). This class of object behaves similarly to objects in object oriented programming languages, such as python. Methods are accessed using `$`

(rather than `.`

as in python).

To sample, use the `$sample()`

method. The data argument accepts a list (as we used in `stan()`

from `rstan`

). However, many of the arguments of `$sample`

have different names than the ones used in `stan()`

from the `rstan`

package:

```
lst_score_data <- list(y = Y, N = length(Y))
fit_normal_cmd <- normal_mod$sample(
data = lst_score_data,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000,
iter_sampling = 1000)
```

To show the posterior summary, access the method `$summary()`

of the object `fit_normal_cmd`

.

Access the samples of `fit_normal_cmd`

using `$draws()`

.

The vignette of https://mc-stan.org/cmdstanr/ shows more use cases, and how the samples can be transformed into other formats (data frame, matrix, etc.) together with the package `posterior`

(https://mc-stan.org/cmdstanr/).

## 10.3 Another simple example: Cloze probability with Stan with the binomial likelihood

Let’s fit a Stan model (`binomial_cloze.stan`

) to estimate the cloze probability of a word given its context: that is, what is the probability of an upcoming word given its previous context; the model that is detailed in 2.2 and was fit in 3.1. We want to estimate the cloze probability of *“umbrella”*, \(\theta\), given the following data: *“umbrella”* was answered \(80\) out of \(100\) trials. We assume a binomial distribution as the likelihood function, and \(Beta(a=4,b=4)\) as a prior distribution for the cloze probability.

```
data {
int<lower = 1> N; // Total number of answers
int<lower = 0, upper = N> k; // Number of times umbrella was answered
}
parameters {
// theta is a probability, it has to be constrained between 0 and 1
real<lower = 0, upper = 1> theta;
}
model {
// Prior on theta:
target += beta_lpdf(theta | 4, 4);
// Likelihood:
target += binomial_lpmf(k | N, theta);
}
```

There is only one parameter in this model, cloze probability represented with the parameter `theta`

, which is a real number constrained between \(0\) and \(1\). Another difference between this and the previous example is that the likelihood function ends with `_lpmf`

rather than with `_lpdf`

. This is because Stan differentiates between distributions of continuous variables, i.e, probability density functions (PDFs), and distributions of discrete variables, i.e., probability mass functions (PMFs).

```
lst_cloze_data <- list(k = 80, N = 100)
binomial_cloze <- system.file("stan_models",
"binomial_cloze.stan",
package = "bcogsci")
fit_cloze <- stan(binomial_cloze, data = lst_cloze_data)
```

Print the summary of the posterior distribution of \(\theta\) below, and show its posterior distribution graphically (see Figure 10.4):

```
## mean 2.5% 97.5% n_eff Rhat
## theta 0.78 0.7 0.85 1482 1
```

```
df_fit_cloze <- as.data.frame(fit_cloze)
mcmc_dens(df_fit_cloze, pars = "theta") +
geom_vline(xintercept = mean(df_fit_cloze$theta))
```

## 10.4 Regression models in Stan

In the following sections, we will revisit and expand on some of the examples that were fit with `brms`

in chapter 4.

### 10.4.1 A first linear regression in Stan: Does attentional load affect pupil size?

As in section 4.1, we focus on the effect of cognitive load on one subject’s pupil size with a subset of the data of Wahn et al. (2016). We use the following likelihood and priors. For details about our decision on priors and likelihood, see section 4.1.

\[\begin{equation} \begin{aligned} p\_size_n &\sim \mathit{Normal}(\alpha + c\_load_n \cdot \beta,\sigma) \\ \alpha &\sim \mathit{Normal}(1000, 500) \\ \beta &\sim \mathit{Normal}(0, 100) \\ \sigma &\sim \mathit{Normal}_+(0, 1000) \end{aligned} \end{equation}\]

The Stan model `pupil_model.stan`

follows:

```
data {
int<lower=1> N;
vector[N] p_size;
vector[N] c_load;
}
parameters {
real alpha;
real beta;
real<lower = 0> sigma;
}
model {
// priors:
target += normal_lpdf(alpha | 1000, 500);
target += normal_lpdf(beta | 0, 100);
target += normal_lpdf(sigma | 0, 1000)
- normal_lccdf(0 | 0, 1000);
// likelihood
target += normal_lpdf(p_size | alpha + c_load * beta, sigma);
}
```

Because we are fitting a regression, we use the location (\(\mu\)) of the likelihood function to regress `p_size`

with the following equation `alpha + c_load * beta`

, where both `p_size`

and `c_load`

are vectors defined in the data block. The following line accumulates the log-likelihood of every observation:

`target += normal_lpdf(p_size | alpha + c_load * beta, sigma);`

This is equivalent to and slightly faster than the following lines:

```
for(n in 1:N)
target += normal_lpdf(p_size[n] | alpha + c_load[n] * beta, sigma);
```

A statement that requires some explanation is the following:

```
target += normal_lpdf(sigma | 0, 1000)
- normal_lccdf(0 | 0, 1000);
```

As in our original example in 4.1, we are assuming a truncated normal distribution as a prior for \(\sigma\). Not only are we setting a lower boundary to the parameter with `lower = 0`

, but we are also “correcting” its prior distribution by subtracting `normal_lccdf(0 | 0, 1000)`

, where `lccdf`

stands for log complement of a cumulative distribution function. Once we add a lower boundary, the probability mass under half of the “regular” normal distribution should be one, that is, when we integrate from zero (rather than from minus infinity) to infinity. As discussed in Box 4.1, we need to normalize the PDF by dividing it by the difference of its CDF evaluated in the new boundaries (\(a = 0\) and \(b = \infty\) in our case):

\[\begin{equation} f_{[a,b]}(x) = \frac{f(x)}{F(b) - F(a)} \tag{10.8} \end{equation}\]

where \(f\) is a PDF and \(F\) a CDF.

This equation in log-space is:

\[\begin{equation} log(f_{[a,b]}(x)) = log(f(x)) - log(F(b) - F(a)) \tag{10.9} \end{equation}\]

In Stan \(\log(f(x))\) corresponds to `normal_lpdf(x |...)`

, and `log(F(x))`

to `normal_lcdf(x|...)`

. Because in our example \(b=\infty\), \(F(b) = 1\), we are dealing with the complement of the log CDF evaluated at \(a =0\), \(\log(1 - F(0))\), that is why we use `normal_lccdf(0 | ...)`

(notice the double `c`

; this symbol represents the complement of the CDF).

To be able to fit the model, Stan requires the data to be input as a list: First, load the data and center the dependent variable in a data frame, then create a list, and finally fit the model.

```
ls_pupil <- list(
p_size = df_pupil$p_size,
c_load = df_pupil$c_load,
N = nrow(df_pupil)
)
pupil_model <- system.file("stan_models",
"pupil_model.stan",
package = "bcogsci"
)
fit_pupil <- stan(pupil_model, data = ls_pupil)
```

Check the traceplots (Figure 10.5).

Examine some summaries of the marginal posterior distributions of the parameters of interest:

```
## mean 2.5% 97.5% n_eff Rhat
## alpha 701.9 661.7 742 3693 1
## beta 33.4 10.3 57 4192 1
## sigma 128.4 102.5 162 3419 1
```

Plot the posterior distributions (Figure 10.6).

To determine the probability that the posterior for beta is larger than zero given the model and data, examine the proportion of samples above zero:

`## [1] 0.997`

To generate prior or posterior predictive distributions, we can create our own functions in R with the `purrr`

function `map_dfr`

(or a for-loop) as we did in section 4.2 with the function `lognormal_model_pred()`

. Alternatively, we can use the `generated quantities`

block in our model:

```
data {
int<lower = 1> N;
vector[N] c_load;
int<lower= 0, upper = 1> onlyprior;
vector[N] p_size;
}
parameters {
real alpha;
real beta;
real<lower = 0> sigma;
}
model {
// priors including all constants
target += normal_lpdf(alpha | 1000, 500);
target += normal_lpdf(beta | 0, 100);
target += normal_lpdf(sigma | 0, 1000)
- normal_lccdf(0 | 0, 1000);
if (!onlyprior)
target += normal_lpdf(p_size | alpha + c_load * beta, sigma);
}
generated quantities {
array[N] real p_size_pred;
p_size_pred = normal_rng(alpha + c_load * beta, sigma);
}
```

For most of the probability functions, there is a matching pseudorandom number generator (PRNG) with the suffix `_rng`

. Here we are using the vectorized function `normal_rng`

. Once `p_size_pred`

is declared as an array of size `N`

, the following statement generates \(N\) predictions (for each iteration of the sampler):

`p_size_pred = normal_rng(alpha + c_load * beta, sigma);`

At the moment not all the PRNG are vectorized, but the ones that are only allow for arrays and, confusingly enough, not vectors. We define an array by indicating `array`

, between brackets, the length of each dimension, then the type, and finaly the name of the variable. For example to define an array of real numbers with three dimension of length 6, 7, and 10 we write `array[6, 7, 10] real var`

.^{34} Vectors and matrices are also valid types for an array. See Box 10.4 for more about the difference between arrays and vectors, and other algebra types. We also included a data variable called `onlyprior`

, this is an integer that can only be set to \(1\) (TRUE) or \(0\) (FALSE). When `onlyprior = 1`

, the likelihood is omitted from the model, `p_size`

is ignored, and `p_size_pred`

is the prior predictive distribution. When `onlyprior = 0`

, the likelihood is incorporated in the model (as it is in the original code `pupil_model.stan`

) using `p_size`

, and `p_size_pred`

is the posterior predictive distribution.

If we want posterior predictive distributions, we fit the model to the data and set `onlyprior = 0`

, if we want prior predictive distributions, we sample from the priors and set `onlyprior = 1`

. Then we use `bayesplot`

functions to visualize predictive checks.

For posterior predictive checks, we would do the following:

```
ls_pupil <- list(
onlyprior = 0,
p_size = df_pupil$p_size,
c_load = df_pupil$c_load,
N = nrow(df_pupil)
)
pupil_gen <- system.file("stan_models",
"pupil_gen.stan",
package = "bcogsci")
fit_pupil <- stan(file = pupil_gen, data = ls_pupil)
```

Store the predicted pupil sizes in `yrep_pupil`

. This variable contains an \(N_{samples} \times N_{observations}\) matrix, that is, each row of the matrix is a draw from the posterior predictive distribution, i.e., a vector with one element for each of the data points in y.

`## [1] 4000 41`

Predictive checks functions in `bayesplot`

(starting with `ppc_`

) require a vector with the observations in the first argument and a matrix with the predictive distribution as its second argument. As an example, in Figure 10.7 we use an overlay of densities and we draw only \(50\) elements (that is \(50\) predicted data sets).

For prior predictive distributions, we simply set `onlyprior = 1`

. The observations (`p_size`

) are ignored by the model, but are required by the data block in Stan. If we haven’t collected data yet, we could include a vector of zeros.

```
ls_pupil_prior <- list(
onlyprior = 1,
p_size = df_pupil$p_size,
# or: p_size = rep(0, nrow(df_pupil)),
c_load = df_pupil$c_load,
N = nrow(df_pupil)
)
prior_pupil <- stan(pupil_gen, data = ls_pupil_prior,
control = list(adapt_delta = 0.9))
```

To avoid divergent transitions, increase the `adapt_delta`

parameter’s default value from \(0.8\) to \(0.9\). It is important to highlight that we cannot safely ignore the warnings of the above model, even if we are not fitting data. This is so because in practice one is still sampling a density using Hamiltonian Monte Carlo, and thus the prior sampling process can break in the same ways as the posterior sampling process. Prior predictive distributions as the one shown in Figure 10.8 can be plot with `ppd_dens_overlay()`

(and in general with functions starting with `ppd_`

which don’t require an argument with data).

**Box 10.4 ****Matrix, vector, or array in Stan?**

Stan contains three basic linear algebra types, `vector`

, `row_vector`

, and `matrix`

. But Stan also allows for building arrays of any dimension from any type of element (integer, real, etc.). This means that there are several ways to define one-dimensional N-sized containers of real numbers,

```
array[N] real a;
vector[N] a;
row_vector[N] a;
```

as well as, two-dimensional N1\(\times\)N2-sized containers of real numbers:

```
array[N1, N2] real m;
matrix[N1, N2] m;
array[N1] vector[N2] b;
array[N1] row_vector[N2] b;
```

These distinctions affect either what we can do with these variables, or the speed of our model, and sometimes are interchangeable. Matrix algebra is only defined for (row) vectors and matrices, that is we cannot multiply arrays. The following line requires all the one-dimensional containers (`p_size`

and `c_load`

) to be defined as vectors (or row_vectors):

`vector[N] mu = alpha + c_load * beta;`

Many “vectorized” operations are also valid for arrays, that is, `normal_lpdf`

, accepts (row) vectors (as we did in our code) or arrays as in the next example. There is of course no point in converting a vector to an array as follows, but this shows that Stan allows both type of one-dimensional containers.

```
array[N] real mu = to_array_1d(alpha + c_load * beta);
target += normal_lpdf(p_size | mu, sigma);
```

By contrast, the outcome of “vectorized” pseudorandom number generator (`_rng`

) functions can only be stored in an array. The following example shows the only way to vectorize this type of function:

```
array[N] real p_size_pred = normal_rng(alpha + c_load * beta,
sigma);
```

Alternatively, one can always use a for-loop, and it won’t matter if `p_size_pred`

is an array or a vector:

```
vector[N] p_size_pred;
for(n in 1:N)
p_size_pred[n] = normal_rng(alpha + c_load[n] * beta, sigma);
```

See also Stan’s manual section on matrices, vector, and arrays (Stan Development Team 2021).

### 10.4.2 Interactions in Stan: Does attentional load interact with trial number affecting pupil size?

We’ll expand the previous model to also include the effect of (centered) trial and its interaction with cognitive load on one subject’s pupil size. Our new likelihood will look as follows:

\[\begin{equation} p\_size_n \sim \mathit{Normal}(\alpha + c\_load_n \cdot \beta_1 + c\_trial \cdot \beta_2 + c\_load \cdot c\_trial \cdot \beta_3, \sigma) \end{equation}\]

Define priors for all the new \(\beta\)s. Since we don’t have more information about the new predictors, they are sampled from identical prior distributions:

\[\begin{equation} \begin{aligned} \alpha &\sim \mathit{Normal}(1000, 500) \\ \beta_1 &\sim \mathit{Normal}(0, 100) \\ \beta_2 &\sim \mathit{Normal}(0, 100) \\ \beta_3 &\sim \mathit{Normal}(0, 100) \\ \sigma &\sim \mathit{Normal}_+(0, 1000) \end{aligned} \end{equation}\]

The following Stan model, `pupil_int1.stan`

, is the direct translation of the new priors and likelihood.

```
data {
int<lower = 1> N;
vector[N] c_load;
vector[N] c_trial;
vector[N] p_size;
}
parameters {
real alpha;
real beta1;
real beta2;
real beta3;
real<lower = 0> sigma;
}
model {
// priors including all constants
target += normal_lpdf(alpha | 1000, 500);
target += normal_lpdf(beta1 | 0, 100);
target += normal_lpdf(beta2 | 0, 100);
target += normal_lpdf(beta3 | 0, 100);
target += normal_lpdf(sigma | 0, 1000)
- normal_lccdf(0 | 0, 1000);
target += normal_lpdf(p_size | alpha + c_load * beta1 +
c_trial * beta2 +
c_load .* c_trial * beta3, sigma);
}
```

When there are matrices or vectors involved, `*`

indicates matrix multiplication whereas `.*`

indicates element-wise multiplication; in R `%*%`

indicates matrix multiplication whereas `*`

indicates element-wise multiplication.

There is, however, an alternative notation that can simplify our code. In the following likelihood, \(p\_size\) is a vector of N observations (in this case 41), \(X\) is the model matrix with a dimension of \(N \times N_{pred}\) (in this case \(41 \times 3\)), and \(\beta\) a vector of \(N_{pred}\) (in this case, 3) rows. Assuming that \(\beta\) is a vector, we indicate with one line that each parameter \(\beta_n\) is sampled from identical prior distributions.

\[\begin{equation} \begin{aligned} p\_size &\sim \mathit{Normal}(\alpha + X \cdot \beta,\sigma)\\ \beta &\sim \mathit{Normal}(0, 100) \\ \sigma &\sim \mathit{Normal}_+(0, 1000) \end{aligned} \end{equation}\]

The translation into Stan code is the following:

```
data {
int<lower = 1> N;
int<lower = 0> K; // number of predictors
matrix[N, K] X; // model matrix
vector[N] p_size;
}
parameters {
real alpha;
vector[K] beta;
real<lower = 0> sigma;
}
model {
// priors including all constants
target += normal_lpdf(alpha | 1000, 500);
target += normal_lpdf(beta | 0, 100);
target += normal_lpdf(sigma | 0, 1000)
- normal_lccdf(0 | 0, 1000);
target += normal_lpdf(p_size | alpha + X * beta, sigma);
}
```

For some likelihood functions, Stan provides a more efficient implementation of the linear regression than the one manually written in the previous code. It’s critical to understand that, in general, a more efficient implementation should not only be faster, but should also achieve the same number of effective samples (or more) than a less efficient implementation (and should also show convergence). In this case, we can achieve that using `_glm`

functions. We can replace the last line with the following statement (the order of the arguments is important):^{35}

`target += normal_id_glm_lpdf(p_size | X, alpha, beta, sigma);`

The most optimized model, `pupil_int.stan`

, includes this last statement. We prepare the data as follows: First create a centered version of trial, `c_trial`

and load `c_load`

, then use the function `model.matrix`

to create the `X`

matrix that contains in each column our predictors and omits the intercept with `0 +`

.

```
df_pupil <- df_pupil %>%
mutate(
c_trial = trial - mean(trial),
c_load = load - mean(load)
)
X <- model.matrix(~ 0 + c_load * c_trial, df_pupil)
ls_pupil_X <- list(
p_size = df_pupil$p_size,
X = X,
K = ncol(X),
N = nrow(df_pupil)
)
```

```
pupil_int <- system.file("stan_models",
"pupil_int.stan",
package = "bcogsci"
)
fit_pupil_int <- stan(pupil_int, data = ls_pupil_X)
```

```
## mean 2.5% 97.5% n_eff Rhat
## alpha 699.50 667.05 731.87 4326 1
## beta[1] 31.27 12.49 50.33 5018 1
## beta[2] -5.81 -8.48 -3.14 4356 1
## beta[3] -1.81 -3.43 -0.20 5080 1
## sigma 104.39 82.74 132.66 3682 1
```

In 10.9, we plot here the 95% CrI of the parameters of interest. We use `regex_pars`

, rather than `pars`

, because we want to capture `beta[1]`

, `beta[2]`

, and `beta[3]`

; `regex_pars`

use regular expressions to select the parameters (for information about regular expressions in R see https://stat.ethz.ch/R-manual/R-devel/library/base/html/regex.html)

```
df_fit_pupil_int <- as.data.frame(fit_pupil_int)
mcmc_intervals(fit_pupil_int,
regex_pars = "beta",
prob_outer = .95,
prob = .8,
point_est = "mean"
)
```

### 10.4.3 Logistic regression in Stan: Does set size and trial affect free recall?

We revisit and expand on the analysis presented in 4.3 of a subset of the data of Oberauer (2019). In this example, we will investigate whether the length of a list and trial number affect the probability of correctly recalling a word.

As in section 4.3, we assume a Bernoulli likelihood with a logit link function, and the following priors (recall that the logistic function is the inverse of the logit).

\[\begin{equation} \begin{aligned} correct_n &\sim \mathit{Bernoulli}( \mathit{logistic}(\alpha + X \cdot \beta))\\ \alpha &\sim \mathit{Normal}(0, 1.5) \\ \beta &\sim \mathit{Normal}(0, 0.1) \end{aligned} \end{equation}\]

Where \(\beta\) is a vector of size \(K = 2\), \(\{\beta_0, \beta_1\}\). Below in `recall.stan`

we present the most efficient way to code this in Stan.

```
data {
int<lower = 1> N;
int<lower=0> K; // number of predictors
matrix[N, K] X; // model matrix
array[N] int correct;
}
parameters {
real alpha;
vector[K] beta;
}
model {
// priors including all constants
target += normal_lpdf(alpha | 0, 1.5);
target += normal_lpdf(beta | 0, .1);
target += bernoulli_logit_glm_lpmf(correct | X, alpha, beta);
}
```

The dependent variable, `correct`

, is an array of integers rather than a vector; this is because vectors are always composed of real numbers, but the Bernoulli likelihood only accepts the integers \(1\) or \(0\). As in the previous example, we are taking advantage of the `_glm`

functions. A less efficient but more transparent option would be to replace the last statement with:

`target += bernoulli_logit_lpmf(correct | alpha + X * beta);`

We might want to use `bernoulli_logit_lpmf`

if we want to define a non-linear relationship between the predictors that are outside the generalized linear model framework. One example would be the following:

`target += bernoulli_logit_lpmf(correct| alpha + exp(X * beta));`

Another more flexible possibility when we want to indicate a Bernoulli likelihood is to use `bernoulli_lpmf`

and add the link manually. The last statement of `recall.stan`

would become the following:

`target += bernoulli_lpmf(correct| inv_logit(alpha + X * beta));`

The function `bernoulli_lpmf`

can be useful if one wants to try other link functions; see exercise 10.4.

Finally, the most transparent form (but less efficient) would be the following for-loop:

```
for (n in 1:N)
target += bernoulli_lpmf(correct[n] | inv_logit(alpha + X[n] * beta));
```

To fit the model as `recall.stan`

, prepare the data by centering the trial number first:

```
df_recall <- df_recall %>%
mutate(c_set_size = set_size - mean(set_size),
c_trial = trial - mean(trial))
```

Next, we create the model matrix, \(X\), and prepare the data as a list. As in section 10.4.2, we exclude the intercept from the matrix `X`

using `0 +...`

. This is because the Stan code that we are using already takes into account that the first column in the model matrix is going to be a vector of ones.

```
X <- model.matrix(~ 0 + c_set_size * c_trial, df_recall)
ls_recall <- list(
correct = df_recall$correct,
X = X,
K = ncol(X),
N = nrow(df_recall)
)
```

```
recall <- system.file("stan_models",
"recall.stan",
package = "bcogsci")
fit_recall <- stan(recall, data = ls_recall)
```

After fitting the model we can print and plot summaries of the posterior distribution.

```
## mean 2.5% 97.5% n_eff Rhat
## alpha 1.99 1.39 2.65 3516 1
## beta[1] -0.19 -0.35 -0.02 4014 1
## beta[2] -0.02 -0.09 0.05 3366 1
## beta[3] 0.00 -0.03 0.04 3886 1
```

In Figure 10.10, we plot the 95% CrI of the parameters of interest.

```
df_fit_recall <- as.data.frame(fit_recall)
mcmc_intervals(df_fit_recall,
regex_pars = "beta",
prob_outer = .95,
prob = .8,
point_est = "mean")
```

As we did in 4.3.4, we might want to communicate the posterior in proportions rather than in log-odds (as seen in the parameters `beta`

). We can do this in R manipulating the data frame `df_fit_recall`

, or extracting the samples with `extract(fit_recall)`

. Another alternative presented here is by using the generated quantities block. To make the code more compact we declare the type of each variable and store its content in the same line in `recall_prop.stan`

.

```
generated quantities {
real average_accuracy = inv_logit(alpha);
vector[K] change_acc = inv_logit(alpha) - inv_logit(alpha - beta);
}
```

Recall that due to the non-linearity of the scale, the effects depend on the average accuracy, and the set size and trial that we start from (in this case we are examining the change of one unit from the average set size and the average trial).

```
recall_prop <- system.file("stan_models",
"recall_prop.stan",
package = "bcogsci")
fit_recall <- stan(recall_prop, data = ls_recall)
```

The plot in Figure 10.11 now shows how the average accuracy deteriorates when the subject is exposed to a set size larger than the average by one, a trial further than the middle one, and the interaction of both.

```
df_fit_recall <- as.data.frame(fit_recall) %>%
rename(
set_size = `change_acc[1]`,
trial = `change_acc[2]`,
interaction = `change_acc[3]`
)
mcmc_intervals(df_fit_recall,
pars = c("set_size", "trial", "interaction"),
prob_outer = .95,
prob = .8,
point_est = "mean"
) +
xlab("Change in accuracy")
```

The plot in Figure 10.11 shows that our model is estimating that by increasing the set size by one unit, the recall accuracy of the single subject is deteriorates by 2%. In contrast, there is hardly any trial effect or interaction between trial and set size.

## 10.5 Summary

This chapter introduced basic Stan syntax for fitting some standard linear models. Example code covered the normal, binomial, Bernoulli, and log-normal likelihoods. We also saw how to express regression models using the matrix model in Stan syntax.

## 10.6 Further reading

For further reading on the Hamiltonian Monte Carlo algorithm, see the rather technical review of Betancourt (2017), or the more conceptual introduction provided by Monnahan, Thorson, and Branch (2017). A useful article with example R code is Neal (2011). A detailed walk-through on its implementation is also provided in Chapter 41 of MacKay (2003). The Stan documentation (Stan Development Team 2021), consisting of a User’s Guide and the Language Reference Manual are important starting points for going deeper into Stan programming: https://mc-stan.org/users/documentation/.

## 10.7 Exercises

**Exercise 10.1 **A very simple model.

In this exercise we revisit the model from 3.2.1. Assume the following:

- There is a true underlying time, \(\mu\), that the subject needs to press the space bar.
- There is some noise in this process.
- The noise is normally distributed (this assumption is questionable given that response times are generally skewed; we fix this assumption later).

That is the likelihood for each observation \(n\) will be:

\[\begin{equation} \begin{aligned} t_n \sim \mathit{Normal}(\mu, \sigma) \end{aligned} \end{equation}\]

- Decide on appropriate priors and fit this model in Stan. Data can be found in
`df_spacebar`

. - Change the likelihood for a log-normal distribution and change the priors. Fit the model in Stan.

**Exercise 10.2 **Incorrect Stan model.

We want to fit both response times and accuracy with the same model. We simulate the data as follows:

We build the following model:

```
data {
int<lower = 1> N;
vector[N] rt;
array[N] int correct;
}
parameters {
real<lower = 0> sigma;
real theta;
}
model {
target += normal_lpdf(mu | 0, 20);
target += lognormal_lpdf(sigma | 3, 1)
for(n in 1:N)
target += lognormal_lpdf(rt[n] | mu, sigma);
target += bernoulli_lpdf(correct[n] | theta);
}
```

Why does this model not work?

```
ls_sim <- list(
rt = df_sim$rt,
correct = df_sim$correct
)
incorrect <- system.file("stan_models",
"incorrect.stan",
package = "bcogsci")
fit_sim <- stan(incorrect, data = ls_sim)
```

```
## Error in stanc(file = file, model_code = model_code, model_name = model_name, : 0
##
## Syntax error in 'string', line 13, column 2 to column 5, parsing error:
##
## Ill-formed expression. Expression followed by ";" expected after "target +=".
```

Try to make it run. (Hint: There are several problems.)

**Exercise 10.3 **Using Stan documentation.

Edit the simple example with Stan from section 10.2, and replace the normal distribution with a skew normal distribution. (Don’t forget to add a prior to the new parameter, and check the Stan documentation or a statistics textbook for more information about the distribution).

Fit the following data:

Does the estimate of the new parameter make sense?

**Exercise 10.4 **The probit link function as an alternative to the logit function.

The probit link function is the inverse of the CDF of the standard normal distribution (\(Normal(0,1)\)). Since the CDF of the standard normal is usually denoted with the Greek letter \(\Phi\) (Phi), the probit is denoted as \(\Phi^{-1}\). Refit the model presented in 10.4.3 changing the logit link function for the probit link (that is transforming the regression to a constrained space using `Phi()`

in Stan).

You will probably see the following as the model runs; this is because the probit link is less numerically stable (i.e., under- and overflows) than the logit link *in Stan*. Don’t worry, it is good enough for this exercise.

```
Rejecting initial value:
Log probability evaluates to log(0), i.e. negative infinity.
Stan can't start sampling from this initial value.
```

- Do the results of the coefficients \(\alpha\) and \(\beta\) change?
- Do the results in probability space change?

**Exercise 10.5 **Examining the position of the queued word on recall.

Refit the model presented in section 10.4.3 and examine whether, set size, trial effects, the position of the queued word (`tested`

in the data set), and their interaction affect free recall. (Tip: You can do this exercise without changing the Stan code.).

How does the accuracy change from position one to position two?

**Exercise 10.6 **The conjunction fallacy.

Paolacci, Chandler, and Ipeirotis (2010) examined whether the results of some classic experiments differ between a university pool population and subjects recruited from Mechanical Turk. We’ll examine whether the results of the conjunction fallacy experiment (or Linda problem: Tversky and Kahneman 1983) are replicated for both groups.

```
## # A tibble: 268 × 2
## source answer
## <chr> <int>
## 1 mturk 1
## 2 mturk 1
## 3 mturk 1
## # … with 265 more rows
```

The conjunction fallacy shows that people often fail to regard a combination of events as less probable than a single event in the combination (Tversky and Kahneman 1983):

*Linda is 31 years old, single, outspoken, and very bright. She majored in philosophy. As a student, she was deeply concerned with issues of discrimination and social justice, and also participated in anti-nuclear demonstrations.*

*Which is more probable?*

*Linda is a bank teller.**Linda is a bank teller and is active in the feminist movement.*

The majority of those asked chose option b even though it’s less probable (\(\Pr(a \land b)\leq \Pr(b)\). The data set is named `df_fallacy`

and it indicates with `0`

option “a” and with `1`

option `b`

.
Fit a logistic regression in Stan and report:

- The estimated overall probability of answering (b) ignoring the group.
- The estimated overall probability of answering (b) for each group.

### References

Betancourt, Michael J. 2016. “Identifying the Optimal Integration Time in Hamiltonian Monte Carlo.”

Betancourt, Michael J. 2017. “A Conceptual Introduction to Hamiltonian Monte Carlo.”

Carpenter, Bob, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael J. Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “Stan: A Probabilistic Programming Language.” *Journal of Statistical Software* 76 (1). Columbia Univ., New York, NY (United States); Harvard Univ., Cambridge, MA (United States).

Gabry, Jonah, and Tristan Mahr. 2019. *bayesplot: Plotting for Bayesian Models*. https://CRAN.R-project.org/package=bayesplot.

Guo, Jiqiang, Jonah Gabry, and Ben Goodrich. 2019. *rstan: R Interface to Stan*. https://CRAN.R-project.org/package=rstan.

Hoffman, Matthew D., and Andrew Gelman. 2014. “The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.” *Journal of Machine Learning Research* 15 (1): 1593–1623. http://dl.acm.org/citation.cfm?id=2627435.2638586.

Lunn, David, Chris Jackson, David J Spiegelhalter, Nicky Best, and Andrew Thomas. 2012. *The BUGS Book: A Practical Introduction to Bayesian Analysis*. Vol. 98. CRC Press.

Lunn, D.J., A. Thomas, N. Best, and D. Spiegelhalter. 2000. “WinBUGS-A Bayesian Modelling Framework: Concepts, Structure, and Extensibility.” *Statistics and Computing* 10 (4). Springer: 325–37.

MacKay, David JC. 2003. *Information Theory, Inference and Learning Algorithms*. Cambridge, UK: Cambridge University Press.

Monnahan, Cole C., James T. Thorson, and Trevor A. Branch. 2017. “Faster Estimation of Bayesian Models in Ecology Using Hamiltonian Monte Carlo.” Edited by Robert B. O’Hara. *Methods in Ecology and Evolution* 8 (3): 339–48. https://doi.org/10.1111/2041-210X.12681.

Neal, Radford M. 2011. “MCMC Using Hamiltonian Dynamics.” In *Handbook of Markov Chain Monte Carlo*, edited by Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Taylor & Francis. https://doi.org/10.1201/b10905-10.

Oberauer, Klaus. 2019. “Working Memory Capacity Limits Memory for Bindings.” *Journal of Cognition* 2 (1): 40. https://doi.org/10.5334/joc.86.

Paolacci, Gabriele, Jesse Chandler, and Panagiotis G Ipeirotis. 2010. “Running Experiments on Amazon Mechanical Turk.” *Judgment and Decision Making* 5 (5): 411–19.

Plummer, Martin. 2016. “JAGS Version 4.2.0 User Manual.”

R Core Team. 2019. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Stan Development Team. 2021. “Stan Modeling Language Users Guide and Reference Manual, Version 2.27.” https://mc-stan.org.

Tversky, Amos, and Daniel Kahneman. 1983. “Extensional Versus Intuitive Reasoning: The Conjunction Fallacy in Probability Judgment.” *Psychological Review* 90 (4). American Psychological Association: 293.

Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2019. “Rank-Normalization, Folding, and Localization: An Improved \(\widehat{R}\) for Assessing Convergence of MCMC.”

Wahn, Basil, Daniel P. Ferris, W. David Hairston, and Peter König. 2016. “Pupil Sizes Scale with Attentional Load and Task Experience in a Multiple Object Tracking Task.” *PLOS ONE* 11 (12): e0168087. https://doi.org/10.1371/journal.pone.0168087.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, and Hiroaki Yutani. 2019. *ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics*. https://CRAN.R-project.org/package=ggplot2.

In the specific case of a model with two parameters, e.g., \(<\mu,\sigma>\), the physical analogy works quite well: The \(<x,y>\) coordinates of the particle would be determined by \(<\mu,\sigma>\), and its \(z\) coordinate would be established by the height of the unnormalized posterior.↩

Incidentally, this \(\log(up(\dot))\) is the variable

`target`

in a Stan model and`lp__`

in its output; see Box 10.1.↩In these output, there are some types that are new to the R user (but they are also used in C++):

`reals`

indicates that any of`real`

,`real[]`

,`vector`

, or`row_vector`

. A return type`R`

with an input type`T`

indicates that the type of the output of the function is the same as type of the argument.↩We simplify the output of

`print`

in the text after this call, by actually calling`summary(fit, pars = pars, probs = c(0.025, 0.975))$summary`

.↩The notation for arrays has changed in recent versions, the previous notation would have been

`real[6, 7, 10] var`

.↩An extra boost in efficiency can be obtained in regular regressions where

`X`

and`Y`

are data (rather than parameters as in cases of missing data or measurement error), since this function can be executed on a GPU.↩