PyStan RuntimeError: Initialization Failed and slow compilation of STAN

I have the same error message as mentioned for example here (PyStan: RuntimeError: Initialization Failed - #5 by patrickjmc or Initialization failed in Pystan), I tried to search, however, I have not found any solution which would work for me. Also, the compilation of my model is very slow.


My stan file is the following:

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1, upper=T> Tsubj[N];
  real<lower=0> gain[N, T];
  real<lower=0> loss[N, T];  // absolute loss amount
  real cert[N, T];
  int<lower=-1, upper=1> gamble[N, T];
}

transformed data {
}

parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; // how many parameters we have
  vector[N] betag_pr;
  vector[N] betal_pr;
  vector[N] tau_pr;
}

transformed parameters {
  vector<lower=0, upper=20>[N] betag;
  vector<lower=0, upper=20>[N] betal;
  vector<lower=0, upper=30>[N] tau;

  for (i in 1:N) {
    betag[i] = Phi_approx(mu_pr[1] + sigma[1] * betag_pr[i]) * 20;
    betal[i] = Phi_approx(mu_pr[2] + sigma[2] * betal_pr[i]) * 20;
    tau[i]    = Phi_approx(mu_pr[3] + sigma[3] * tau_pr[i]) * 30;
  }
}

model {
  // linear fit based on Tom et all. without lambda (lambda is betag/betal)
  // hyper parameters
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);

  // individual parameters w/ Matt trick
  betag_pr ~ normal(0, 1.0);
  betal_pr ~ normal(0, 1.0);
  tau_pr    ~ normal(0, 1.0);

  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
      real evGamble;  // they are left as arrays as an example for RL models.
      real pGamble;

      // loss[i, t]=absolute amount of loss (pre-converted in R)
      evSafe   = cert[i, t];
      evGamble = 0.5 * (betag[i] * gain[i, t] - betal[i] * loss[i, t]);
      pGamble  = inv_logit(tau[i] * (evGamble - evSafe));
      gamble[i, t] ~ bernoulli(pGamble);
    }
  }
}
generated quantities {
  real<lower=0, upper=20> mu_betag;
  real<lower=0, upper=20> mu_betal;
  real<lower=0, upper=30> mu_tau;

  real log_lik[N];

  // For posterior predictive check
  real y_pred[N, T];

  // Set all posterior predictions to 0 (avoids NULL values)
  for (i in 1:N) {
    for (t in 1:T) {
      y_pred[i, t] = -1;
    }
  }

  mu_betag = Phi_approx(mu_pr[1]) * 20;
  mu_betal = Phi_approx(mu_pr[2]) * 20;
  mu_tau    = Phi_approx(mu_pr[3]) * 30;

  { // local section, this saves time and space
    for (i in 1:N) {
      log_lik[i] = 0;
      for (t in 1:Tsubj[i]) {
        real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
        real evGamble;  // they are left as arrays as an example for RL models.
        real pGamble;

        evSafe   = cert[i, t];
        evGamble = 0.5 * (betag[i] * gain[i, t] - betal[i] * loss[i, t]);
        pGamble    = inv_logit(tau[i] * (evGamble - evSafe));
        log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);

        // generate posterior prediction for current trial
        y_pred[i, t] = bernoulli_rng(pGamble);
      }
    }
  }
}

and I think that the problem is not in data (like in the mentioned post), since I have tested the same data both on ra_prospect and ra_noRA (hBayesDM Getting Started) and this code is a direct derivate of those two. The data are also of the same format specified/required by these.


The compilation of the code works fine:

import time
from datetime import timedelta
# Compile the model

start_compile = time.time()

stan_model_linearfit = pystan.StanModel(file=model_path)

end_compile = time.time()

print('Duration of compilation: ', timedelta(seconds=end_compile - start_compile))
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_7b9ec17863952415fa17979fa588c8d5 NOW.

Duration of compilation:  0:03:21.283095

However, also takes much longer than the original two STAN files mentioned above (any idea why?). The comparative is compilation of the no_RA.stan which took 0:00:48.568954 and the ra_prospect took 0:01:01.220933.

The code for running the model is then:

# run the model
start_model = time.time()
fit_linearfit = stan_model_linearfit.sampling(data=LA_data_prep, 
                    chains=4, iter=2000, warmup=1000, thin=1, init='random', verbose=True, 
                    control = {"adapt_delta":0.95, "stepsize":1, "max_treedepth":10}, n_jobs=-1)
end_model = time.time()
print('Duration of the model run is: ', timedelta(seconds=end_model - start_model))

and the error message is:

RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
“”"
Traceback (most recent call last):
File “C:\Users\user\AppData\Local\Continuum\anaconda3\lib\multiprocessing\pool.py”, line 121, in worker
result = (True, func(*args, **kwds))
File “C:\Users\user\AppData\Local\Continuum\anaconda3\lib\multiprocessing\pool.py”, line 44, in mapstar
return list(map(*args))
File “stanfit4anon_model_7b9ec17863952415fa17979fa588c8d5_2527243203695064452.pyx”, line 373, in stanfit4anon_model_7b9ec17863952415fa17979fa588c8d5_2527243203695064452._call_sampler_star
File “stanfit4anon_model_7b9ec17863952415fa17979fa588c8d5_2527243203695064452.pyx”, line 406, in stanfit4anon_model_7b9ec17863952415fa17979fa588c8d5_2527243203695064452._call_sampler
RuntimeError: Initialization failed.
“”"

The above exception was the direct cause of the following exception:

RuntimeError Traceback (most recent call last)
in
3 fit_linearfit = stan_model_linearfit.sampling(data=LA_data_prep,
4 chains=4, iter=2000, warmup=1000, thin=1, init=‘random’, verbose=True,
----> 5 control = {“adapt_delta”:0.95, “stepsize”:1, “max_treedepth”:10}, n_jobs=-1)
6 end_model = time.time()
7 print('Duration of the model run is: ', timedelta(seconds=end_model - start_model))

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pystan\model.py in sampling(self, data, pars, chains, iter, warmup, thin, seed, init, sample_file, diagnostic_file, verbose, algorithm, control, n_jobs, **kwargs)
811 call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars))
812 call_sampler_star = self.module._call_sampler_star
→ 813 ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
814 samples = [smpl for _, smpl in ret_and_samples]
815

~\AppData\Local\Continuum\anaconda3\lib\site-packages\pystan\model.py in _map_parallel(function, args, n_jobs)
83 try:
84 pool = multiprocessing.Pool(processes=n_jobs)
—> 85 map_result = pool.map(function, args)
86 finally:
87 pool.close()

~\AppData\Local\Continuum\anaconda3\lib\multiprocessing\pool.py in map(self, func, iterable, chunksize)
266 in a list that is returned.
267 ‘’’
→ 268 return self._map_async(func, iterable, mapstar, chunksize).get()
269
270 def starmap(self, func, iterable, chunksize=None):

~\AppData\Local\Continuum\anaconda3\lib\multiprocessing\pool.py in get(self, timeout)
655 return self._value
656 else:
→ 657 raise self._value
658
659 def _set(self, i, obj):

RuntimeError: Initialization failed.

Hi, maybe try with another set of initial values

Here are some examples:

init = '1'
init = '0.5'
init = '0'
init = {"mu_pr" : [0.1]*3, "sigma" : [0.1]*3, "betag_pr" : [1]*N, "betal": [1]*N, "tau" : [1]*N}.get

Thanks for the reply. I tried all those but none worked.

This may suffer from rounding errors. It’s more stable when written as

gamble[i, t] ~ bernoulli_logit(tau[i] * (evGamble - evSafe));
1 Like

It seems to be working! Thanks. Can you explain a bit more what the problem was and how this solves it?

Should the same also apply to here? I did not change this, I just wonder what the difference is.

pGamble = inv_logit(tau[i] * (evGamble - evSafe));
log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);