#' Sensitivity analysis of parameters with SMC
#'
#' @export
#' @param object \code{biips} model object as returned by
#' \code{\link{biips_model}}.
#' @param param_values named list of numeric values. Each member is named after
#' a parameter for which we want to study sensitivity. Each member value is a
#' vector or array whose last dimension is equal to the number of values for
#' which we want to study sensitivity.
#' @param n_part integer. Number of particles.
#' @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.
#'
#' @return \code{biips_smc_sensitivity} returns a list with the following
#' members:
#' \item{log_marg_like}{vector of log marginal likelihood estimates
#' at the different values of the parameters.}
#' \item{log_marg_like_pen}{vector of penalized log marginal likelihood
#' estimates at the different values of the parameters.}
#' \item{max_param}{parameters value with maximum \code{log_marg_like}.}
#' \item{max_log_marg_like}{maximum \code{log_marg_like} value.}
#' \item{max_param_pen}{parameters value with maximum \code{log_marg_like_pen}.}
#' \item{max_log_marg_like_pen}{maximum \code{log_marg_like_pen} value.}
#'
#' @seealso \code{\link{biips_model}}, \code{\link{biips_smc_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), logtau = log(1))
#' model <- biips_model(modelfile, data)
#'
#' n_part <- 50
#' logtau_val <- -10:10
#' out_sens <- biips_smc_sensitivity(model, list(logtau = logtau_val), n_part)
#'
#' par(mfrow = c(2, 1))
#' plot(logtau_val, out_sens$log_marg_like, type = 'l', col = 'blue',
#' xlab = 'logtau', ylab = 'log p(y|logtau) ', main = 'SMC sensitivity')
#' points(data$logtau, min(out_sens$log_marg_like), pch = 17, col = 'green')
#'
#' plot(logtau_val, out_sens$log_marg_like_pen, type = 'l', col = 'blue',
#' xlab = 'logtau', ylab = 'log p(y|logtau) + log p(logtau)')
#' plml <- out_sens$log_marg_like_pen
#' ymin <- min(plml[is.finite(plml)])
#' points(data$logtau, ymin, pch = 17, col = 'green')
biips_smc_sensitivity <- function(object, param_values, n_part, ...) {
stopifnot(is.biips(object))
n_param <- length(param_values)
stopifnot(is.list(param_values), n_param > 0)
stopifnot(sapply(param_values, is.numeric))
varnames <- names(param_values)
stopifnot(!is.null(varnames), nchar(varnames) > 0)
pn <- parse_varnames(varnames)
## Get the sizes of parameter values
sizes <- lapply(param_values, function(x) dim(as.array(x)))
## Get the number of parameter values to evaluate (the last dimension)
n_val <- lapply(sizes, function(x) x[length(x)])
n_values <- n_val[[1]]
## Check if number of parameter values match
if (any(n_val != n_values))
stop("Invalid param_values argument: number of values to evaluate does not match")
## Get the dimensions of parameters
dimen <- lapply(sizes, function(x) {
if (length(x) == 1)
1 else x[-length(x)]
})
len <- lapply(dimen, prod)
## stop biips verbosity
verb <- rbiips("verbosity", 0)
on.exit(rbiips("verbosity", verb)) # set verbosity on again when function terminates
## Clone console
model <- clone_model(object)
## initialize -----------
log_marg_like <- rep(NA, times = n_values)
log_marg_like_pen <- rep(NA, times = n_values)
max_log_marg_like <- -Inf
max_log_marg_like_pen <- -Inf
max_param <- NULL
max_param_pen <- NULL
rbiips("message", paste("Analyzing sensitivity with", n_part, "particles"))
## progress bar
bar <- rbiips("progress_bar", n_values, "*", "iterations")
## Iterate over different values --------
for (k in 1:n_values) {
## Update value of each parameter
value <- list()
for (i in 1:n_param) {
var <- varnames[[i]]
ind <- seq(from = (k - 1) * len[[i]] + 1, to = k * len[[i]])
value[[var]] <- param_values[[i]][ind]
dim(value[[var]]) <- dimen[[i]]
## change param value
if (!.Call("change_data", model$ptr(), pn$names[[i]], pn$lower[[i]],
pn$upper[[i]], value[[i]], FALSE, PACKAGE = "rbiips"))
stop("data change failed: invalid parameter. ", var, " = ", paste0(value[i],
collapse = " "))
}
## get log prior
log_prior <- 0
for (i in 1:n_param) {
log_p <- .Call("get_log_prior_density", model$ptr(), pn$names[[i]], pn$lower[[i]],
pn$upper[[i]], PACKAGE = "rbiips")
## FIXME: why?
if (is.na(log_p))
next
if (is.nan(log_p))
stop("Failure evaluating parameter log prior for variable ", varnames[[i]],
" = ", value[i])
log_prior <- log_prior + log_p
}
## run smc sampler
ok <- smc_forward_algo(model, n_part = n_part, ...)
if (!ok)
stop("Failure running SMC forward sampler. parameter: ", paste(names(value),
"=", value, sep = "", collapse = ";"))
## log marginal likelihood
log_marg_like[k] <- rbiips("get_log_norm_const", model$ptr())
if (log_marg_like[k] > max_log_marg_like) {
max_log_marg_like <- log_marg_like[k]
max_param <- value
}
log_marg_like_pen[k] <- log_marg_like[k] + log_prior
if (log_marg_like_pen[k] > max_log_marg_like_pen) {
max_log_marg_like_pen <- log_marg_like_pen[k]
max_param_pen <- value
}
## advance progress bar
rbiips("advance_progress_bar", bar, 1)
}
## output --------
out <- list(log_marg_like = log_marg_like, max_param = max_param, max_log_marg_like = max_log_marg_like,
log_marg_like_pen = log_marg_like_pen, max_param_pen = max_param_pen, max_log_marg_like_pen = max_log_marg_like_pen)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.