# Statistical Rethinking Simple Models Crashing

I’m working through Statistical rethinking right now in order to learn Stan, and having issues where fitting most MCMC models crash my R session. In some cases they’ll fit and the session crashes after running a few more code blocks, in some it crashes when fitting. I assume this is a memory issue, but I can’t imagine straightforward models would crash my session like this. I’m using a 2015 dual-core 8gb macbook pro. Is this common?

This is a persisting problem, but the most recent example is model 14.1:

``````m14.1 <- ulam(
alist(
wait ~ normal( mu , sigma ),
mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
c(a_cafe,b_cafe)[cafe] ~ multi_normal( c(a,b) , Rho , sigma_cafe ),
a ~ normal(5,2),
b ~ normal(-1,0.5),
sigma_cafe ~ exponential(1),
sigma ~ exponential(1),
Rho ~ lkj_corr(2)
) , data=d , chains=2 , cores=1)
``````

This crashes the session on fitting. I also tried fitting with pure Stan, not the ulam overlay:

`````````{r}
try01 <- as.character("
data{
vector wait;
int afternoon;
int cafe;
}
parameters{
vector b_cafe;
vector a_cafe;
real a;
real b;
vector<lower=0> sigma_cafe;
real<lower=0> sigma;
corr_matrix Rho;
}
model{
vector mu;
Rho ~ lkj_corr( 2 );
sigma ~ exponential( 1 );
sigma_cafe ~ exponential( 1 );
b ~ normal( -1 , 0.5 );
a ~ normal( 5 , 2 );
{
vector YY;
vector MU;
MU = [ a , b ]';
for ( j in 1:20 ) YY[j] = [ a_cafe[j] , b_cafe[j] ]';
YY ~ multi_normal( MU , quad_form_diag(Rho , sigma_cafe) );
}
for ( i in 1:200 ) {
mu[i] = a_cafe[cafe[i]] + b_cafe[cafe[i]] * afternoon[i];
}
wait ~ normal( mu , sigma );
}
")
m14.1 <- stan( model_code = try01, data=d, chains=2)
``````

which also crashes for me. I’ve tried running from Rstudio and from command line. Full disclosure, I’m doing it within the context of an Rmarkdown script, so I’m running as:

``````Rscript -e "rmarkdown::render('ch14.Rmd', clean=TRUE)
``````

Relevant session info:

`````` version  R version 4.0.1 (2020-06-06)
os       macOS Catalina 10.15.4
system   x86_64, darwin19.5.0
ui       RStudio

rethinking_2.01
rstan        * 2.21.2    2020-07-27  CRAN (R 4.0.1)
``````

Welcome to the Stan community!

The below works for me with `rethinking v2.12`, `StanHeaders v2.21.0.5`, and `rstan v2.21.2` on `OS X` (I’ve cut’n pasted the code from the book, i.e., edition 2). It takes some time to compile but then it samples for, like, a second. Like a warm knife through butter :)

``````a <- 3.5            # average morning wait time
b <- (-1)           # average difference afternoon wait time
sigma_a <- 1        # std dev in intercepts
sigma_b <- 0.5      # std dev in slopes
rho <- (-0.7)       # correlation between intercepts and slopes

Mu <- c( a , b )

## R code 14.3
cov_ab <- sigma_a*sigma_b*rho
Sigma <- matrix( c(sigma_a^2,cov_ab,cov_ab,sigma_b^2) , ncol=2 )

## R code 14.4
matrix( c(1,2,3,4) , nrow=2 , ncol=2 )

## R code 14.5
sigmas <- c(sigma_a,sigma_b) # standard deviations
Rho <- matrix( c(1,rho,rho,1) , nrow=2 ) # correlation matrix

# now matrix multiply to get covariance matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)

## R code 14.6
N_cafes <- 20

## R code 14.7
library(MASS)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm( N_cafes , Mu , Sigma )

## R code 14.8
a_cafe <- vary_effects[,1]
b_cafe <- vary_effects[,2]

set.seed(22)
N_visits <- 10
afternoon <- rep(0:1,N_visits*N_cafes/2)
cafe_id <- rep( 1:N_cafes , each=N_visits )
mu <- a_cafe[cafe_id] + b_cafe[cafe_id]*afternoon
sigma <- 0.5  # std dev within cafes
wait <- rnorm( N_visits*N_cafes , mu , sigma )
d <- data.frame( cafe=cafe_id , afternoon=afternoon , wait=wait )

set.seed(867530)
m14.1 <- ulam(
alist(
wait ~ normal( mu , sigma ),
mu <- a_cafe[cafe] + b_cafe[cafe]*afternoon,
c(a_cafe,b_cafe)[cafe] ~ multi_normal( c(a,b) , Rho , sigma_cafe ),
a ~ normal(5,2),
b ~ normal(-1,0.5),
sigma_cafe ~ exponential(1),
sigma ~ exponential(1),
Rho ~ lkj_corr(2)
) , data=d , chains=4 , cores=4 )
``````
3 Likes

Yep, looks like you’re right… It appears all of this runs fine if typed straight into an R console, or in a standard R script, it just dies in the context of Rmd, for whatever reason. Frustrating, I’m surprised the execution is different on memory, but nice to have a workable solution.

Thanks for checking! Sorry for the static with this being ~trivial.

1 Like

Perhaps using the latest support for `cmdstanr` in `rethinking` could be an option?

``````> devtools::install_github("stan-dev/cmdstanr")
> install_cmdstan()
``````

restart, and then run `ulam()` with option `cmdstan=TRUE`.

2 Likes

Just wanted to come back and thank you for this! I was bashing my head against my computer today running into the problem again, and finally got cmdstan running and everything seems to be working like a charm.

Thanks so much!

2 Likes