Statistical Rethinking Simple Models Crashing

Interfaces RStan
1

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[200] wait;
  int afternoon[200];
  int cafe[200];
}
parameters{
  vector[20] b_cafe;
  vector[20] a_cafe;
  real a;
  real b;
  vector<lower=0>[2] sigma_cafe;
  real<lower=0> sigma;
  corr_matrix[2] Rho;
}
model{
  vector[200] mu;
  Rho ~ lkj_corr( 2 );
  sigma ~ exponential( 1 );
  sigma_cafe ~ exponential( 1 );
  b ~ normal( -1 , 0.5 );
  a ~ normal( 5 , 2 );
  {
    vector[2] YY[20];
    vector[2] 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 [1] CRAN (R 4.0.1)
2

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
3

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
4

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
5

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