@bbbales2 & @wds15 have lamented for quite a while the lack of application of ode_adams in their modeling & development practice. This is partially my fault in being terse in documentation and later procrastination, and partially due to the bias of models used in Stan’s testing & benchmark, which is really a bias of fields that Stan has touched so far. I’ve started amend the issue in documentation, and here would like to give an example in which ode_adams triumphs.
So far Stan’s ODE development routinely involves only a handful of examples: predator-prey model by @Bob_Carpenter , epiomology model as found in posteriordb, a PKPD model,, a textbook oscillation model, and Lorenz attractor model. The model below is not really more sophisticated, in fact, one can argue it is simpler for being linear:
This is Airy’s equation, and in code it’s just
functions {
vector airy(real t, vector y) {
return [ y[2], -t * y[1] ]';
}
}
Assuming we want to infer y_0 given y observed at t=90,91,\dots,100:
data {
int n; // 11 data points from t=90 to 100
real ts[n];
vector[n] yobs;
real rtol; // 1e-10
real atol; // 1e-10
}
and also want to ensure the accuracy at these points is at scale of 1e-6. Consider the equation has frequent applied in fields optics and electromagnetics, this is not too much to ask for instrumentation. Since it has analytical solution, we can backtrace to determine that tolerance required should be at least 1e-10 in ODE solvers. This happens to be the default of BDF & Adams solvers.
parameters {
real<lower=0> y0;
real<lower=0> sigma;
}
transformed parameters {
vector[2] yinit = [y0, 0.0]';
vector[2] y[n] = ode_adams_tol(airy, yinit, 0.0, ts, rtol, atol, 1000000); // or bdf or rk45
}
model {
y0 ~ cauchy(0, 0.5);
sigma ~ normal(0, 0.01);
for (i in 1:n) {
yobs[i] ~ normal(y[i][1], sigma);
}
}
Simulated data based on y0=1.0 sigma=0.01:
n <- 11
ts <- c(90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100)
yobs <- c(
-0.27204644, 0.26198631, -0.25835151, 0.22486253, -0.16963251, 0.09098462,
0.02549337, -0.14729033, 0.23571404, -0.27632096, 0.25950532)
rtol <- 1.e-10
atol <- 1.e-10
Fitting from init y0=1.1; sigma=0.1, regardless integrator used, we get something like
However, the time performance among the three integrators are very different:
with adams takes less than half of the time than bdf. To see why, we print out the order used in each Adams time step(keep in mind BDF & Adams are variant-order, with BDF order 1-5 and Adams order 1-12), in a fixed_param run:
> fit <- mod$sample(data = "airy.data.R", init="airy.init.R", iter_sampling=1, fixed_param=TRUE, seed=987)
Running MCMC with 1 chain...
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 9
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
Chain 1 Last order: 7
This shows Adams is constantly using much higher orders than rk45 and the max order of BDF, so it can take less steps to reach a same accuracy. It is this capability of bumping up to high orders in smooth non-stiff problems that makes Adams suitable to problems that requires higher accuracy, such as in orbital dynamics and oceanography. Some simulations in these fields require a long-duration time(where is Mercury in year 2050?) and one must reduce tolerance to avoid accumulation of local errors.
The model & data in this post can be found in Torsten’s example models.