pmmh-object: Manipulate PMMH objects.
In biips/rbiips: Bayesian Inference with Interacting Particle Systems

Description Usage Arguments Details Value See Also Examples

Manipulate PMMH objects.

The method biips_pmmh_update performs adaptation and burn-in iterations for the PMMH algorithm.

The method biips_pmmh_samples performs iterations for the PMMH algorithm and returns samples.

is.pmmh(object)

## S3 method for class 'pmmh'
biips_pmmh_update(object, n_iter, n_part, thin = 1,
  max_fail = 0, rw_adapt = TRUE, output = "p", ...)

## S3 method for class 'pmmh'
biips_pmmh_samples(object, n_iter, n_part, thin = 1,
  max_fail = 0, output = "p", ...)

`object`	a `pmmh` object as returned by `biips_pmmh_init`.
`n_iter`	integer. Number of burn-in iterations.
`n_part`	integer. Number of particles used in SMC algorithms.
`thin`	integer. Thinning interval. Returns samples every `thin` iterations (default = 1)
`max_fail`	integer. maximum number of failed SMC algorithms allowed. (default=0).
`rw_adapt`	logical. Activate adaptation of the proposal (default=`TRUE`).
`output`	string. Select additional members to be returned in the `mcmcarray.list` output. The string can contain several characters in `('p', 'l', 'a', 's', 'f')`. See details. (default = `'p'`)
`...`	Additional arguments to be passed to the SMC algorithm such as `rs_thres` and `rs_type`. See `biips_smc_samples` for more details.

The output string arguments can be used to query additional members in the mcmcarray.list output. If output contains:

p: returns member log_marg_like_pen. mcmcarray with penalized log marginal likelihood estimates over iterations.
l: returns member log_marg_like. mcmcarray with log marginal likelihood estimates over iterations.
a: returns member info$accept_rate. mcmcarray with acceptance rate over iterations.
s: returns member info$rw_step. mcmcarray with standard deviations of the random walk over iterations.
f: returns member info$n_fail. number of failed SMC algorithms.

The function is.pmmh returns TRUE if the object is of class pmmh.

The methods biips_pmmh_update and biips_pmmh_update return an object of class mcmcarray.list.

biips_pmmh_samples output contains one mcmcarray member for each monitored variable returned by the param_names() and latent_names() member functions of the pmmh object.

The members of the mcmcarray.list object are mcmcarray objects for different variables. Assuming dim is the dimension of the monitored variable, the mcmcarray object is an array of dimension c(dim, n_iter) with the following attributes (accessible with attr):

`name`	string with the name of the variable.
`lower`	vector with the lower bounds of the variable.
`upper`	vector with the upper bounds of the variable.

If the output argument is not empty, the output contains additional members. See details.

biips_pmmh_init

modelfile <- system.file('extdata', 'hmm.bug', package = 'rbiips')
stopifnot(nchar(modelfile) > 0)
cat(readLines(modelfile), sep = '\n')

data <- list(tmax = 10, p = c(.5, .5), logtau_true = log(1))
model <- biips_model(modelfile, data)

n_part <- 50
obj_pmmh <- biips_pmmh_init(model, 'logtau', latent_names = c('x', 'c[2:10]'),
                            inits = list(logtau = -2))  # Initialize
is.pmmh(obj_pmmh)
out_pmmh_burn <- biips_pmmh_update(obj_pmmh, 100, n_part)  # Burn-in
out_pmmh <- biips_pmmh_samples(obj_pmmh, 100, n_part, thin = 1)  # Samples

dens_pmmh_lt <- biips_density(out_pmmh$logtau)
summ_pmmh_x <- biips_summary(out_pmmh$x, order = 2, probs = c(0.025, 0.975))
dens_pmmh_x <- biips_density(out_pmmh$x)
summ_pmmh_c <- biips_summary(out_pmmh[['c[2:10]']])
table_pmmh_c <- biips_table(out_pmmh[['c[2:10]']])

par(mfrow = c(2, 2))
plot(c(out_pmmh_burn$log_marg_like_pen, out_pmmh$log_marg_like_pen), type = 'l',
     col = 'blue', xlab = 'PMMH iteration', ylab = 'log p(y|logtau) + log p(logtau)')

plot(out_pmmh$logtau[1, ], type = 'l', col = 'blue',
     xlab = 'PMMH iteration', ylab = 'logtau')
points(0, model$data()$logtau_true, pch = 17, col = 'green')

plot(dens_pmmh_lt, col = 'blue', ylab = 'posterior density')
points(model$data()$logtau_true, 0, pch = 17, col = 'green')

biips_hist(out_pmmh$logtau, col = 'blue', ylab = 'posterior density')
points(model$data()$logtau_true, 0, pch = 17, col = 'green')

par(mfrow = c(2, 2))
plot(model$data()$x_true, type = 'l', col = 'green', xlab = 't', ylab = 'x[t]')
lines(summ_pmmh_x$mean, col = 'blue')
matlines(matrix(unlist(summ_pmmh_x$quant), data$tmax), lty = 2, col = 'blue')
legend('topright', leg = c('true', 'PMMH estimate'), lty = c(2, 1),
       col = c('green', 'blue'), bty = 'n')

barplot(.5*(model$data()$c_true==1), col = 'green', border = NA, space = 0, offset = 1,
        ylim=c(0,2), xlab='t', ylab='c[t]==1', axes = FALSE)
axis(1, at=1:data$tmax-.5, labels=1:data$tmax)
axis(2, line = 1, at=c(0,2), labels=NA)
text(data$tmax/2, 1.75, 'true')
barplot(.5*c(NA, summ_pmmh_c$mode==1), col = 'blue', border = NA, space = 0,
        axes = FALSE, add = TRUE)
text(data$tmax/2, .75, 'PMMH mode')

t <- 5
plot(dens_pmmh_x[[t]], col='blue', ylab = 'posterior density')
points(model$data()$x_true[t], 0, pch = 17, col = 'green')

plot(table_pmmh_c[[t-1]], col='blue', ylab = 'posterior probability mass')
points(model$data()$c_true[t], 0, pch = 17, col = 'green')