Create PMMH objects

The function biips_pmmh_init initializes the Particle Marginal Metropolis-Hastings (PMMH) algorithm.

The PMMH algorithm is a particle MCMC algorithm using SMC within a MCMC algorithm. It splits the variables in the graphical model into two sets:

The parameters in param_names are sampled with a MH proposal (multivariate Gaussian random walk). They must be continuous and must be defined as single stochastic nodes in the model.
The other latent variables are sampled using a SMC algorithm. They can be monitored using the latent_names parameter.

biips_pmmh_init(model, param_names, latent_names = c(), inits = list(),
  rw_step = list(), transform = TRUE, n_rescale = 400, alpha = 0.99,
  beta = 0.05)

Arguments

model	`biips` model object as returned by `biips_model`.
param_names	character vector. The names of the variables updated with MH proposal. The 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]')`.
latent_names	character vector. The names of the variables to be updated with SMC proposal that need to be monitored.
inits	list of numeric values. Initial values for the parameters in `param_names`. (default = samples from the prior distribution). `inits` can either be a named list with no unnamed member or an unnamed list of the same length as `param_names`.
rw_step	list of numeric values. Random walk standard deviations. `rw_step` can either be a named list with no unnamed member or an unnamed list of the same length as `param_names`. If `transform` is `TRUE` (default), the given standard deviations apply to the transformed variables. Set `transform` to `FALSE` if you wish to set the random walk standard deviation for the untransformed random variables.
transform	logical. Activate automatic parameters transformation (default = `TRUE`). Transformations apply independently to each component of the parameters depending on their support: unbounded (-\(\infty\), +\(\infty\)): f(x) = x lower bounded [L, +\(\infty\)): f(x) = log(x-L) upper bounded (-\(\infty\), U]: f(x) = log(U-x) lower-upper bounded [L, U]: f(x) = log((x-L)/(U-x)) so that we apply a random walk on the real-valued transformed vector.
n_rescale	integer. Number of iterations for rescaling (default = 400).
alpha	real in [0,1]. Tuning parameter of the rescaling (default = 0.99).
beta	real in [0,1]. Tuning parameter of the proposal (default = 0.05).

Value

The function biips_pmmh_init returns an object of class pmmh which can be used to generate samples from the posterior distribution of the monitored variables in param_names and latent_names.

An object of class pmmh 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.

param_names()

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

latent_names()

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

sample_param(sample)

Get and set the current sample value of the parameters.

sample_param_tr(sample)

Get and set the current sample value of the transfromed parameters.

sample_latent(sample)

Get and set the current sample value of the latents.

log_prior(log_prior)

Get and set the current value of the log prior of the parameters.

log_marg_like(log_marg_like)

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

n_iter()

Get and set the current curent number of iterations.

n_rescale()

Get the nb of iterations for rescaling.

rw_dim()

Get the dimensions of the parameters.

rw_len()

Get the total length of the parameters.

rw_step()

Get the random walk standard deviations for the rescaling.

rw_cov()

Get the empirical covariance matrix of the samples.

rw_alpha()

Get the tuning parameter of the rescaling.

rw_beta()

Get the tuning parameter of the proposal.

rw_transform(sample, funtype = 'transform')

Applies transformation function defined by the funtype string argument. Possible values are 'transform' (direct transformation), 'inverse' (inverse transformation) and 'lderiv' (log derivative of the transformation).

rw_rescale(ar)

rescales the random walk step to target an optimal acceptance rate.

rw_learn_cov()

update the empirical covariance matrix of the 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))
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: 168
#>   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_pmmh <- biips_pmmh_init(model, 'logtau', latent_names = c('x', 'c[2:10]'),
                            inits = list(logtau = -2))  # Initialize
#> * Initializing PMMH
is.pmmh(obj_pmmh)
#> [1] TRUE
out_pmmh_burn <- biips_pmmh_update(obj_pmmh, 100, n_part)  # Burn-in
#> * Adapting PMMH with 50 particles
#>   |--------------------------------------------------| 100%
#>   |++++++++++++++++++++++++++++++++++++++++++++++++++| 100 iterations in 0.34 s
out_pmmh <- biips_pmmh_samples(obj_pmmh, 100, n_part, thin = 1)  # Samples
#> * Generating 100 PMMH samples with 50 particles
#>   |--------------------------------------------------| 100%
#>   |**************************************************| 100 iterations in 0.27 s

Arguments

Value

See also

Examples

Contents