Version: 0.4.0
Title: 'rbi' Helper Functions
Author: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Maintainer: Sebastian Funk <sebastian.funk@lshtm.ac.uk>
Imports: rbi (≥ 0.10.0), data.table, lubridate, reshape2, Matrix
Suggests: markdown, testthat, covr (≥ 3.2.0), stringi, knitr, rmarkdown
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.
License: GPL-3
URL: https://libbi.org, https://github.com/sbfnk/rbi, https://github.com/sbfnk/rbi.helpers
BugReports: https://github.com/sbfnk/rbi.helpers/issues
SystemRequirements: libbi (>= 1.4.2)
LazyLoad: no
RoxygenNote: 7.2.3
NeedsCompilation: no
Encoding: UTF-8
VignetteBuilder: knitr
Packaged: 2023-08-19 07:27:44 UTC; eidesfun
Repository: CRAN
Date/Publication: 2023-08-24 11:00:05 UTC

Compute Deviance Information Criterion (DIC) for a libbi model

Description

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.

Usage

## S3 method for class 'libbi'
DIC(x, bootstrap = 0, ...)

Arguments

x

a libbi object

bootstrap

number of bootstrap samples to take, 0 to just take data

...

any parameters to be passed to 'bi_read' (e.g., 'burn')

Value

DIC

Author(s)

Sebastian Funk

Examples

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)

Compute acceptance rate

Description

This function takes the provided libbi object which has been run, or a bi file, and returns a the acceptance rate

Usage

acceptance_rate(...)

Arguments

...

parameters to get_traces (especially 'x')

Value

acceptance rate

Author(s)

Sebastian Funk

Examples

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)

Adapt the number of particles

Description

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

Usage

adapt_particles(
  x,
  min = 1,
  max = 1024,
  target_variance = 1,
  quiet = FALSE,
  target.variance,
  ...
)

Arguments

x

a libbi object

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 target_variance instead

...

parameters for libbi$run

Value

a libbi with the desired proposal distribution

Examples

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)

Adapt the proposal distribution of MCMC using the covariance of samples

Description

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

Usage

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,
  ...
)

Arguments

x

link{libbi} object

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 {adapt} instead) if TRUE (default: FALSE), the size of the (diagonal multivariate normal) proposal distribution will be adapted

correlations

(deprecated, use {adapt} instead) if TRUE (default: FALSE), the shape of the (diagonal multivariate normal) proposal distribution will be adapted according to the empirical covariance

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 sample

Value

a libbi with the desired proposal distribution

Examples

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)

Construct a covariance matrix

Description

This function takes the provided libbi which has been run and returns the square root of the covariance matrix, which can be used for proposal distributions

Usage

get_mvn_params(x, scale = 1, correlations = TRUE, fix)

Arguments

x

a libbi which has been run

scale

a factor by which to scale all elements of the covariance matrix

correlations

logical; if TRUE, correlations are taken into account when constructing the parameters

fix

if this is set, all elements of the covariance matrix will be set to it

Value

the updated bi model

Convert numeric times to actual times or dates

Description

This function converts from numeric times (i.e., 0, 1, 2, ...) to actual times or dates

Usage

numeric_to_time(x, origin, unit, ...)

Arguments

x

a libbi object which has been run, or a list of data frames containing state trajectories (as returned by bi_read)

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 lubridate::period, optionally with a number in advance, e.g. "day" or "2 weeks" or "3 seconds"

...

any arguments for bi_read (e.g., file)

Value

a list of data frames as returned by bi_read, but with real times

Split a unit string

Description

Splits a unit string (e.g., "2 weeks") into the amount (2) and unit ("weeks")

Usage

split_unit(unit_string)

Arguments

unit_string

the string to split

Value

a list with two elements, "num" (the amount) and "unit", , for use with lubridate::period

Author(s)

Sebastian Funk

Convert actual times or dates to numeric times

Description

This function converts from real times/dates to numeric times (0, 1, 2, ...)

Usage

time_to_numeric(x, origin, unit)

Arguments

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 lubridate::period, optionally with a number in advance, e.g. "day" or "2 weeks" or "3 seconds"

Value

a list of data frames that can be passed to libbi

Construct a proposal from run results

Description

This function takes the provided bi_model and adds a generic Gaussian proposal distribution.

Usage

update_proposal(model, truncate = TRUE, blocks = c("parameter", "initial"))

Arguments

model

a bi_model object

truncate

truncate the multivariate normal proposals according to the used priors, e.g. truncating a parameter with beta prior at 0 and 1

blocks

blocks to use (out of "parameter" and "initial")

Value

the updated bi model