Hi,
I recently learned about Stan’s map_rect and I’m seriously excited about it. Since I work a lot with ODE models, I think I might just be the target group for this feature.
Below, there’s an example of how I would like to use an ODE integrator inside a function (unit_logLik) passed to map_rect. The data argument x_r contains the data, observation times, initial value, initial time and the dose data used in the ODE.
When I wanted to pass the data to integrate_ode_bdf I came upon slightly awkward behaviour. Since the initial value and the parameter I wanted to pass on to my ODE right hand side are only single real values, I first tried passing in
integrate_ode_bdf(pk_rhs, x_r[17], x_r[18], x_r[9:16], theta, x_r[19], x_i)
That however, doesn’t work of course since x_r[17] and x_r[19] evaluate to doubles and not to vectors, as expected. My next try was
integrate_ode_bdf(pk_rhs, { x_r[17] }, x_r[18], x_r[9:16], theta, { x_r[19] }, x_i)
which works fine for stanc but failes during compilation. As far as my limited knowledge of C++ could take me, the problem is that the Stan code { x_r[19] } gets compiled to
stan::math::array_builder<local_scalar_t__>(...)
where local_scalar_t__ is stan::math::var, therefore returning a std::vector<stan::math::var> instead of a std::vector<double> as expected by integrate_ode_bdf. This wasn’t a problem for the initial value, since those can be parameters, I assume. The only solution that worked in the end was to replace { x_r[17] } and { x_r[19] } by x_r[17:17] and x_r[19:19].
It’s not a big deal to do this, but it feels slightly awkward that the first way doesn’t work (even though I understand from the C++ code why it cannot). I’m not sure if stanc could be capable of capturing a mistake like this and warning the user about it? Or I guess it should at least be found in the manual. Which it maybe is? I didn’t find anything…
functions {
real[] pk_rhs(real t,
real[] y,
real[] theta,
real[] x_r,
int[] x_i) {
real ka = theta[1];
real f = theta[2];
real Vmax = theta[3];
real Km = theta[4];
real V = theta[5];
real dose = x_r[1];
return { (ka * f * dose * exp(-ka * (t + 2))
- Vmax * y[1] / (Km + y[1])) / V };
}
vector unit_logLik(vector shared,
vector unit_spec,
real[] x_r,
int[] x_i) {
real theta[5] = to_array_1d(shared[1:5]);
real sigma = shared[6];
int n_obs = x_i[1];
real y_hat[n_obs, 1] =
integrate_ode_bdf(pk_rhs, x_r[17:17],
x_r[18], x_r[9:16],
theta, x_r[19:19], x_i);
vector[n_obs] lp;
for (i in 1:n_obs) {
lp[i] = normal_lpdf(x_r[i] | log(y_hat[i,1]), sigma);
}
return lp;
}
}