Package 'rbi.helpers'

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-08 03:16:46 UTC
Source: https://github.com/sbfnk/rbi.helpers

Help Index


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)

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)

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


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