Initialization failed in PyStan

Modeling
1

I’m trying to fit a hierarchical model for data distributed lognormally, with a normal prior for mu and uniform prior for sigma. I’m getting an initialization failed error without much other helpful output - any help would be greatly appreciated. Please note that my real sample size is much bigger :)

Apologies if this is a dumb question - I’m fairly new to this.

test_data = [22.51,
  48.34,
  19.43,
  28.81,
  23.69,
  9.76,
]

model_string = """
data {
    int<lower=0> control_n;
    vector<lower=0>[control_n] pA;
}

parameters {
  real muA;
  real sigmaA;
}

model {  
  muA ~ normal(3, 4);
  sigmaA ~ uniform(3, 5);
  
  pA ~ lognormal(muA, sigmaA);
}
"""

stan_data = {"pA": test_data, "control_n": len(test_data)}
sm = pystan.StanModel(model_code=model_string)
fit = sm.sampling(data=stan_data, iter=1000, chains=1)

The output is:

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_d24d16f55b52b2947d279a1230e8d96c NOW.
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-55-fee1dfe5f54b> in <module>
     28 stan_data = {"pA": test_data, "control_n": len(test_data)}
     29 sm = pystan.StanModel(model_code=model_string)
---> 30 fit = sm.sampling(data=stan_data, iter=1000, chains=1)

/Users/Shared/anaconda3/envs/stats/lib/python3.7/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)
    776         call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars))
    777         call_sampler_star = self.module._call_sampler_star
--> 778         ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
    779         samples = [smpl for _, smpl in ret_and_samples]
    780 

/Users/Shared/anaconda3/envs/stats/lib/python3.7/site-packages/pystan/model.py in _map_parallel(function, args, n_jobs)
     88             pool.join()
     89     else:
---> 90         map_result = list(map(function, args))
     91     return map_result
     92 

stanfit4anon_model_d24d16f55b52b2947d279a1230e8d96c_4922423782970638396.pyx in stanfit4anon_model_d24d16f55b52b2947d279a1230e8d96c_4922423782970638396._call_sampler_star()

stanfit4anon_model_d24d16f55b52b2947d279a1230e8d96c_4922423782970638396.pyx in stanfit4anon_model_d24d16f55b52b2947d279a1230e8d96c_4922423782970638396._call_sampler()

RuntimeError: Initialization failed.
2

Hi,

You could try with smaller init area

.sampling(..., init=0.5)

But I would try to use a bit more relaxed prior for sigma. Maybe try half-normal or half-student t (with some degree of freedom, e.g. 3)

sigma ~ normal(0,1);

Or

sigma ~ student_t(3, 0, 1);

I’m not totally sure what would the data look like with that kind of sigma for log-normal, maybe you could try to do some prior predictive checks and see.

3

Uniform priors are not recommended but if you do use them then you must also declare the bounds in parameter block

parameters {
  ...
  real<lower=3,upper=5> sigmaA;
}
model {
  sigmaA ~ uniform(3, 5);
  ...
}
4

Thank you! This solved the problem :)