Define transformed parameters

[lamprosbouranis]

Hello,

Congrats on this very interesting work. Up to now I have written my own Stan program for Bayesian analysis of compartmental epidemic models, like the SEIR. A typical ODE system involves 18 states, two fixed parameters and 10 parameters to sample from.

Typically, I solve the ODE with a user-defined trapezoidal rule, which I have adapted to the needs of my work, and is faster than an application with RK45/BDF. The posterior outputs agree irrespectively of the solver used. Extending the time horizon of my analysis creates a computational bottleneck, with an HMC run with the trapezoidal rule requiring almost 40 hours with 210 daily time points and your work could offer a faster alternative with good accuracy of posterior estimates.

By construction, I use the transformed parameters block to perform calculations and implement the ODE solver, for calculating expected values.

Looking at the vignette of this package, it remains unclear to me:

1. If I can create a Functions block where I define the ODE systems, its solver and various utility functions.
2. If I can create a Transformed Parameters block.

Edit 25/6: I see that some functions to generate extra blocks are found here; any advice would be helpful.

Thanks, Lampros

[jtimonen]

Hi, so how this works is you define your ODE system, model parameters, etc as arguments to the ode_model() function, and the Stan code blocks are then generated using those functions you mentioned. Those internal functions are not for users, the only public API to create the model is the ode_model() function.

If I can create a Functions block where I define the ODE systems, its solver and various utility functions.

The functions block is generated internally, the ODE system is defined using the odefun_body argument to ode_model(). The Stan code that gets written includes all solvers that are currently supported (euler, midpoint, rk4, rk45, bdf) and which one you actually use is selected when you call sample. There is unfortunately currently no way to use a user-specified solver or other user-specified functions.

If I can create a Transformed Parameters block.

The transformed parameters block is always there and it will always include at least something like

transformed parameters {
  array[N] vector[D] y_sol_tpar;
  real log_lik_tpar;
  y_sol_tpar = solve_ode(solver, rel_tol, abs_tol, max_num_steps, num_steps,
                         y0, t0, t, ...);
  log_lik_tpar = log_likelihood(y_sol_tpar, ...);
}

and additionally every StanTransformation whose "origin" is "parameters" will be added there. You can generate these with stan_transform(..., origin = "parameters", ...), see the documentation of stan_transform() etc.

[jtimonen]

There should probably be more vignettes. If you wish to have some specific solver included in the package then let me know.