Silently ignored error when specifying init values


I’ve been trying to fix initialisation errors in a wiener family model, where there are multiple constraints (such as the non-decision time being lower than the smallest observed outcome).

However, I noticed that my initialisation values were seemingly not passed to Stan. I’ve identified this as resulting from a silent error when passing initial values for real Stan parameters as lists, rather than single values. The passed values are simply ignored. Making this mistake the other way around (specifying single inits for vector parameters, even if of length 1) results in a dimension error.

To try and replicate this I set up a toy problem:

b <- c(2, 0.75)
x <- rnorm(200)
group <- as.integer(rbernoulli(200))
y <- rnorm(200, mean = 1 + group + x)
dat1 <- data.frame(x, y, group)

formula <- bf(y ~ 1 + group, sigma ~ group, family=brmsfamily("gaussian", link_sigma="identity"))

write(make_stancode(formula, data = dat1), file = "test-stancode.stan")
standat <- make_standata(formula, dat1)

This produces the following Stan code:

// generated with brms 2.13.5
functions {
data {
  int<lower=1> N;  // number of observations
  vector[N] Y;  // response variable
  int<lower=1> K;  // number of population-level effects
  matrix[N, K] X;  // population-level design matrix
  int<lower=1> K_sigma;  // number of population-level effects
  matrix[N, K_sigma] X_sigma;  // population-level design matrix
  int prior_only;  // should the likelihood be ignored?
transformed data {
  int Kc = K - 1;
  matrix[N, Kc] Xc;  // centered version of X without an intercept
  vector[Kc] means_X;  // column means of X before centering
  int Kc_sigma = K_sigma - 1;
  matrix[N, Kc_sigma] Xc_sigma;  // centered version of X_sigma without an intercept
  vector[Kc_sigma] means_X_sigma;  // column means of X_sigma before centering
  for (i in 2:K) {
    means_X[i - 1] = mean(X[, i]);
    Xc[, i - 1] = X[, i] - means_X[i - 1];
  for (i in 2:K_sigma) {
    means_X_sigma[i - 1] = mean(X_sigma[, i]);
    Xc_sigma[, i - 1] = X_sigma[, i] - means_X_sigma[i - 1];
parameters {
  vector[Kc] b;  // population-level effects
  real Intercept;  // temporary intercept for centered predictors
  vector[Kc_sigma] b_sigma;  // population-level effects
  real Intercept_sigma;  // temporary intercept for centered predictors
transformed parameters {
model {
  // initialize linear predictor term
  vector[N] mu = Intercept + Xc * b;
  // initialize linear predictor term
  vector[N] sigma = Intercept_sigma + Xc_sigma * b_sigma;
  // priors including all constants
  target += student_t_lpdf(Intercept | 3, 1.4, 2.5);
  target += student_t_lpdf(Intercept_sigma | 3, 0, 2.5);
  // likelihood including all constants
  if (!prior_only) {
    target += normal_lpdf(Y | mu, sigma);
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept - dot_product(means_X, b);
  // actual population-level intercept
  real b_sigma_Intercept = Intercept_sigma - dot_product(means_X_sigma, b_sigma);

If I now try to run this while providing initial values that I know will violate constraints, I noticed that initial values for real parameters are ignored if they are lists:

init_fun <- function() {
    Intercept = 0,
    b = list(0),
    Intercept_sigma = list(-20),
    b_sigma = list(0)
rstan::stan("test-stancode.stan", init = init_fun, data = standat)

Sampling will proceed as normal, even though I would have expected it to either:

  • throw an error stating that real parameters can not be passed lists as initial values
  • somehow use the value provided, which would lead sampling to fail

This is in contrast to when supplying values that should be a list but are passed as single values instead:

init_fun <- function() {
    Intercept = 0,
    b = 0,
    Intercept_sigma = list(-20),
    b_sigma = list(0)
rstan::stan("test-stancode.stan", init = init_fun, data = standat)

Which produces the output:

Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b; dims declared=(1); dims found=()
[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                           
[2] "  mismatch in number dimensions declared and found in context; processing stage=parameter initialization; variable name=b; dims declared=(1); dims found=()"

It would be great if a similar error could be thrown for the first example I posted, as silently ignoring mismatching value dimensions makes it very hard to debug a program.

PS: For reference, I initially incorrectly raised this as an issue at the brms repository:

Sorry - I just noticed that list values for vector parameters are also ignored 🙃 they should be atomic vectors or matrices, I guess?

While this is probably because of my use of lists in a context where a more experienced R user wouldn’t use them, I wonder if it would be possible to raise an error anyways?

For inits you need a list of list, where the first “list” needs as many elements as you have chains.
The default number of chains in 4, so that your stan call with inits should work with the following:

            init = lapply(1:4,init_fun()),
            data = standat)

Thanks - passing the inits down to Stan is actually working using the code I posted (passing init_fun as the init argument).

My problem is with how RStan deals with init values inside that list that aren’t the right data type. If for example a value inside the init list is a list itself, this gets silently ignored.

OK thanks for the pointer ;-).

I think I also understand where the misunderstanding is.
The documentation sais that inits can be a list of lists.
What this means is that if you have 4 chains, you need a list with 4 entries, each of which is a list.

my_intits = list(
              list(Intercept = 0, b = rep(0, K-1),...),
              list(Intercept = 0, b = rep(0, K-1),...),
              list(Intercept = 0, b = rep(0, K-1),...),
              list(Intercept = 0, b = rep(0, K-1),...))

But each list in my_intits above must not contain lists, only reals, integers, vectors, matrices and arrays.
As you can see from my example for paramter b above, you need to submit the init values with the correct dimensions.

Thanks! That makes sense (although I would say it’s not 100% clear from the documentation).

Is it possible to raise a warning if values in the named list are not valid init values? At the moment, this is silently ignored. Some sort of assertion on the content of the init list would be very useful.

I was trying to figure out where the best place to do this might be - maybe somewhere here?

I agree that instructions could be a bit clearer.
I’m not involved in the development of brms/rstan/cmndstanr, but I guess that if you make a new issue (in rstan and cmdstanr, the brms backends) and even make a small pull request, this will be welcomed.
(I think the most useful thing would be examples in the documentation.)

Also, good luck with the Wiener. A key problem that cropped up from time to time was that one needs to be careful that model/prior for the NDT do not allow numbers larger than minimum observed RT.

Thanks! I will open those issues and give my examples, that’s a great idea.

Yes, the NDT is giving me a real headache :) hopefully I’ll get there eventually.

Thanks @Guido_Biele for helping out here!

@janfreyberg You’re right that the warnings (or lack thereof) and doc could be better here, so yeah definitely open those issues and we’ll try to make improvements. Thank you!