Initialization failed in PyStan during covid-19 effective reproduction estimation model

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.

I set n_jobs = 1, but still got this error…


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-335-ec35bdf6f457> in <module>()
      5                 n_jobs=1,
      6                 pars = ['Rt', 's_smooth', 'domestic_infect', 'imported_infect', 'cases_infect'],
----> 7                 seed = 1234)

1 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)
     88             pool.join()
     89     else:
---> 90         map_result = list(map(function, args))
     91     return map_result
     92 

stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985.pyx in stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985._call_sampler_star()

stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985.pyx in stanfit4anon_model_41af6b195d0a37a4a461f78cd03c2699_5376968664909052985._call_sampler()

RuntimeError: Initialization failed.

Hi,
your questions seems to have slipped through, sorry for that.

I am not certain how PyStan error messages map to Stan error messages, but initialization failures usually mean that the log density cannot be computed (e.g. it ends up infinite or NaN) for the randomly chosen initial values of the parameters. It is however weird you don’t get a more detailed message maybe @ahartikainen can help with the PyStan part?

I however don’t see an immediate culprit. Maybe some of the poisson or continuous_poisson parameters (results of the convolution calls) end up negative? Note that you can add print statements to your program so that you can check both that all parts of the model are executed and inspect intermediate values

The model generally looks complex and a bit weird to me. I would suggest you start with a very simple model, fit it to simulated data and then add other layers in small increments, otherwise it is almost impossible to debug any issues.

Also you don’t have any prior on s_smooth and Rt which is quite likely to lead to problems down the line.

Best of luck with your model!

Thanks for your reply, actually I checked my code carefully and found that p_blah_blah_rev were wrong. I changed them and everything works perfectly. Here is my code, so that you guys run them in google colab.
This code is for open science, so you guys improve it to estimate Rt more precisely and sorry for some Japanese tho.