Sensitivity analysis of parameters with SMC

Source: R/sensitivity.r

Sensitivity analysis of parameters with SMC

biips_smc_sensitivity(object, param_values, n_part, ...)

Arguments

object	biips model object as returned by biips_model.
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.
n_part	integer. Number of particles.
...	Additional arguments to be passed to the SMC algorithm such as rs_thres and rs_type. See biips_smc_samples for more details.

Value

biips_smc_sensitivity returns a list with the following members:

log_marg_like

vector of log marginal likelihood estimates at the different values of the parameters.

log_marg_like_pen

vector of penalized log marginal likelihood estimates at the different values of the parameters.

max_param

parameters value with maximum log_marg_like.

max_log_marg_like

maximum log_marg_like value.

max_param_pen

parameters value with maximum log_marg_like_pen.

max_log_marg_like_pen

maximum log_marg_like_pen value.

See also

biips_model, biips_smc_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
logtau_val <- -10:10
out_sens <- biips_smc_sensitivity(model, list(logtau = logtau_val), n_part)
#> * Analyzing sensitivity with 50 particles
#>   |--------------------------------------------------| 100%
#>   |**************************************************| 21 iterations in 0.05 s

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')