Title: | 'rbi' Helper Functions |
---|---|
Description: | Contains a collection of helper functions to use with 'rbi', the R interface to 'LibBi', described in Murray et al. (2015) <doi:10.18637/jss.v067.i10>. It contains functions to adapt the proposal distribution and number of particles in particle Markov-Chain Monte Carlo, as well as calculating the Deviance Information Criterion (DIC) and converting between times in 'LibBi' results and R time/dates. |
Authors: | Sebastian Funk <[email protected]> |
Maintainer: | Sebastian Funk <[email protected]> |
License: | GPL-3 |
Version: | 0.4.0 |
Built: | 2024-10-31 04:05:29 UTC |
Source: | https://github.com/sbfnk/rbi.helpers |
This function takes the provided libbi
object which has been
run, or a bi file, and returns a the acceptance rate
acceptance_rate(...)
acceptance_rate(...)
... |
parameters to |
acceptance rate
Sebastian Funk
example_run <- rbi::bi_read( system.file(package = "rbi.helpers", "example_run.nc") ) example_model_file <- system.file(package = "rbi", "PZ.bi") example_bi <- rbi::attach_data( rbi::libbi(example_model_file), "output", example_run ) acceptance_rate(example_bi)
example_run <- rbi::bi_read( system.file(package = "rbi.helpers", "example_run.nc") ) example_model_file <- system.file(package = "rbi", "PZ.bi") example_bi <- rbi::attach_data( rbi::libbi(example_model_file), "output", example_run ) acceptance_rate(example_bi)
This function takes the provided libbi
and runs
MCMC at a single point (i.e., repeatedly proposing the same parameters),
adapting the number of particles distribution until the variance of the
log-likelihood crosses the value given as target.variance
(1 by
default).
adapt_particles( x, min = 1, max = 1024, target_variance = 1, quiet = FALSE, target.variance, ... )
adapt_particles( x, min = 1, max = 1024, target_variance = 1, quiet = FALSE, target.variance, ... )
x |
a |
min |
minimum number of particles |
max |
maximum number of particles |
target_variance |
target log-likelihood variance; once this is crossed, the current number of particles will be used |
quiet |
if set to TRUE, will not provide running output of particle numbers tested |
target.variance |
deprecated; use |
... |
parameters for libbi$run |
a libbi
with the desired proposal distribution
example_obs <- rbi::bi_read(system.file(package="rbi", "example_dataset.nc")) example_model <- rbi::bi_model(system.file(package="rbi", "PZ.bi")) example_bi <- rbi::libbi(model = example_model, obs = example_obs) obs_states <- rbi::var_names(example_model, type = "obs") max_time <- max(vapply(example_obs[obs_states], function(x) { max(x[["time"]]) }, 0)) ## Not run: adapted <- adapt_particles(example_bi, nsamples = 128, end_time = max_time) ## End(Not run)
example_obs <- rbi::bi_read(system.file(package="rbi", "example_dataset.nc")) example_model <- rbi::bi_model(system.file(package="rbi", "PZ.bi")) example_bi <- rbi::libbi(model = example_model, obs = example_obs) obs_states <- rbi::var_names(example_model, type = "obs") max_time <- max(vapply(example_obs[obs_states], function(x) { max(x[["time"]]) }, 0)) ## Not run: adapted <- adapt_particles(example_bi, nsamples = 128, end_time = max_time) ## End(Not run)
This function takes the provided libbi
object and
runs MCMC, adapting the proposal distribution until the desired
acceptance rate is achieved. If a scale is given, it will be used
to adapt the proposal at each iteration
adapt_proposal( x, min = 0, max = 1, scale = 2, max_iter = 10, adapt = c("size", "shape", "both"), size = FALSE, correlations = TRUE, truncate = TRUE, quiet = FALSE, ... )
adapt_proposal( x, min = 0, max = 1, scale = 2, max_iter = 10, adapt = c("size", "shape", "both"), size = FALSE, correlations = TRUE, truncate = TRUE, quiet = FALSE, ... )
x |
|
min |
minimum acceptance rate |
max |
maximum acceptance rate |
scale |
scale multiplier/divider for the proposal. If >1 this will be inverted. |
max_iter |
maximum of iterations (default: 10) |
adapt |
what to adapt; if "size" (default), the width of independent proposals will be adapted; if "shape", proposals will be dependent (following a multivariate normal) taking into account empirical correlations; if "both", the size will be adapted before the shape |
size |
(deprecated, use |
correlations |
(deprecated, use |
truncate |
if TRUE, the proposal distributions will be truncated according to the support of the prior distributions |
quiet |
if set to TRUE, will not provide running output of particle numbers tested |
... |
parameters for |
a libbi
with the desired proposal distribution
example_obs <- rbi::bi_read(system.file(package="rbi", "example_dataset.nc")) example_model <- rbi::bi_model(system.file(package="rbi", "PZ.bi")) example_bi <- rbi::libbi(model = example_model, obs = example_obs) obs_states <- rbi::var_names(example_model, type="obs") max_time <- max(vapply(example_obs[obs_states], function(x) { max(x[["time"]]) }, 0)) # adapt to acceptance rate between 0.1 and 0.5 ## Not run: adapted <- adapt_proposal(example_bi, nsamples = 100, end_time = max_time, min = 0.1, max = 0.5, nparticles = 256, correlations = TRUE ) ## End(Not run)
example_obs <- rbi::bi_read(system.file(package="rbi", "example_dataset.nc")) example_model <- rbi::bi_model(system.file(package="rbi", "PZ.bi")) example_bi <- rbi::libbi(model = example_model, obs = example_obs) obs_states <- rbi::var_names(example_model, type="obs") max_time <- max(vapply(example_obs[obs_states], function(x) { max(x[["time"]]) }, 0)) # adapt to acceptance rate between 0.1 and 0.5 ## Not run: adapted <- adapt_proposal(example_bi, nsamples = 100, end_time = max_time, min = 0.1, max = 0.5, nparticles = 256, correlations = TRUE ) ## End(Not run)
Computes the DIC of a libbi object containing Monte-Carlo samples. The effective number of parameters is calculated following Gelman et al., Bayesian Data Analysis: Second Edition, 2004, p. 182.
## S3 method for class 'libbi' DIC(x, bootstrap = 0, ...)
## S3 method for class 'libbi' DIC(x, bootstrap = 0, ...)
x |
a |
bootstrap |
number of bootstrap samples to take, 0 to just take data |
... |
any parameters to be passed to 'bi_read' (e.g., 'burn') |
DIC
Sebastian Funk
example_run <- rbi::bi_read( system.file(package = "rbi", "example_output.nc") ) example_model_file <- system.file(package = "rbi", "PZ.bi") example_bi <- rbi::attach_data( rbi::libbi(example_model_file), "output", example_run ) DIC(example_bi)
example_run <- rbi::bi_read( system.file(package = "rbi", "example_output.nc") ) example_model_file <- system.file(package = "rbi", "PZ.bi") example_bi <- rbi::attach_data( rbi::libbi(example_model_file), "output", example_run ) DIC(example_bi)
This function converts from numeric times (i.e., 0, 1, 2, ...) to actual times or dates
numeric_to_time(x, origin, unit, ...)
numeric_to_time(x, origin, unit, ...)
x |
a |
origin |
the time origin, i.e. the date or time corresponding to time 0 |
unit |
the unit of time that each time step corresponds to; this must be
a unit understood by |
... |
any arguments for |
a list of data frames as returned by bi_read
, but with real
times
This function converts from real times/dates to numeric times (0, 1, 2, ...)
time_to_numeric(x, origin, unit)
time_to_numeric(x, origin, unit)
x |
a data frame containing a "time" column, or a list containing such data frames |
origin |
the time origin, i.e. the date or time corresponding to time 0 |
unit |
the unit of time that each time step corresponds to; this must be
a unit understood by |
a list of data frames that can be passed to libbi