Need help. Keep getting initialization error

Modeling
1

Hi everyone,

I am new to Stan (with Pystan for the interface). I am currently trying to solve a simple inference problem involving differential equation for modeling bacterial growth rate, which code as attached below. I keep getting “RuntimeError: Initialization failed.” error.

functions {
        real[] growth(real t, real[] y, real[] p, real[] x_r, int[] x_i) {
            real dydt[1];
            dydt[1] = (p[1] * (1 - (y[1]/p[2]))) * y[1];
            return dydt;
        }
    }
    data {
        int<lower=1> T;
        real y[T, 1];
        real t0;
        real ts[T];
    }
    transformed data {
        real x_r[0];
        int x_i[0];
    }
    parameters {
        real<lower=0> sigma;
        real y0[1]; //infer the initial state
        real p[2]; //model parameters
    }
    model {
        real y_hat[T, 1];
        sigma ~ uniform(0.1, 0.5);
        p[1] ~ uniform(0, 1);
        p[2] ~ uniform(0, 2);
        y0[1] ~ uniform(0, 0.2);
        
        y_hat = integrate_ode_rk45(growth, y0, t0, ts, p, x_r, x_i);
        for (t in 1:T)
            y[t] ~ normal(y_hat[t], sigma);
    }
}

What did I do wrong here? Here is the input dictionary in Python:

data = {
    'T': len(od),
    'ts': od.index, #timeframe as the index for the data
    'y': od, #data to be fitted
    't0': 0
}

Also, when I was trying to find a solution for ODE inference from the internet, I found some people put the integrate_ode_rk45 function in the transformed parameters block, and some others put it in the model block. Which one is the correct way to do this?

Thanks in advance for the help and suggestions!

Here is the full error message, if it helps:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-54-7445a20b3355> in <module>
     51 
     52 # Train the model and generate samples
---> 53 fit = sm.sampling(data=data, iter=100, chains=2, n_jobs=1)
     54 print(fit)

~/Documents/Notebook/env/lib/python3.9/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 

~/Documents/Notebook/env/lib/python3.9/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_4b7224b9f462e46cebdf72f52a562750_685902710779344041.pyx in stanfit4anon_model_4b7224b9f462e46cebdf72f52a562750_685902710779344041._call_sampler_star()

stanfit4anon_model_4b7224b9f462e46cebdf72f52a562750_685902710779344041.pyx in stanfit4anon_model_4b7224b9f462e46cebdf72f52a562750_685902710779344041._call_sampler()

RuntimeError: Initialization failed.
1 Like
2

Hm, not sure! try adding print statements to discern which variables are taking unexpected values leading to a NaN or -/+Inf for the target.

The only difference between the TP & model blocks when it comes to declaring/computing variables is that variables in the TP block:

  1. Will be saved to file in the output (model-block-declared variables will not be saved)
  2. Are permitted to have bounds in their declaration, in which case the sampler will throw an error when their computation falls outside the declared bounds. This may not halt the sampler, but will do so if too many samples lead to a TP falling outside its bounds.
1 Like
3

I haven’t looked closely, but biologically it think you need y0 to be positive, and maybe the growth equation sails off to negative infinity if y0 is negative? Try declaring y0 with the appropriate constraint and see what happens.