nimble-dev / nimble

The base NIMBLE package for R
http://R-nimble.org
BSD 3-Clause "New" or "Revised" License
154 stars 23 forks source link

Request: Microsimulation Markov Chain Monte Carlo #1305

Closed stephematician closed 1 week ago

stephematician commented 1 year ago

It looks like it would be relatively trivial to support Metropolis-Hastings random-walk algorithms that can be used to calibrate 'microsimulation' models. The canonical approach to this follows that of Rutter et al (2009). You can read the paper for an example, but in brief the modification to the random walk is:

Why the "additional" calculation of the log probabilities you ask? Because in these problems one or more parameters of the distributions are estimated by invoking a microsimulation model that produces values that contain small random errors. If the log probabilities aren't re-evaluated, then the chain is doomed whenever the errors in the microsimulation approximation lead to a spuriously high probability of the data.

Anyway, my proposal would be to add an option to the Metropolis-Hastings RW samplers, something like recycleLogProb (default TRUE for the usual behaviour, otherwise FALSE). Then within the run() method of the RW sampler:

     # >>>
        if (!recycleLogProb) invisible(model$calculate(target) + model$calculate(calcNodesNoSelf))
     # <<<
        model[[target]] <<- propValue
        logMHR <- model$calculateDiff(target)
        if(logMHR == -Inf) {
...

With a similar modification for the RW_block sampler.

This does slow down the sampler a lot but it is considered (for now) a necessary evil when calibrating these types of models using MCMC.

Thoughts welcome!

Rutter CM, Miglioretti DL, Savarino JE. Bayesian Calibration of Microsimulation Models. J Am Stat Assoc. 2009 Dec 1;104(488):1338-1350. doi: 10.1198/jasa.2009.ap07466

perrydv commented 1 year ago

Hi @stephematician, thanks for the suggestion! A few thoughts for you:

  1. You can easily prototype the idea by copying the source code for the built-in samplers into a working file, change the name argument, modify as you wish, and use them as user-defined samplers.
  2. We'd have to think about whether it is a general enough need to put into the main samplers vs separate samplers. Granted, the change is very small.
  3. It seems somewhat related to Particle MCMC, which we have implemented and describe here, in case that is relevant. A possible difference is that the particle filter simulated log likelihood of PMCMC uses the model object itself. We use a sampler that allows an arbitrary likelihood calculation to call the particle filter.