pimp my pomp

Tips, tricks, advice, and examples.

• Keeping a database of parameter-space explorations
• Reproducibility on a multicore machine via bake, stew, and freeze
• How can I include a vector of variables in a Csnippet?
• How to handle missing data
• How to deal with accumulator variables when t0 is much less than t[1]

Keeping a database of parameter-space explorations

Likelihood surfaces for dynamic models can be very complex and the computations needed to explore them can be expensive. By keeping a record of all parameter points visited, along with the computed likelihood at each point, is a good way to ensure that you continually improve your picture of the likelihood surface.

Doing this can be as simple as maintaining a CSV file with one column for each parameter, plus the likelihood (and s.e.). It can be useful to supplement this with an indication of the name of the model and any other qualifying information.

Reproducibility on a multicore machine via bake, stew, and freeze

It is often the case that heavy pomp computations are best performed in parallel on a cluster or multi-core machine. This poses some challenges in trying to ensure reproducibility and avoiding repetition of expensive calculations. The bake, stew, and freeze functions provide some useful facilities in this regard.

For example:

require(pomp)
pompExample(ricker)

require(foreach)
require(doMC)
registerDoMC(5)

bake(file="pfilter1.rds",seed=607686730,kind="L'Ecuyer",{
  foreach (i=1:10, .combine=c, 
           .options.multicore=list(set.seed=TRUE)) %dopar% {
     pf <- pfilter(ricker,Np=1000)
     logLik(pf)
   }
}) -> ll

In the above bake first checks to see whether the file pfilter1.rds exists. If it does, it then loads it (using readRDS) and stores the result in ll. If it does not, it evaluates the expression embraced in the brackets, stores the result in pfilter1.rds, and returns it. While the expression is evaluated, the R session's pseudorandom number generator (RNG) is temporarily set to the state specified by the seed and kind arguments (see ?set.seed). In this case, the expression to be evaluated makes use of the foreach and doMC packages to run 10 particle filtering operations in parallel on a multicore machine. We therefore use a parallel RNG ("L'Ecuyer").

The bake function stores or retrieves and returns a single R object. If one wants to produce multiple objects in a reproducible way, use stew. For example:

require(pomp)
pompExample(ricker)

stew(file="pfilter2.rda",seed=607686730,kind="L'Ecuyer",{
  te <- system.time(
  foreach (i=1:10, .combine=c, 
           .options.multicore=list(set.seed=TRUE)) %dopar% {
      pf <- pfilter(ricker,Np=1000)
      logLik(pf)
   } -> ll2
  )
})

In the above, stew again temporarily sets the RNG state before evaluating the expression. The objects te and ll2 are created during this evaluation; these are stored in pfilter2.rda to be retrieved if the snippet is run a second time.

Like bake and stew, the freeze function temporarily sets the state of the RNG and evaluates an arbitrary R expression, and finally returns the RNG state to its status quo ante. Unlike bake and stew, freeze neither stores nor retrieves results.

How can I include a vector of variables in a Csnippet?

See this entry in the FAQ.

How to handle missing data

See FAQ 3.3. Thanks to @christelpei for the question and the answer.

How to deal with accumulator variables when t0 is much less than t[1]

See FAQ 3.4. Again, thanks to @christelpei for the question and the answer.