Matching the results of a STAN ODE model with a different ODE solver?

[David Manheim]

Hi,

I'm working on an ODE-based model, and wanted to validate it. The best practice is to create a data generator, and so I attempted to do so, but I can't replicate the dynamics of an ODE model.

As a basic example, I used DeSolve to try to replicate the basic dynamics of Michael Betancourt's model, posted in another thread; it is attached. I run the model, and then use the predicted values to populate the DESolve model, and it doesn't match.

Is there a reason that I see large differences in the models? Is it just numeric inaccuracy, and if so, how can I reduce the STAN model's numeric error? (Or is there a bug in my code?)

Thanks,

David

[Michael Betancourt]

I have no idea of the validity of DESolve, but I simulated data in Stan and recovered

it very accurately using this exact model.  I’m guessing there’s a misspecification of

the model between the two programs.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<Replicate Outbreak in R.R><fit_outbreak.stan>

[David Manheim]

I'd like to at least be able to validate the STAN ODE solver with some external ODE solver; has that been done anywhere?

[Daniel Lee]

Hi David,

A few comments:

- The script you attached doesn't provide enough pieces to work through your example. If you'd really like practical help, please provide something that runs.

- As Michael pointed out, this is probably due to a misspecification across programs. What's with this line? "Fit_Outbreak_Replicate_Thetas['kappa'] = 1000000"

- To answer your question, we've validated the ode integration against analytic solutions. They're in our unit tests of the math library.

Daniel

[Bob Carpenter]

And it's not our ODE solver --- it's Boost's. The
thing to test is that we're computing derivatives the
right way.

I've also matched it against the SoilR soil-carbon
modeling package which solves diff eqs using deSolve.

We've also validated in other non-public projects vs. deSolve
itself.

- Bob

[David Manheim]

1) If you run the R script I'm now attaching before running the R script I included earlier, it will run.

2) It matches the stan program line 47; "kappa <- 1000000;"
3) Thanks. I'll try to look more closely at why I can't get this example to replicate more closely, and see if there is something I can do to make the fit closer. 

[David Manheim]

A bit of progress; 

The initialization for STAN is set explicitly (here to 0,) but that for DESOLVE is set to the first time period listed in the time vector (here, 1.)

Fixing that reduces, but does not eliminate, the discrepancy. I'm still looking for more reasons for the mismatch.

[Bob Carpenter]

I'm still not clear on what you're trying to do.

You sent us a Stan model that involves an ODE, whereas deSolve is
a diff eq solver (among other things). What are you trying
to match?

To just run Stan's diff eq solver, you can have an empty model and
parameters and just run the generated quantities (there's a special
flag for models with no parameters). That will let you just run the
solver without the associated derivatives and gradients and Bayesian
posterior fitting.

- Bob

[David Manheim]

I'm trying to match the time series of measurements by using the parameters; if the input parameters match, the output time series should match.

In this case, I can generate whatever STAN results I want, perhaps using an empty model as suggested, but when I input the same initial state and transition parameters into a different ODE solver, I should get the same time series as a result. I don't right now, which tells me that I'm specifying something wrong. Since I'm more comfortable with the (simpler) DeSolve setup, I'm guessing it's the STAN model that I'm mis-understanding, and therefore mis-specifying. I'd like to be able to match the results.

Once I can do this, as a next step I'm planning on using the data generator specified in R to make sure my model converges when excluding certain features in the actual data that cause identifiability issues. (Silly simple example: If I have three quantities that are multiplied to generate the data, I should be able to .estimate their product, even if I can't estimate the individual values.

Is that clearer? Am I confused somewhere about how I should go about this?

[Daniel Lee]

Hi David,

To start, send over the current set of scripts that you expect to match. Given how much time that's passed since the first version, maybe there's a difference.

On Mon, Nov 30, 2015 at 4:04 PM, David Manheim <davidm...@gmail.com> wrote:

I'm trying to match the time series of measurements by using the parameters; if the input parameters match, the output time series should match.

In this case, I can generate whatever STAN results I want, perhaps using an empty model as suggested, but when I input the same initial state and transition parameters into a different ODE solver, I should get the same time series as a result. I don't right now, which tells me that I'm specifying something wrong. Since I'm more comfortable with the (simpler) DeSolve setup, I'm guessing it's the STAN model that I'm mis-understanding, and therefore mis-specifying. I'd like to be able to match the results.

That's a good start. Have you tried a simpler example?

 

Once I can do this, as a next step I'm planning on using the data generator specified in R to make sure my model converges when excluding certain features in the actual data that cause identifiability issues. (Silly simple example: If I have three quantities that are multiplied to generate the data, I should be able to .estimate their product, even if I can't estimate the individual values.

Is that clearer? Am I confused somewhere about how I should go about this?

No, but what you're describing sounds like a misspecification between Stan and deSolve.

I would go to a simple example, maybe even a toy example, and build up. Have you successfully matched deSolve and Stan on any model? Or an analytic model? We've done that exercise, so start there.

[Bob Carpenter]

I don't think there's anything at all tricky about it. You
just need to write the diff eqs as a Stan function and provide
the data and parameters as data. Then you'll get the same answers
to within tolerance as any other working ODE solver (at least for
non-stiff cases Stan can solve).

- Bob

P.S. It's "Stan", not "STAN", unless you're shouting.