#' @keywords internal
pmmh_set_param <- function(object, ...) UseMethod("pmmh_set_param")
#' @keywords internal
pmmh_set_param.biips <- function(object, param_names, pn_param, values) {
stopifnot(is.biips(object))
console <- object$ptr()
n_param <- length(pn_param$names)
sample_param <- list()
if (length(values) > 0) {
stopifnot(length(values) == n_param)
# take sample values in values parameter
for (i in 1:n_param) {
var <- param_names[[i]]
stopifnot(!is.null(values[[var]]))
ok <- rbiips("change_data", console, pn_param$names[[i]], pn_param$lower[[i]],
pn_param$upper[[i]], values[[var]], TRUE)
if (!ok)
stop("Data change failed: invalid value for variable ", var)
sample_param[[var]] <- values[[var]]
}
} else {
# or sample init value
for (i in 1:n_param) {
# tryCatch({
var <- param_names[[i]]
sample_param[[var]] <- rbiips("sample_data", console, pn_param$names[[i]],
pn_param$lower[[i]], pn_param$upper[[i]], get_seed())
# }, error = function(e) {
# print(e)
# stop("Can not sample variable ", var, ": BUG TO BE FIXED")
# })
}
#### FIXME what if the variable is observed ?
}
return(invisible(sample_param))
}
#' @title Create PMMH objects
#'
#' @description The function \code{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:
#' \itemize{
#' \item The parameters in \code{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.
#' \item The other latent variables are sampled using a SMC algorithm. They can
#' be monitored using the \code{latent_names} parameter.
#' }
#'
#' @export
#' @param model \code{biips} model object as returned by
#' \code{\link{biips_model}}.
#' @param 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.: \code{c('var1',} \code{'var2[1]',}
#' \code{'var3[1:10]',} \code{'var4[1, 5:10, 3]')}.
#' @param latent_names character vector. The names of the variables to be
#' updated with SMC proposal that need to be monitored.
#' @param inits list of numeric values. Initial values for the parameters
#' in \code{param_names}. (default = samples from the prior distribution).
#' \code{inits} can either be a named list with no unnamed member or an
#' unnamed list of the same length as \code{param_names}.
#' @param transform logical. Activate automatic parameters transformation
#' (default = \code{TRUE}). Transformations apply independently to each
#' component of the parameters depending on their support: \itemize{ \item
#' unbounded (-\eqn{\infty}, +\eqn{\infty}): f(x) = x \item lower bounded [L,
#' +\eqn{\infty}): f(x) = log(x-L) \item upper bounded (-\eqn{\infty}, U]:
#' f(x) = log(U-x) \item 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.
#' @param rw_step list of numeric values. Random walk standard deviations.
#' \code{rw_step} can either be a named list with no unnamed member or an
#' unnamed list of the same length as \code{param_names}. If \code{transform}
#' is \code{TRUE} (default), the given standard deviations apply to the transformed variables.
#' Set \code{transform} to \code{FALSE} if you wish to set the random walk standard deviation for the
#' untransformed random variables.
#' @param n_rescale integer. Number of iterations for rescaling (default =
#' 400).
#' @param alpha real in [0,1]. Tuning parameter of the rescaling
#' (default = 0.99).
#' @param beta real in [0,1]. Tuning parameter of the proposal (default
#' = 0.05).
#'
#' @return The function \code{biips_pmmh_init} returns an object of class
#' \code{pmmh} which can be used to generate samples from the posterior
#' distribution of the monitored variables in \code{param_names} and
#' \code{latent_names}.
#'
#' An object of class \code{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.
#' \item{model()}{Get the \code{\link{biips}} model object.}
#' \item{param_names()}{Get a character vector with the names of the monitored
#' parameters.}
#' \item{latent_names()}{Get a character vector with the names of
#' the monitored latents.}
#' \item{sample_param(sample)}{Get and set the current sample value of the parameters.}
#' \item{sample_param_tr(sample)}{Get and set the current sample value of the
#' transfromed parameters.}
#' \item{sample_latent(sample)}{Get and set the current sample value of the
#' latents.}
#' \item{log_prior(log_prior)}{Get and set the current value of the
#' log prior of the parameters.}
#' \item{log_marg_like(log_marg_like)}{Get and set the current value of
#' the log marginal likelihood.}
#' \item{n_iter()}{Get and set the current curent number of iterations.}
#' \item{n_rescale()}{Get the nb of iterations for rescaling.}
#' \item{rw_dim()}{Get the dimensions of the parameters.}
#' \item{rw_len()}{Get the total length of the parameters.}
#' \item{rw_step()}{Get the random walk standard deviations for the
#' rescaling.}
#' \item{rw_cov()}{Get the empirical covariance matrix of the samples.}
#' \item{rw_alpha()}{Get the tuning parameter of the rescaling.}
#' \item{rw_beta()}{Get the tuning parameter of the proposal.}
#' \item{rw_transform(sample, funtype = 'transform')}{Applies transformation
#' function defined by the \code{funtype} string argument. Possible values are
#' \code{'transform'} (direct transformation), \code{'inverse'} (inverse
#' transformation) and \code{'lderiv'} (log derivative of the transformation).}
#' \item{rw_rescale(ar)}{rescales the random walk step to target an optimal
#' acceptance rate.}
#' \item{rw_learn_cov()}{update the empirical covariance matrix of the samples.}
#'
#' @seealso \code{\link{biips_model}}, \code{\link{biips_pmmh_update}},
#' \code{\link{biips_pmmh_samples}}
#'
#' @examples
#' 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
biips_pmmh_init <- function(model, param_names, latent_names = c(), inits = list(),
rw_step = list(), transform = TRUE, n_rescale = 400, alpha = 0.99, beta = 0.05) {
## check arguments
stopifnot(is.biips(model))
stopifnot(is.logical(transform), length(transform) == 1)
stopifnot(is.numeric(n_rescale), length(n_rescale) == 1, n_rescale >= 1, is.finite(n_rescale))
stopifnot(is.numeric(alpha), length(alpha) == 1, alpha >= 0, alpha <= 1)
stopifnot(is.numeric(beta), length(beta) == 1, beta >= 0, beta <= 1)
## TODO check param_names and latent_names are valid
## check param_names
stopifnot(is.character(param_names), length(param_names) > 0)
pn_param <- parse_varnames(param_names)
n_param <- length(param_names)
if (any(duplicated(sapply(param_names, to_biips_vname))))
stop("duplicated names in param_names")
## check latent_names
if (length(latent_names) > 0) {
stopifnot(is.character(latent_names))
parse_varnames(latent_names)
if (any(duplicated(sapply(latent_names, to_biips_vname))))
stop("duplicated names in latent_names")
}
## check inits
stopifnot(is.list(inits))
if (length(inits) > 0) {
stopifnot(length(inits) == n_param, sapply(inits, is.numeric))
stopifnot(sapply(inits, function(x) all(is.finite(x))))
init_names <- names(inits)
if (is.null(init_names))
names(inits) <- param_names else {
stopifnot(init_names %in% param_names)
if (any(duplicated(sapply(init_names, to_biips_vname))))
stop("duplicated names in inits")
}
inits <- inits[param_names]
}
## check rw_step
stopifnot(is.list(rw_step))
if (length(rw_step) > 0) {
stopifnot(length(rw_step) == n_param, sapply(rw_step, is.numeric))
stopifnot(sapply(rw_step, function(x) all(is.finite(x) && x > 0)))
step_names <- names(rw_step)
if (is.null(step_names)) {
names(rw_step) <- param_names
} else {
stopifnot(step_names %in% param_names)
if (any(duplicated(sapply(step_names, to_biips_vname))))
stop("duplicated names in rw_step")
}
rw_step <- rw_step[param_names]
}
## stop biips verbosity
verb <- rbiips("verbosity", 0)
on.exit(rbiips("verbosity", verb))
## display message
rbiips("message", "Initializing PMMH")
## Clone console
model <- clone_model(model)
## Init the parameters of the random walk
sample_param <- pmmh_set_param(model, param_names, pn_param, inits)
state <- list(sample_param = sample_param, sample_latent = list(), log_marg_like = -Inf,
log_prior = -Inf)
## store the dimensions of the variables
fun_get_dim <- function(x) dim(as.array(x))
sample_dim <- lapply(sample_param, fun_get_dim)
sample_len <- sum(sapply(sample_dim, prod))
# target_prob = target acceptance probability. The default seems to be a fairly
# robust optimal value. Comes from Rosenthal 2009, Examples of Adaptives MCMC, p
# 16. In one-dimensional case, we take 0.44 which is optimal in certain settings.
target_prob <- if (sample_len == 1)
0.44 else 0.234
## Init random walk stepsize for the part with diagonal covariance matrix
if (length(rw_step) == 0) {
## default values
rw_step <- lapply(sample_dim, function(x) 0.1/sqrt(sample_len) * array(1,
dim = x))
} else {
## check dimensions
rw_step_dim <- lapply(rw_step, fun_get_dim)
stopifnot(mapply(identical, rw_step_dim, sample_dim))
}
## Concatenate all log value
log_step <- log(unlist(rw_step))
# local initial parameters of the algorithm of adaptation of the step of random
# walk in the PMMH algorithm
rw <- list(n_rescale = n_rescale, n_cov = n_rescale/2, n_iter = 0, ar_mean = 0,
log_step = log_step, beta = beta, alpha = alpha, target_prob = target_prob,
dim = sample_dim, len = sample_len, mean = c(), cov = c(), transform = list(),
transform_inv = list(), transform_lderiv = list())
# we start learning the covariance matrix after n_cov iterations
## helper transformation functions
# lower-upper bounded support: logit transform
tr_logit <- function(l,u){ l; u;
function(x) log((x - l)/(u - x))
} # direct transformation
tr_logit_inv <- function(l,u){ l; u;
function(y) l + (u - l)/(1 + exp(-y))
} # inverse transformation
tr_logit_lderiv <- function(l,u){ l; u;
function(x) log(u - l) - log(x - l) - log(u - x)
} # log derivative
# lower bounded support: log transform
tr_log <- function(l){ l; function(x) log(x - l) }
tr_log_inv <- function(l){ l; function(y) l + exp(y) }
tr_log_lderiv <- function(l){ l; function(x) -log(x - l) }
# upper bounded support: -log transform
tr_mlog <- function(u){ u; function(x) log(u - x) }
tr_mlog_inv <- function(u){ u; function(y) u - exp(y) }
tr_mlog_lderiv <- function(u){ u; function(x) -log(u - x) }
## Define variable transformations
k <- 0
for (i in 1:n_param) {
var <- param_names[[i]]
len <- prod(sample_dim[[i]])
if (transform) {
# tryCatch({
support <- rbiips("get_fixed_support", model$ptr(), pn_param$names[[i]],
pn_param$lower[[i]], pn_param$upper[[i]])
# }, error = function(e) {
# print(e)
# stop("CANNOT GET FIXED SUPPORT OF VARIABLE ", var, ": BUG TO BE FIXED")
# })
if (nrow(support) != len)
stop("support size does not match")
} else {
support <- matrix(c(rep(-Inf, len), rep(+Inf, len)), nrow = len)
}
for (j in 1:len) {
k <- k + 1
L <- support[j, 1]
U <- support[j, 2]
if (is.finite(L)) {
if (is.finite(U)) {
# lower-upper bounded support: logit transform
rw$transform[[k]] <- tr_logit(L, U) # direct transformation
rw$transform_inv[[k]] <- tr_logit_inv(L, U) # inverse transformation
rw$transform_lderiv[[k]] <- tr_logit_lderiv(L, U) # log derivative
} else {
if (U < 0)
stop("upper can not be -Inf")
# lower bounded support: log transform
rw$transform[[k]] <- tr_log(L)
rw$transform_inv[[k]] <- tr_log_inv(L)
rw$transform_lderiv[[k]] <- tr_log_lderiv(L)
}
} else {
if (L > 0)
stop("lower can not be +Inf")
if (is.finite(U)) {
# upper bounded support: -log transform
rw$transform[[k]] <- tr_mlog(U)
rw$transform_inv[[k]] <- tr_mlog_inv(U)
rw$transform_lderiv[[k]] <- tr_mlog_lderiv(U)
} else {
if (U < 0)
stop("upper can not be -Inf")
# upper bounded support: identity
rw$transform[[k]] <- function(x) x
rw$transform_inv[[k]] <- function(y) y
rw$transform_lderiv[[k]] <- function(x) 0
}
}
}
}
# function that applies transformation functions to samples
rw_transform <- function(sample, funtype = "transform") {
funlist <- switch(funtype, transform = rw$transform, inverse = rw$transform_inv,
lderiv = rw$transform_lderiv)
# concatenate all variables in a vector always in the order of rw$dim
sample_vec <- as.list(c(unlist(sample)))
out_vec <- mapply(FUN = function(x, fun) fun(x), sample_vec, funlist)
# rearrange vectorized parameter to list of arrays with appropriate dimensions
out <- list()
from <- 1
for (n in names(rw$dim)) {
len <- prod(rw$dim[[n]])
to <- from + len - 1
out[[n]] <- array(out_vec[from:to], dim = rw$dim[[n]])
from <- from + len
}
return(out)
}
state$sample_param_tr <- rw_transform(sample_param)
## output object
obj_pmmh <- list(model = function() {
model
}, param_names = function() {
param_names
}, latent_names = function() {
latent_names
}, sample_param = function(sample) {
if (!missing(sample)) state$sample_param <<- sample
invisible(state$sample_param)
}, sample_param_tr = function(sample_tr) {
if (!missing(sample_tr)) state$sample_param_tr <<- sample_tr
invisible(state$sample_param_tr)
}, sample_latent = function(sample) {
if (!missing(sample)) state$sample_latent <<- sample
invisible(state$sample_latent)
}, log_prior = function(log_prior) {
if (!missing(log_prior)) state$log_prior <<- log_prior
invisible(state$log_prior)
}, log_marg_like = function(log_marg_like) {
if (!missing(log_marg_like)) state$log_marg_like <<- log_marg_like
invisible(state$log_marg_like)
}, n_iter = function(n_iter) {
if (!missing(n_iter)) rw$n_iter <<- n_iter
invisible(rw$n_iter)
}, n_rescale = function() {
rw$n_rescale
}, rw_dim = function() {
rw$dim
}, rw_len = function() {
rw$len
}, rw_step = function() {
exp(rw$log_step)
}, rw_cov = function() {
rw$cov
}, rw_alpha = function() {
rw$alpha
}, rw_beta = function() {
rw$beta
}, rw_transform = rw_transform, rw_rescale = function(ar) {
ar <- min(1, ar)
rw$ar_mean <<- rw$ar_mean + 1 * (ar - rw$ar_mean)/rw$n_iter
rw$log_step <<- rw$log_step + rw$alpha^(rw$n_iter) * (rw$ar_mean - rw$target_prob)
return(invisible())
}, rw_learn_cov = function() {
n_iter <- rw$n_iter
n_cov <- rw$n_cov
## Concatenate all variables in a column vector
sample_vec_tr <- unlist(state$sample_param_tr)
if (n_iter == n_cov + 1) {
rw$mean <<- sample_vec_tr
rw$cov <<- outer(sample_vec_tr, sample_vec_tr)
} else if (n_iter > n_cov + 1) {
## Recursive update of the empirical mean and covariance matrix
q <- (n_iter - n_cov - 1)/(n_iter - n_cov)
q2 <- (n_iter - n_cov - 1)/(n_iter - n_cov)^2
z <- sample_vec_tr - rw$mean
rw$cov <<- q * rw$cov + q2 * outer(z, z)
rw$mean <<- q * rw$mean + 1/(n_iter - n_cov) * sample_vec_tr
}
return(invisible())
})
class(obj_pmmh) <- "pmmh"
return(obj_pmmh)
}
#' Manipulate PMMH objects.
#'
#' @name pmmh-object
#' @aliases pmmh biips_pmmh_update biips_pmmh_samples
#' @param object a \code{pmmh} object as returned by
#' \code{\link{biips_pmmh_init}}.
#' @seealso \code{\link{biips_pmmh_init}}
NULL
#' @rdname pmmh-object
#' @export
#' @return The function \code{is.pmmh} returns \code{TRUE} if the object is of
#' class \code{pmmh}.
is.pmmh <- function(object) {
return(class(object) == "pmmh")
}
#' @keywords internal
#' @importFrom stats runif rnorm
pmmh_rw_proposal <- function(object) {
# concatenate all variables in a vector always in the order of rw_dim
sample_vec_tr <- c(unlist(object$sample_param_tr()))
## Check dimension
rw_len <- object$rw_len()
stopifnot(length(sample_vec_tr) == rw_len)
n_iter <- object$n_iter()
n_rescale <- object$n_rescale()
beta <- object$rw_beta()
cov <- object$rw_cov()
if (length(cov) == 0 || (n_iter < n_rescale) || (runif(1) < beta)) {
# Proposal with diagonal covariance
prop_vec_tr <- sample_vec_tr + object$rw_step() * rnorm(rw_len)
} else {
# proposal with learnt covariance
epsilon <- 1e-05 # For numerical stability
cov_chol <- t(chol(cov + epsilon * diag(1, nrow = rw_len)))
## TODO use pivot in chol ? check positive semi-definite ?
prop_vec_tr <- sample_vec_tr + 2.38/sqrt(rw_len) * cov_chol %*% rnorm(rw_len)
}
# rearrange vectorized parameter to list of arrays with appropriate dimensions
rw_dim <- object$rw_dim()
prop_tr <- list()
from <- 1
for (n in names(rw_dim)) {
len <- prod(rw_dim[[n]])
to <- from + len - 1
prop_tr[[n]] <- array(prop_vec_tr[from:to], dim = rw_dim[[n]])
from <- from + len
}
return(prop_tr)
}
#' Realize one step of the PMMH algorithm.
#' @keywords internal
#' @importFrom stats runif
pmmh_one_update <- function(object, pn_param, n_part, rw_rescale, rw_learn, ...) {
console <- object$model()$ptr()
param_names <- object$param_names()
latent_names <- object$latent_names()
sample_param <- object$sample_param()
sample_param_tr <- object$sample_param_tr()
sample_latent <- object$sample_latent()
log_prior <- object$log_prior()
log_marg_like <- object$log_marg_like()
n_fail <- 0
n_param <- length(param_names)
n_latent <- length(latent_names)
## Increment iterations counter
n_iter <- object$n_iter() + 1
object$n_iter(n_iter)
## Random walk proposal
prop_tr <- pmmh_rw_proposal(object)
prop <- object$rw_transform(prop_tr, "inverse")
lderiv <- object$rw_transform(prop, "lderiv")
### TODO: check NA ?
# for (var in names(prop)) { if (any(is.na(prop[[var]]))) { stop('PMMH proposal
# have NA: ', var, ' = ', prop[[var]]) } }
## compute log prior density
log_prior_prop <- 0
for (i in 1:n_param) {
## change model data
var <- param_names[[i]]
stopifnot(!is.null(prop[[var]]))
ok <- rbiips("change_data", console, pn_param$names[[i]], pn_param$lower[[i]],
pn_param$upper[[i]], prop[[var]], TRUE)
if (!ok) {
## DATA CHANGE FAILED: proposed parameter value might be out of bounds ? TODO:
## double check this is the case and there is no other reason.
log_prior_prop <- -Inf
# warning('Data change failed:', var, '=', prop[[var]], '\n')
break
}
log_p <- rbiips("get_log_prior_density", console, pn_param$names[[i]], pn_param$lower[[i]],
pn_param$upper[[i]])
log_p <- log_p - sum(abs(lderiv[[var]]))
if (is.na(log_p))
stop("Failed to get log prior density: node ", var, " is not stochastic.")
log_prior_prop <- log_prior_prop + log_p
}
if (is.nan(log_prior_prop)) {
# accept_rate = 0
stop("Failed to compute log prior density : ", log_prior_prop)
}
if (log_prior_prop == -Inf) {
## If proposal is not in the support of the prior
accept_rate <- 0
} else {
## compute log marginal likelihood: run smc sampler
ok <- smc_forward_algo(object$model(), n_part = n_part, ...)
if (!ok) {
log_marg_like_prop <- -Inf
n_fail <- n_fail + 1
warning("Failure running smc forward sampler\n")
} else {
log_marg_like_prop <- rbiips("get_log_norm_const", console)
if (is.nan(log_marg_like_prop) || log_marg_like_prop == Inf)
### TODO error or n_fail increment ?
stop("Failed to compute log marginal likelihood: ", log_marg_like_prop)
}
## Acceptance rate
accept_rate <- min(1, exp(log_marg_like_prop - log_marg_like + log_prior_prop -
log_prior))
if (is.nan(accept_rate))
stop("Failed to compute acceptance rate: ", accept_rate)
## Accept-Reject step
if (runif(1) < accept_rate) {
sample_param <- prop
sample_param_tr <- prop_tr
log_prior <- log_prior_prop
log_marg_like <- log_marg_like_prop
if (n_latent > 0) {
## sample one particle for the latent variables
sampled_value <- rbiips("sample_gen_tree_smooth_particle", console,
get_seed())
for (var in latent_names) {
var_in <- to_biips_vname(var)
sample_latent[[var]] <- sampled_value[[var_in]]
}
}
}
}
## Update PMMH object with current state
object$sample_param(sample_param)
object$sample_param_tr(sample_param_tr)
object$sample_latent(sample_latent)
object$log_prior(log_prior)
object$log_marg_like(log_marg_like)
## rescale random walk step
if (rw_rescale)
object$rw_rescale(accept_rate)
## update random walk covariance matrices
if (rw_learn)
object$rw_learn_cov()
return(list(accept_rate = accept_rate, n_fail = n_fail))
}
pmmh_algo <- function(object, ...) UseMethod("pmmh_algo")
#' @keywords internal
pmmh_algo.pmmh <- function(object, n_iter, n_part, return_samples, thin = 1, max_fail = 0,
rw_adapt = FALSE, output = "p", proposal = "prior", ...) {
stopifnot(is.pmmh(object))
stopifnot(is.numeric(n_iter), length(n_iter) == 1, n_iter > 0, is.finite(n_iter))
n_iter <- as.integer(n_iter)
stopifnot(is.numeric(n_part), length(n_part) == 1, n_part >= 1, is.finite(n_part))
n_part <- as.integer(n_part)
stopifnot(is.logical(return_samples), length(return_samples) == 1)
stopifnot(is.numeric(thin), length(thin) == 1, thin >= 1, thin <= n_iter)
thin <- as.integer(thin)
stopifnot(is.numeric(max_fail), length(max_fail) == 1, max_fail >= 0)
max_fail <- as.integer(max_fail)
stopifnot(is.logical(rw_adapt), length(rw_adapt) == 1)
stopifnot(is.character(output))
output <- unlist(strsplit(output, NULL))
output <- match.arg(output, c("p", "l", "a", "f", "s"), several.ok = TRUE)
output <- unique(output)
proposal <- match.arg(proposal, c("auto", "prior"))
## stop biips verbosity
verb <- rbiips("verbosity", 0)
# reset verbosity when function terminates
on.exit(rbiips("verbosity", verb))
## Initialization -------------------------
console <- object$model()$ptr()
# monitor variables
param_names <- object$param_names()
pn_param <- parse_varnames(param_names)
latent_names <- object$latent_names()
pn_latent <- parse_varnames(latent_names)
n_param <- length(param_names)
n_latent <- length(latent_names)
if (n_latent > 0)
monitor(object$model(), latent_names, "s")
# build smc sampler
if (!rbiips("is_sampler_built", console))
rbiips("build_smc_sampler", console, proposal=="prior")
# toggle rescaling adaptation
rw_rescale <- rw_adapt && object$n_iter() < object$n_rescale()
## FIXME if (!rw.adapt) { ## check adaptation if (object$.rw_adapt() &&
## !object$.rw_check_adapt()) message('NOTE: Stopping adaptation of the PMMH
## random walk.\n') ## turn off adaptation object$.rw_adapt_off() }
# set current param value to the model
sample_param <- object$sample_param()
pmmh_set_param(object$model(), param_names, pn_param, sample_param)
# Initialize counters
n_samples <- floor(n_iter/thin)
ind_sample <- 0
n_fail <- 0
# Output structure with MCMC samples
out <- list()
if ("p" %in% output)
out$log_marg_like_pen <- mcmcarray(dim = c(1, n_samples), name = "log_marg_like_pen")
if ("l" %in% output)
out$log_marg_like <- mcmcarray(dim = c(1, n_samples), name = "log_marg_like")
if ("a" %in% output)
out$info$accept_rate <- mcmcarray(dim = c(1, n_samples), name = "accept_rate")
if ("s" %in% output) {
len <- length(object$rw_step())
out$info$rw_step <- mcmcarray(dim = c(len, n_samples), name = "rw_step")
}
if (return_samples) {
for (ivar in 1:length(param_names)) {
var <- param_names[[ivar]]
dimen <- dim(sample_param[[var]])
out[[var]] <- mcmcarray(dim = c(dimen, n_samples),
name = pn_param$names[[ivar]],
lower = pn_param$lower[[ivar]],
upper = pn_param$upper[[ivar]])
}
}
## display message and progress bar
mess <- if (return_samples)
paste("Generating", n_samples, "PMMH samples with") else {
if (rw_adapt)
"Adapting PMMH with" else "Updating PMMH with"
}
symbol <- if (rw_adapt)
"+" else "*"
rbiips("message", paste(mess, n_part, "particles"))
bar <- rbiips("progress_bar", n_iter, symbol, "iterations")
### TODO: display expected time of run
## Metropolis-Hastings iterations -------------------------------
for (i in 1:n_iter) {
out_step <- pmmh_one_update(object, pn_param = pn_param, n_part = n_part,
rw_rescale = rw_rescale, rw_learn = rw_adapt, ...)
n_fail <- n_fail + out_step$n_fail
## Check nb of failures FC: MODIFY THIS EVENTUALLY
if (n_fail > max_fail)
stop("Number of failures exceeds max_fail: ", n_fail, " failures.")
# Stop rescale
if (rw_rescale && (object$n_iter() >= object$n_rescale())) {
rw_rescale <- FALSE
}
## Store output
if ((i%%thin) == 0) {
ind_sample <- ind_sample + 1
log_marg_like <- object$log_marg_like()
log_prior <- object$log_prior()
if ("p" %in% output)
out$log_marg_like_pen[ind_sample] <- log_marg_like + log_prior
if ("l" %in% output)
out$log_marg_like[ind_sample] <- log_marg_like
if ("a" %in% output)
out$info$accept_rate[ind_sample] <- out_step$accept_rate
if ("s" %in% output)
out$info$rw_step[, ind_sample] <- object$rw_step()
if (return_samples) {
sample_param <- object$sample_param()
sample_latent <- object$sample_latent()
for (var in param_names) {
len <- length(sample_param[[var]])
from <- (ind_sample - 1) * len + 1
to <- (ind_sample - 1) * len + len
out[[var]][from:to] <- sample_param[[var]]
}
if (ind_sample == 1) {
## pre-allocation here to be sure that sample is not empty
for (ivar in 1:length(latent_names)) {
var <- latent_names[[ivar]]
dimen <- dim(sample_latent[[var]])
out[[var]] <- mcmcarray(dim = c(dimen, n_samples),
name = pn_latent$names[[ivar]],
lower = pn_latent$lower[[ivar]],
upper = pn_latent$upper[[ivar]])
}
}
for (var in latent_names) {
len <- length(sample_latent[[var]])
from <- (ind_sample - 1) * len + 1
to <- (ind_sample - 1) * len + len
out[[var]][from:to] <- sample_latent[[var]]
}
}
}
## advance progress bar
rbiips("advance_progress_bar", bar, 1)
}
## number of failures
if ("f" %in% output)
out$info$n_fail <- n_fail
class(out) <- "mcmcarray.list"
return(out)
}
#' @export
biips_pmmh_update <- function(object, ...) UseMethod("biips_pmmh_update")
#' The method \code{biips_pmmh_update} performs adaptation and burn-in
#' iterations for the PMMH algorithm.
#'
#' @export
#' @rdname pmmh-object
#' @param n_iter integer. Number of burn-in iterations.
#' @param thin integer. Thinning interval. Returns samples every \code{thin}
#' iterations (default = 1)
#' @param n_part integer. Number of particles used in SMC algorithms.
#' @param ... Additional arguments to be passed to the SMC
#' algorithm such as \code{rs_thres} and \code{rs_type}.
#' See \code{\link{biips_smc_samples}} for more details.
#' @param output string. Select additional members to be returned in the
#' \code{\link{mcmcarray.list}} output. The string can contain several
#' characters in \code{('p', 'l', 'a', 's', 'f')}. See details. (default = \code{'p'})
#'
#' @details The \code{output} string arguments can be used to query additional
#' members in the \code{\link{mcmcarray.list}} output. If \code{output}
#' contains:
#' \describe{
#' \item{p}{returns member \code{log_marg_like_pen}. \code{\link{mcmcarray}} with penalized
#' log marginal likelihood estimates over iterations.}
#' \item{l}{returns member \code{log_marg_like}. \code{\link{mcmcarray}}
#' with log marginal likelihood estimates over iterations.}
#' \item{a}{returns member \code{info$accept_rate}. \code{\link{mcmcarray}} with acceptance
#' rate over iterations.}
#' \item{s}{returns member \code{info$rw_step}. \code{\link{mcmcarray}}
#' with standard deviations of the random walk over iterations.}
#' \item{f}{returns member \code{info$n_fail}. number of failed SMC algorithms.}
#' }
#'
#' @param max_fail integer. maximum number of failed SMC algorithms allowed.
#' (default=0).
#' @param rw_adapt logical. Activate adaptation of the proposal
#' (default=\code{TRUE}).
#'
#' @return The methods \code{biips_pmmh_update} and \code{biips_pmmh_update}
#' return an object of class \code{\link{mcmcarray.list}}.
#'
#' @return \code{biips_pmmh_samples} output contains one \code{\link{mcmcarray}}
#' member for each monitored variable returned by the \code{param_names()} and
#' \code{latent_names()} member functions of the \code{pmmh} object.
#'
#' @return The members of the \code{\link{mcmcarray.list}} object are
#' \code{\link{mcmcarray}} objects for different variables. Assuming
#' \code{dim} is the dimension of the monitored variable, the
#' \code{\link{mcmcarray}} object is an array of dimension \code{c(dim,
#' n_iter)} with the following attributes (accessible with
#' \code{\link[base]{attr}}):
#' \item{name}{string with the name of the variable.}
#' \item{lower}{vector with the lower bounds of the variable.}
#' \item{upper}{vector with the upper bounds of the variable.}
#'
#' @return If the \code{output} argument is not empty, the output contains
#' additional members. See details.
#'
#' @examples
#' 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')
biips_pmmh_update.pmmh <- function(object, n_iter, n_part, thin = 1, max_fail = 0,
rw_adapt = TRUE, output = "p", ...) {
out <- pmmh_algo(object, n_iter, n_part, thin = thin, return_samples = FALSE,
max_fail = max_fail, rw_adapt = rw_adapt, output = output, ...)
return(invisible(out))
}
#' @export
biips_pmmh_samples <- function(object, ...) UseMethod("biips_pmmh_samples")
#' The method \code{biips_pmmh_samples} performs iterations for the PMMH
#' algorithm and returns samples.
#'
#' @export
#' @rdname pmmh-object
biips_pmmh_samples.pmmh <- function(object, n_iter, n_part, thin = 1, max_fail = 0,
output = "p", ...) {
out <- pmmh_algo(object, n_iter, n_part, thin = thin, return_samples = TRUE,
max_fail = max_fail, output = output, ...)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.