C stack error occurs when using rstan in a for-loop or functions of purrr
I have the very same issue as the thread starter had. In my case, I’m using brms::brm() with rstan backend in a for-loop (actually purrr::walk2). Moreover, although @andrjohns suggests assigning a unique name for different models is helpful, my case may be incurable by that method.
MWE
What I intend to do
I computed a lot of models using brm() with cmdstanr backend with different prior settings for a parameter. Now, I want to compute the log marginal likelihood of each model. To obtain the likelihood, I need to compile brm() with rstan backend first, as discussed on here. But compiling brm() with rstan backend multiple times cause the C stack error, even if (1) rstan appears in a loop environment and (2) every rstan model has no name.
Attach packages
library(lme4) ## For the dataset called sleepstudy
library(tidyverse)
library(magrittr)
library(brms)
library(cmdstanr)
library(rstan)
## Parallelize the chains using all the cores:
options(mc.cores = parallel::detectCores())
## Save compiled models:
rstan_options(auto_write = TRUE)
Set priors
I provided here 9 different combinations of the prior mean and sd for a explanatory variable called Days.
priors <- data.frame(
prior_mean = prior_mean <- c(-5, 0, 5) %>%
rep(., each = length(c(0.01, 0.1, 1))),
prior_sd = prior_sd <- c(0.01, 0.1, 1) %>%
rep(., length(c(-5, 0, 5))),
id = case_when(
prior_mean == 0 ~ paste0(prior_mean, "_", prior_sd),
prior_mean < 0 ~ paste0("m", abs(prior_mean), "_", prior_sd),
prior_mean > 0 ~ paste0("p", prior_mean, "_", prior_sd)
)
)
Sample the model with various parameter settings
Since I want to compute the models with cmdstanr using 9 sets of priors, I made an original function called computation, which takes the prior mean and sd and compose the brm(). Then I compute the 9 models in purrr::walk2, which is one of dplyr's loop function.
computation = function(prior_mean, prior_sd) {
file_name = case_when(
prior_mean == 0 ~ paste0("D:/test/c/c_", prior_mean, "_", prior_sd, "_fit.Rds"),
prior_mean < 0 ~ paste0("D:/test/c/c_m", abs(prior_mean), "_", prior_sd, "_fit.Rds"),
prior_mean > 0 ~ paste0("D:/test/c/c_p", prior_mean, "_", prior_sd, "_fit.Rds")
)
# pass formula, prior, data, options, etc.
brms_object = brm(
Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
set_prior(
paste0(
"normal(",
prior_mean,
", ",
prior_sd,
")"
),
class = "b",
coef = "Days"
),
prior(normal(0, 1), class = sd),
prior(lkj(2), class = cor)
),
warmup = 2000,
iter = 10000,
chains = 4,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = TRUE),
backend = "cmdstanr"
)$fit
saveRDS(brms_object, file_name)
}
walk2(
priors$prior_mean,
priors$prior_sd,
computation
)
So far, I had no problem to sample the models with cmdstanr backend and to save the results as .Rds file.
Compute log marginal likelihood with various prior settings
By creating a original function ml, I set the brms models with rstan backend as iter = 0 (and chain = 0) and pass them bridgesampling::bridge_sampler(..., stanfit_model = .) to calculate log marginal likelihood in different prior parameter settings. However, this flow hangs in a loop environment, i.e. walk2(), at third or fourth iteration. Note that I do not assign any object name for the stanfits.
ml = function(prior_mean, prior_sd) {
brm(
Reaction ~ 0 + Intercept + Days + (0 + Intercept + Days|Subject),
data = sleepstudy,
family = "normal",
prior =
c(
prior(normal(0, 1), class = b, coef = Intercept),
set_prior(
paste0(
"normal(",
prior_mean,
", ",
prior_sd,
")"
),
class = "b",
coef = "Days"
),
prior(normal(0, 1), class = sd),
prior(lkj(2), class = cor)
),
iter = 0,
chains = 0,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = TRUE),
backend = "rstan"
)$fit %>%
bridgesampling::bridge_sampler(
samples = case_when(
prior_mean == 0 ~ paste0("D:/test/c/c_", prior_mean, "_", prior_sd, "_fit.Rds"),
prior_mean < 0 ~ paste0("D:/test/c/c_m", abs(prior_mean), "_", prior_sd, "_fit.Rds"),
prior_mean > 0 ~ paste0("D:/test/c/c_p", prior_mean, "_", prior_sd, "_fit.Rds")
) %>% readRDS(),
stanfit_model = .
) %>%
saveRDS(
.,
case_when(
prior_mean == 0 ~ paste0("D:/test/ml/ml_r1_", prior_mean, "_", prior_sd, ".Rds"),
prior_mean < 0 ~ paste0("D:/test/ml/ml_r1_m", abs(prior_mean), "_", prior_sd, ".Rds"),
prior_mean > 0 ~ paste0("D:/test/ml/ml_r1_p", prior_mean, "_", prior_sd, ".Rds")
)
)
}
purrr::walk2(
priors$prior_mean,
priors$prior_sd,
ml
) ##### The loop hangs at fourth iteration
Session Info
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
rstan_2.21.2
StanHeaders_2.21.0-7
cmdstanr_0.3.0.9000
brms_2.14.4
