Create a PIMH object.

Source: R/pimh.r

The function biips_pimh_init initializes the Particle Independent Metropolis-Hastings (PIMH) algorithm.

The PIMH algorithm provides MCMC samples of the variables in variable_names, using a SMC algorithm as proposal distribution in an independent Metropolis-Hastings (MH) algorithm.

biips_pimh_init(model, variable_names)

Arguments

model	biips model object as returned by biips_model.
variable_names	character vector. The names of the unobserved variables to monitor. Names can contain subset indices which must define a valid subset of the variables of the model, e.g.: c('var1', 'var2[1]', 'var3[1:10]', 'var4[1, 5:10, 3]').

Value

The function biips_pimh_init returns an object of class pimh which can be used to generate samples from the posterior distribution of the monitored variables in variable_names.

An object of class pimh is a list of functions that share a common environment. These functions are meant for internal purpose only. They are used to query information on the current state of the algorithm.

model()

Get the biips model object.

variable_names()

Get a character vector with the names of the monitored variables.

sample(sample)

Get and set the current sample.

log_marg_like(log_marg_like)

Get and set the current value of the log marginal likelihood.

See also

biips_model, biips_pimh_update, biips_pimh_samples

Examples

modelfile <- system.file('extdata', 'hmm.bug', package = 'rbiips')
stopifnot(nchar(modelfile) > 0)
cat(readLines(modelfile), sep = '\n')
#> var c_true[tmax], x_true[tmax], c[tmax], x[tmax], y[tmax]
#> 
#> data
#> {
#>   x_true[1] ~ dnorm(0, 1/5)
#>   y[1] ~ dnorm(x_true[1], exp(logtau_true))
#>   for (t in 2:tmax)
#>   {
#>     c_true[t] ~ dcat(p)
#>     x_true[t] ~ dnorm(0.5*x_true[t-1]+25*x_true[t-1]/(1+x_true[t-1]^2)+8*cos(1.2*(t-1)), ifelse(c_true[t]==1, 1/10, 1/100))
#>     y[t] ~ dnorm(x_true[t]/4, exp(logtau_true))
#>   }
#> }
#> 
#> model
#> {
#>   logtau ~ dunif(-3, 3)
#>   x[1] ~ dnorm(0, 1/5)
#>   y[1] ~ dnorm(x[1], exp(logtau))
#>   for (t in 2:tmax)
#>   {
#>     c[t] ~ dcat(p)
#>     x[t] ~ dnorm(0.5*x[t-1]+25*x[t-1]/(1+x[t-1]^2)+8*cos(1.2*(t-1)), ifelse(c[t]==1, 1/10, 1/100))
#>     y[t] ~ dnorm(x[t]/4, exp(logtau))
#>   }
#> }

data <- list(tmax = 10, p = c(.5, .5), logtau_true = log(1), logtau = log(1))
model <- biips_model(modelfile, data)
#> * Parsing model in: /home/adrien/Dropbox/workspace/rbiips/inst/extdata/hmm.bug
#> * Compiling data graph
#>   Declaring variables
#>   Resolving undeclared variables
#>   Allocating nodes
#>   Graph size: 169
#>   Sampling data
#>   Reading data back into data table
#> * Compiling model graph
#>   Declaring variables
#>   Resolving undeclared variables
#>   Allocating nodes
#>   Graph size: 180

n_part <- 50
obj_pimh <- biips_pimh_init(model, c('x', 'c[2:10]'))  # Initialize
#> * Initializing PIMH
is.pimh(obj_pimh)
#> [1] TRUE
out_pimh_burn <- biips_pimh_update(obj_pimh, 100, n_part)  # Burn-in
#> * Updating PIMH with 50 particles
#>   |--------------------------------------------------| 100%
#>   |**************************************************| 100 iterations in 0.24 s
out_pimh <- biips_pimh_samples(obj_pimh, 100, n_part)  # Samples
#> * Generating PIMH samples with 50 particles
#>   |--------------------------------------------------| 100%
#>   |**************************************************| 100 iterations in 0.26 s