Weibull model that crashes RStudio

Modeling
1

Hi
I’m trying to run a model but it crashes RStudio with a message I can’t read in full (see below) and also don’t understand. The model is probably wrong, but based on this error message(s) I can’t figure out what I did wrong.

Data simulation

black<-round(rweibull(1000,1.1,65))
other<-round(rweibull(1000,1.1,55))

d<-data.frame(days_to_event=c(black,other),color_id=c(rep(1,1000),rep(2,1000)))
d<- d %>%  mutate(adopted = if_else(days_to_event >= 200, 0, 1))


dat <- list(
  N=nrow(d),
  days_to_event = d$days_to_event,
  color_id = d$color_id ,
  adopted = d$adopted
)

Stan model

data{
    int N;
    int adopted[N];
    vector[N] days_to_event;
    int color_id[N];

}
parameters{
    vector[2] a;
    real<lower=0> k;
    real beta;
}
transformed parameters{
  vector[N] lambda;
         for ( i in 1:N ) {
         lambda[i] = exp((a[color_id[i]]*beta)*k);
    }
}
model{
   a~exponential(0.4);
   beta~exponential(0.4);
   k~exponential(1);
    for ( i in 1:N ) 
        if ( adopted[i] == 0 ) target += weibull_lccdf(days_to_event[i] | k,lambda[i]);
    for ( i in 1:N )
        if ( adopted[i] == 1 ) days_to_event[i] ~ weibull( k,lambda[i]);
}
generated quantities{
  real pred[N];
   for ( i in 1:N ) 
   if ( adopted[i] == 1 ) pred[i] = weibull_rng( k,lambda[i] );
   for ( i in 1:N ) 
   if ( adopted[i] == 0 ) pred[i] = weibull_lccdf(days_to_event[i] | k,lambda[i]);
}

Run the model:

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
parallel:::setDefaultClusterOptions(setup_strategy = "sequential")
fit<-stan("model_weibull.stan",data=dat)

RStudio crashes shortly after with this message. Unfortunately once it crashes I can’t read the full error message so I had to take a screenshot:

weib
weib1031×555 127 KB

2

Rstudio has been having a lot of issues with Stan lately. Maybe post your code/experience over in this thread so the Rstudio folks can help?

3

Take out the rounding in black and other.

@rok_cesnovar I think this is a bug since the declaration of a vector should promote the data to reals.

Other than that you have another issue. The parameters a and beta should be constrained to the positive reals since an exponential prior is only defined for positive reals. However, this is unlikely what you want since you’re creating lambda using exp, which means lambda > 0.

Efficiency stuff. You’re doubling loops. Instead of having to loop through N twice just do

  for ( i in 1:N ) {
        if ( adopted[i] == 0 ) target += weibull_lccdf( days_to_event[i] | k, lambda[i] );
          else days_to_event[i] ~ weibull( k, lambda[i] );
     }

The parameter lambda can be declared in one line as the exp() function is vectorized.

transformed parameters{
  vector[N] lambda = exp( (a[color_id] * beta) * k);
}

You can do better by removing the loops completely. To do this you need to create indexes of when adopted = 0 and 1.

 target += weibull_lccdf( days_to_event[index_zero] | k, lambda[index_zero] );
 days_to_event[index_one] ~ weibull( k, lambda[index_one] );
2 Likes
4

@spinkney I run the corrected model but I’m getting the same crash. I also tried to run it directly from R (Ubuntu 20.04) and R crashes with no error messages. Maybe I missed something?

data{
    int N;
    int adopted[N];
    vector[N] days_to_event;
    int color_id[N];

}
parameters{
    vector[2] a;
    real<lower=0> k;
    real<lower=0> beta;
}
transformed parameters{
vector[N] lambda = exp( (a[color_id] * beta) * k);
}

model{
   a~exponential(0.4);
   beta~exponential(0.4);
   k~exponential(1);
   for ( i in 1:N ) {
        if (adopted[i] == 0 ) target += weibull_lccdf( days_to_event[i] | k, lambda[i] );
        else days_to_event[i] ~ weibull( k, lambda[i] );
     }
}
generated quantities{
  real pred[N];
   for ( i in 1:N ) 
   if ( adopted[i] == 1 ) pred[i] = weibull_rng( k,lambda[i] );
   for ( i in 1:N ) 
   if ( adopted[i] == 0 ) pred[i] = weibull_lccdf(days_to_event[i] | k,lambda[i]);
}
5

did you remove the rounding? Also you need a lower = 0 on a

# not this
# black<-round(rweibull(1000,1.1,65))
# other<-round(rweibull(1000,1.1,55))
# this
black <- rweibull(1000, 1.1, 65)
other <- rweibull(1000, 1.1, 55)

I was able to run it in cmdstanr

6

Simulation:

black<-rweibull(1000,1.1,65)
other<-rweibull(1000,1.1,55)
d<-data.frame(days_to_event=c(black,other),color_id=c(rep(1,1000),rep(2,1000)))
d<- d %>%  mutate(adopted = if_else(days_to_event >= 200, 0, 1))

Stan model:

data{
    int N;
    int adopted[N];
    vector[N] days_to_event;
    int color_id[N];

}
parameters{
    vector<lower=0>[2] a;
    real<lower=0> k;
    real<lower=0> beta;
}
transformed parameters{
vector[N] lambda = exp( (a[color_id] * beta) * k);
}

model{
   a~exponential(0.4);
   beta~exponential(0.4);
   k~exponential(1);
   for ( i in 1:N ) {
        if (adopted[i] == 0 ) target += weibull_lccdf( days_to_event[i] | k, lambda[i] );
        else days_to_event[i] ~ weibull( k, lambda[i] );
     }
}
generated quantities{
  real pred[N];
   for ( i in 1:N ) 
   if ( adopted[i] == 1 ) pred[i] = weibull_rng( k,lambda[i] );
   for ( i in 1:N ) 
   if ( adopted[i] == 0 ) pred[i] = weibull_lccdf(days_to_event[i] | k,lambda[i]);
}

With vector<lower=0>[2] a; I get the same crash. Does it run on R/RStudio…? Maybe the model is so poorly specified that it crashes R but not cmdstanr?

7

Weird. I’m in rstudio and this runs

library(cmdstanr)
library(data.table)
mod <- cmdstan_model("weibull_crash.stan")

black <- rweibull(1000, 1.1, 65)
other <- rweibull(1000, 1.1, 55)

d <-data.table(days_to_event=c(black,other),color_id=c(rep(1,1000),rep(2,1000)))
d[, adopted := ifelse(days_to_event >= 200, 0, 1)]

dat <- list(
  N=nrow(d),
  days_to_event = d$days_to_event,
  color_id = d$color_id ,
  adopted = d$adopted
)
out <- mod$sample(data = dat,
                  parallel_chains = 4)
1 Like
8

RStan crashes are commonly caused by a couple of things:

  • Having -march=native in your Makevars.win
  • Recently upgrading from R3.x to R4.x and having the old .rds model files get loaded
  • Not having enough RAM to store all parameter estimates for the model

If none of those are your issue, then restart R (making sure rstan doesn’t get loaded) and try reinstalling rstan and StanHeaders from source via:

# Compile packages using all cores
Sys.setenv(MAKEFLAGS = paste0("-j",parallel::detectCores()))

install.packages(c("StanHeaders","rstan"),type="source")
3 Likes
9

Thanks @spinkney @andrjohns , reinstalling rstan solved the problem.

3 Likes
10

I observed the same problem agains shortly after. Removing the .rds file solved it.