continueEstimation: Add extra iterations to burnin or output.
In StatisticsNZ/demest: Bayesian Demographic Estimation and Forecasting

continueEstimation

R Documentation

Add extra iterations to burnin or output.

Description

Continue estimation process, retaining current settings, but extending the burnin or sampling from the posterior distribution. continueEstimation is called after estimateModel, estimateCounts, or estimateAccount, and can be called multiple times.

Usage

continueEstimation(
  filename,
  nBurnin = NULL,
  nSim = 1000,
  nThin = NULL,
  scaleNoise = 0,
  parallel = NULL,
  outfile = NULL,
  verbose = FALSE,
  useC = TRUE
)

Arguments

`filename`	The filename used by the original call.
`nBurnin`	Number of iteration discarded before recording begins.
`nSim`	Number of additional iterations.
`nThin`	Thinning interval.
`scaleNoise`	Governs noise added to Metropolis-Hastings ratio when updating accounts. Should only be used non-zero when generating initial values. Currently experimental, and may change.
`parallel`	Logical. If `TRUE` (the default), parallel processing is used.
`outfile`	Where to direct the ‘stdout’ and ‘stderr’ connection output from the workers when parallel processing. Passed to function `[parallel]{makeCluster}`.
`verbose`	Logical. If `TRUE` (the default) a message is printed at the end of the calculations.
`useC`	Logical. If `TRUE` (the default), the calculations are done in C. Setting `useC` to `FALSE` may be useful for debugging.

Details

The treatment of output from previous calls to the estimation functions or to continueEstimation depends on nBurnin in the current call to continueEstimation. If nBurnin is NULL, then continueEstimation adds iterations to any existing posterior sample. If nBurnin is a number, then iterations from any previous calls to the estimate function or to continueEstimation are treated as burnin, and the construction of the posterior sample begins again from scratch. See below for an example.

If nThin is set to a different value from previous calls to the estimate function or continueEstimation, then the thinning ratio, and hence correlations between successive iterations, will change. This is safe if the previous iterations are being used as burnin, but needs to be done with care if the previous iterations will form part of the posterior sample.

Because model output includes the state of the random number generator, it should be possible to obtain identical results by (i) calling an estimation function followed by one or more calls to continueEstimation, and (ii) doing all the calculations in one call to an estimaton function. For instance,

estimateModel with nBurnin = 2000 and nSim = 0
continueEstimation with nSim = 1000
continueEstimation with nSim = 1000

should give the same results as

estimateModel with nBurnin = 2000 and nSim = 2000.

Note that the total size of the posterior sample depends not just on nSim, but also on nChain and nThin. In the simplest case (ie a single call to an estimation function, and nSim * nChain divisible by nThin) the number of iterations in the sample equals nSim * nChain / nThin.

References

continueEstimation was inspired by the discussion of checkpointing in Geyer, C. 2011. Introduction to Markov chain Monte Carlo. Brooks, S., Gelman, A., Jones, G. L., and Meng, X-L. (eds.) Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC.

Examples

## prepare data
deaths <- demdata::VADeaths2
popn <- demdata::VAPopn
deaths <- round(deaths)
deaths <- Counts(deaths)
popn <- Counts(popn)

## Show that estimation using a single call to
## 'estimateModel' gives the same results as
## a call to 'estimateModel' followed by a call
## to 'continueEstimation'

## estimate all at once
set.seed(1)
filename.all.at.once <- tempfile()
estimateModel(Model(y ~ Poisson(mean ~ age)),
              y = deaths,
              exposure = popn,
              filename = filename.all.at.once,
              nBurnin = 20,
              nSim = 20,
              nThin = 2,
              nChain = 2,
              parallel = FALSE)
rates.all.at.once <- fetch(filename.all.at.once,
                           where = c("model", "likelihood", "rate"))

## estimate using 'continueEstimation'
set.seed(1)
filename.with.continue <- tempfile()
estimateModel(Model(y ~ Poisson(mean ~ age)),
              y = deaths,
              exposure = popn,
              filename = filename.with.continue,
              nBurnin = 20,
              nSim = 0,
              nThin = 2,
              nChain = 2,
              parallel = FALSE)
continueEstimation(filename = filename.with.continue,
                   nSim = 20)
rates.with.continue <- fetch(filename.all.at.once,
                           where = c("model", "likelihood", "rate"))

## the two approaches give the same answer
all.equal(rates.all.at.once, rates.with.continue)


## Demonstrate the differences between nBurnin = 0
## and nBurnin > 0.

## nBurnin is NULL in call to continueEstimation, so keep 
## iterations from original call to estimateModel
filename.keep.original <- tempfile()
estimateModel(Model(y ~ Poisson(mean ~ age)),
              y = deaths,
              exposure = popn,
              filename = filename.keep.original,
              nBurnin = 20,
              nSim = 10,
              nThin = 2,
              nChain = 2,
              parallel = FALSE)
continueEstimation(filename = filename.keep.original,
                   nSim = 20)
fetchSummary(filename.keep.original) # see 'nBurnin' and 'nSim'

## nBurnin non-NULL in call to continueEstimation, so treat
## iterations from original call to estimateModel
## as part of burnin
filename.discard.original <- tempfile()
estimateModel(Model(y ~ Poisson(mean ~ age)),
              y = deaths,
              exposure = popn,
              filename = filename.discard.original,
              nBurnin = 20,
              nSim = 10,
              nThin = 2,
              nChain = 2,
              parallel = FALSE)
continueEstimation(filename = filename.discard.original,
                   nBurnin = 10,
                   nSim = 10)
fetchSummary(filename.discard.original) # see 'nBurnin' and 'nSim'

StatisticsNZ/demest documentation built on Nov. 2, 2023, 7:56 p.m.