Initialization failed

I am trying to fit the hierarchical bayesian logistic model using pystan:
#construct the data block
bmtl_code = “”"
data {
int<lower=0> N; //num observations
int<lower=1> J; // num predicators
int<lower=1> K; //num tasks
int<lower=1, upper=K> kk[N]; //task label for each observation
matrix[N,J] x; //observation predicators
int y[N]; //observation outcomes
matrix[J,K] Zero;
}

parameters {
vector[K] alpha; //intercept alpha
matrix[J, K] beta; //coefficients beta
corr_matrix[K] Omega; //prior correlation
vector<lower=0>[K] sigma; //prior scale
real<lower=0> tau;
vector<lower=0>[J] psi;
}

transformed parameters {
vector[J] r; //shrinkage scalar
vector[J] r_square;
r = tau*psi;
for (j in 1:J)
r_square[j]=r[j]*r[j];
}

model {
psi ~ cauchy(0,1);
tau ~ cauchy(0,1);

alpha ~ cauchy(0,10);
sigma ~ cauchy(0,2.5);
Omega ~ lkj_corr(1);
for (j in 1:J)
beta[j] ~ multi_normal(Zero[j], r_square[j] * quad_form_diag(Omega, sigma));

for (n in 1:N)
y[n] ~ bernoulli(inv_logit(alpha[kk[n]] + dot_product(x[n],beta[,kk[n]])));
}
“”"
bmtl_data = dict(N=X_train.shape[0], J=num_features, K=num_tasks,
kk=df_train[‘label’].values, x=X_train, y=y_train,
Zero=np.zeros((num_features,num_tasks)))

fit = pystan.stan(model_code=bmtl_code, data=bmtl_data)

I ran the code above and it returned the following message:

RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
“”"
Traceback (most recent call last):
File “/anaconda3/lib/python3.6/multiprocessing/pool.py”, line 119, in worker
result = (True, func(*args, **kwds))
File “/anaconda3/lib/python3.6/multiprocessing/pool.py”, line 44, in mapstar
return list(map(*args))
File “stanfit4anon_model_6f6086b057a6e799617a8be0e0bc8dce_4816030848099204205.pyx”, line 368, in stanfit4anon_model_6f6086b057a6e799617a8be0e0bc8dce_4816030848099204205._call_sampler_star
File “stanfit4anon_model_6f6086b057a6e799617a8be0e0bc8dce_4816030848099204205.pyx”, line 401, in stanfit4anon_model_6f6086b057a6e799617a8be0e0bc8dce_4816030848099204205._call_sampler
RuntimeError: Initialization failed.
“”"

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

RuntimeError Traceback (most recent call last)
in ()
3 Zero=np.zeros((num_features,num_tasks)))
4
----> 5 fit = pystan.stan(model_code=bmtl_code, data=bmtl_data)

/anaconda3/lib/python3.6/site-packages/pystan/api.py in stan(file, model_name, model_code, fit, data, pars, chains, iter, warmup, thin, init, seed, algorithm, control, sample_file, diagnostic_file, verbose, boost_lib, eigen_lib, n_jobs, **kwargs)
400 sample_file=sample_file, diagnostic_file=diagnostic_file,
401 verbose=verbose, algorithm=algorithm, control=control,
–> 402 n_jobs=n_jobs, **kwargs)
403 return fit

/anaconda3/lib/python3.6/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)
724 call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars))
725 call_sampler_star = self.module._call_sampler_star
–> 726 ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
727 samples = [smpl for _, smpl in ret_and_samples]
728

/anaconda3/lib/python3.6/site-packages/pystan/model.py in _map_parallel(function, args, n_jobs)
79 try:
80 pool = multiprocessing.Pool(processes=n_jobs)
—> 81 map_result = pool.map(function, args)
82 finally:
83 pool.close()

/anaconda3/lib/python3.6/multiprocessing/pool.py in map(self, func, iterable, chunksize)
264 in a list that is returned.
265 ‘’’
–> 266 return self._map_async(func, iterable, mapstar, chunksize).get()
267
268 def starmap(self, func, iterable, chunksize=None):

/anaconda3/lib/python3.6/multiprocessing/pool.py in get(self, timeout)
642 return self._value
643 else:
–> 644 raise self._value
645
646 def _set(self, i, obj):

RuntimeError: Initialization failed.

I am new to pystan and really not sure how to debug.

Any help would be greatly appreciated!
Jasmine

Use bernoulli_logit(eta) rather than bernoulli(inv_logit(eta)). Also, you would be much better served to fill in a vector[N] eta like

for (n in 1:N) {
  eta[n] = alpha[kk[n]] + dot_product(x[n],beta[,kk[n]];
}

and then evaluate the log-likelihood once with

y ~ bernoulli_logit(eta);

Thanks for the quick reply.
I have modified my code as you suggest and it seems initializing successfully. However, the code has run for 3 hours without any error or output. There are 64612 observations and 24 features in my model. Without any prompt information, I don’t know how far does my code run. I am wondering if something wrong in my model.

My code:
#construct the data block
bmtl_code = “”"
data {
int<lower=0> N; //num observations
int<lower=1> J; // num predicators
int<lower=1> K; //num tasks
int<lower=1, upper=K> kk[N]; //task label for each observation
matrix[N,J] x; //observation predicators
int y[N]; //observation outcomes
matrix[J,K] Zero;
}

parameters {
vector[K] alpha; //intercept alpha
matrix[J, K] beta; //coefficients beta
corr_matrix[K] Omega; //prior correlation
vector<lower=0>[K] sigma; //prior scale
real<lower=0> tau;
vector<lower=0>[J] psi;
}

transformed parameters {
vector[J] r; //shrinkage scalar
vector[J] r_square;
vector[N] eta;
r = tau*psi;
for (j in 1:J){
r_square[j]=r[j]*r[j];
}
for (n in 1:N){
eta[n]=alpha[kk[n]] + dot_product(x[n],beta[,kk[n]]);
}
}

model {
psi ~ cauchy(0,1);
tau ~ cauchy(0,1);

alpha ~ cauchy(0,10);
sigma ~ cauchy(0,2.5);
Omega ~ lkj_corr(1);
for (j in 1:J){
beta[j] ~ multi_normal(Zero[j], r_square[j] * quad_form_diag(Omega, sigma));
}
y ~ bernoulli_logit(eta);
}
“”"
bmtl_data = dict(N=X_train.shape[0], J=num_features, K=num_tasks,
kk=df_train[‘label’].values, x=X_train, y=y_train,
Zero=np.zeros((num_features,num_tasks)))

fit = pystan.stan(model_code=bmtl_code, data=bmtl_data)

Running message:
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_4a5092a673047c964a13bcadc3583b74 NOW.
In file included from /var/folders/vn/sc2r7p_1557_lflrw2r6jg740000gn/T/tmpol8xa1rx/stanfit4anon_model_4a5092a673047c964a13bcadc3583b74_2693544027111664653.cpp:599:
In file included from /anaconda3/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4:
In file included from /anaconda3/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18:
In file included from /anaconda3/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1816:
/anaconda3/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: "Using deprecated NumPy API, disable it by " “#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION” [-W#warnings]
#warning "Using deprecated NumPy API, disable it by "
^
/var/folders/vn/sc2r7p_1557_lflrw2r6jg740000gn/T/tmpol8xa1rx/stanfit4anon_model_4a5092a673047c964a13bcadc3583b74_2693544027111664653.cpp:9155:30: warning: comparison of integers of different signs: ‘Py_ssize_t’ (aka ‘long’) and ‘std::__1::vector<std::__1::basic_string, std::__1::allocator<std::__1::basic_string > >::size_type’ (aka ‘unsigned long’) [-Wsign-compare]
__pyx_t_12 = ((__pyx_t_9 != __pyx_v_fitptr->param_names_oi().size()) != 0);
~~~~~~~~~ ^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /var/folders/vn/sc2r7p_1557_lflrw2r6jg740000gn/T/tmpol8xa1rx/stanfit4anon_model_4a5092a673047c964a13bcadc3583b74_2693544027111664653.cpp:603:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan_fit.hpp:22:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/src/stan/services/diagnose/diagnose.hpp:10:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/src/stan/model/test_gradients.hpp:7:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/src/stan/model/log_prob_grad.hpp:4:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/stan/math/rev/mat.hpp:12:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/stan/math/prim/mat.hpp:298:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/stan/math/prim/arr.hpp:38:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/stan/math/prim/arr/functor/integrate_ode_rk45.hpp:17:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/numeric/odeint.hpp:61:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/numeric/odeint/util/multi_array_adaption.hpp:29:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array.hpp:21:
In file included from /anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array/base.hpp:28:
/anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array/concept_checks.hpp:42:43: warning: unused typedef ‘index_range’ [-Wunused-local-typedef]
typedef typename Array::index_range index_range;
^
/anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array/concept_checks.hpp:43:37: warning: unused typedef ‘index’ [-Wunused-local-typedef]
typedef typename Array::index index;
^
/anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array/concept_checks.hpp:53:43: warning: unused typedef ‘index_range’ [-Wunused-local-typedef]
typedef typename Array::index_range index_range;
^
/anaconda3/lib/python3.6/site-packages/pystan/stan/lib/stan_math/lib/boost_1.64.0/boost/multi_array/concept_checks.hpp:54:37: warning: unused typedef ‘index’ [-Wunused-local-typedef]
typedef typename Array::index index;
^
6 warnings generated.
/anaconda3/lib/python3.6/site-packages/pystan/misc.py:399: FutureWarning: Conversion of the second argument of issubdtype from float to np.floating is deprecated. In future, it will be treated as np.float64 == np.dtype(float).type.
elif np.issubdtype(np.asarray(v).dtype, float):

As you’ve figured out, there’s probably some funky with your model or data :P. The way to figure out what is going on is to start with small amounts of simulated data.

Simulate from your model and see if you recover the parameters you used to generate things. Start with really small bits of data – enough so that you can run your model in less than a minute. Waiting three hours to realize there is a bug is pretty slow :P. You’ll have to take your time on this.

Non-identified parameters or giant correlations in your output can cause these slowdowns. Pairplots (like these: seaborn.pairplot — seaborn 0.13.0 documentation) can be really helping for finding problems.