I’m trying to make a model about COVID-19 in Pystan on published data. But, I got the initialization failed error. I have checked all the previous questions and the data format for me is changed to np.array.astype(int).
Error
---------------------------------------------------------------------------
RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
File "/usr/lib/python3.6/multiprocessing/pool.py", line 119, in worker
result = (True, func(*args, **kwds))
File "/usr/lib/python3.6/multiprocessing/pool.py", line 44, in mapstar
return list(map(*args))
File "stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985.pyx", line 373, in stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985._call_sampler_star
File "stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985.pyx", line 406, in stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985._call_sampler
RuntimeError: Initialization failed.
"""
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
<ipython-input-332-5ad30f3913a4> in <module>()
4 thin = 5,
5 pars = ['Rt', 's_smooth', 'domestic_infect', 'imported_infect', 'cases_infect'],
----> 6 seed = 1234)
3 frames
/usr/local/lib/python3.6/dist-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
/usr/local/lib/python3.6/dist-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()
/usr/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):
/usr/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.
Here is the block of stan model,
functions {
// calculating the convolutions
vector convolution(vector x, vector y_rev) {
int T = num_elements(x);
vector[T-1] res;
for (t in 2:T) {
res[t-1] = dot_product(x[1:(t-1)], y_rev[(T-t+2):T]);
}
return res;
}
// continuous Poisson distribution
real continuous_poisson_lpdf(vector y, vector lam) {
real lp = sum(y .* log(lam) - lam - lgamma(y + 1));
return lp;
}
}
data {
int<lower=1> T; // number of days
int<lower=0> domestic_onset[T];
int<lower=0> imported_onset[T];
int<lower=0> domestic_report[T];
int<lower=0> imported_report[T];
vector[T] p_gt_rev; // generation time distribution (reverse order)
vector[T] p_inc_rev; // incubation period distribution (reverse order)
vector[T] p_rep_rev; // report time distribution (reverse order)
}
parameters {
vector<lower=0>[T-1] Rt;
real<lower=0> s_smooth; // smoothing factor
simplex[T] domestic_infect_from_onset_raw; // distribution of domestic infected people from onset data
simplex[T] imported_infect_from_onset_raw;
simplex[T] domestic_infect_from_report_raw;
simplex[T] imported_infect_from_report_raw;
}
transformed parameters {
vector[T] domestic_infect_from_onset = domestic_infect_from_onset_raw * sum(domestic_onset); // Number of domestic infected people by using back projection from onset data.
vector[T] imported_infect_from_onset = imported_infect_from_onset_raw * sum(imported_onset);
vector[T] domestic_infect_from_report = domestic_infect_from_report_raw * sum(domestic_report);
vector[T] imported_infect_from_report = imported_infect_from_report_raw * sum(imported_report);
vector[T] domestic_infect = domestic_infect_from_onset + domestic_infect_from_report;
vector[T] imported_infect = imported_infect_from_onset + imported_infect_from_report;
vector[T] cases_infect = domestic_infect + imported_infect;
}
model {
// smoothing on the raw scale
domestic_infect_from_onset_raw[2:T] ~ normal(domestic_infect_from_onset_raw[1:(T-1)], s_smooth);
imported_infect_from_onset_raw[2:T] ~ normal(imported_infect_from_onset_raw[1:(T-1)], s_smooth);
domestic_infect_from_report_raw[2:T] ~ normal(domestic_infect_from_report_raw[1:(T-1)], s_smooth);
imported_infect_from_report_raw[2:T] ~ normal(imported_infect_from_report_raw[1:(T-1)], s_smooth);
// back projection
domestic_onset[2:T] ~ poisson(convolution(domestic_infect_from_onset, p_inc_rev));
imported_onset[2:T] ~ poisson(convolution(imported_infect_from_onset, p_inc_rev));
domestic_report[2:T] ~ poisson(convolution(domestic_infect_from_report, p_rep_rev));
imported_report[2:T] ~ poisson(convolution(imported_infect_from_report, p_rep_rev));
// estimating effective reproduction number
domestic_infect[2:T] ~ continuous_poisson(Rt .* convolution(cases_infect, p_gt_rev));
}
And data is like this,
# Define distributions
epsilon = 1e-5
def weib(x,n,a):
return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)
def lnorm(x,m,s):
return 1 / (np.sqrt(2 * math.pi) * s * (x + epsilon)) * np.exp(- (np.log(np.maximum(x-m,0) + epsilon))**2 / (2 * s**2))
# Discretesized distribution
##------ predetermined distributions ------
## discretesized distribution of generation time (Nishiura, et al, 2020)
T = len(df_cases)
c_gt = weib(np.array(range(0,T)), 2.305, 5.452)
p_gt_rev = list(np.sort(c_gt)[::-1])
## discretesized distribution of incubation period (from infection to onset) (Linton et al 2020 - with right truncation excl. Wuhan residents)
c_inc = lnorm(np.array(range(0,T)), 1.519, 0.615)
p_inc_rev = list(np.sort(c_inc)[::-1])
## discretesized distribution of the delay from infection to report
# right-truncated reporting delay from onset to repot (estimated values using the MHLW data): dweibull(t, shape = 1.741, scale = 8.573)
def val(t):
def integrand(x):
return weib(t-x,1.741,8.573) * lnorm(x,1.519,0.615)
val, _= integrate.quad(integrand,0,t)
return val
c_rep = [val(t) for t in range(T)]
p_rep_rev = list(np.sort(c_rep)[::-1])
where p_gt_rev, p_inc_rev and p_rep_rev are 113 length-list type data without ‘NaN’. And others are like below;
do_onset.astype(int).values
array([ 0, 1, 3, 0, 1, 1, 1, 1, 0, 20, 2, 0, 2,
1, 1, 1, 5, 9, 3, 5, 1, 6, 4, 10, 11, 13,
1, 1, 6, 7, 16, 7, 15, 18, 19, 24, 34, 35, 28,
22, 36, 29, 36, 26, 27, 33, 33, 39, 30, 38, 41, 40,
28, 32, 34, 44, 59, 82, 105, 178, 84, 159, 206, 324, 28,
90, 170, 474, 28, 311, 539, 321, 25, 535, 446, 407, 414, 383,
472, 411, 385, 336, 371, 356, 324, 293, 283, 280, 226, 225, 216,
157, 219, 155, 162, 115, 148, 106, 105, 120, 94, 69, 79, 72,
61, 39, 30, 18, 19, 8, 5, 0, 0])
im_onset.astype(int).values
array([ 1, 1, 0, 1, 1, 1, 1, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0,
2, 1, 1, 0, 0, 0, 0, 2, 1, 0, 0, 1, 0, 4, 0, 0, 1,
1, 3, 0, 1, 6, 1, 0, 0, 5, 2, 5, 2, 7, 4, 5, 16, 9,
2, 7, 3, 6, 5, 3, 13, 24, 11, 18, 13, 12, 16, 13, 11, 12, 11,
8, 6, 8, 1, 5, 5, 1, 3, 0, 0, 2, 2, 2, 0, 2, 0, 3,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
do_reported.astype(int).values
array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 7, 2, 1, 1, 0, 1, 1, 1, 2, 0, 2, 0,
0, 1, 0, 0, 4, 0, 2, 4, 2, 2, 1, 4, 6,
0, 4, 0, 0, 4, 2, 0, 11, 1, 5, 0, 4, 5,
8, 1, 20, 55, 20, 61, 24, 60, 72, 56, 26, 21, 69,
147, 90, 171, 71, 28, 31, 40, 69, 83, 115, 111, 124, 59,
89, 120, 58, 76, 87, 50, 22, 40, 105, 63, 123, 55, 37,
48, 43, 33, 24, 17, 10, 15, 16, 4])
im_reported.astype(int).values
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1,
6, 1, 1, 1, 1, 0, 1, 2, 2, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0])
All length are T = 113.
In summary, i couldn’t find any error in data, so maybe i guess problems are in bounds setting or any other blocks. I’d really appreciate it if you some guys helped me.