Title: | Flexible Markov Chain Monte Carlo via Reparameterization |
---|---|
Description: | drjacoby is an R package for performing Bayesian inference via Markov chain monte carlo (MCMC). In addition to being highly flexible it implements some advanced techniques that can improve mixing in tricky situations. |
Authors: | Bob Verity [aut, cre], Pete Winskill [aut] |
Maintainer: | Bob Verity <[email protected]> |
License: | MIT + file LICENSE |
Version: | 1.5.4 |
Built: | 2024-10-24 03:15:32 UTC |
Source: | https://github.com/mrc-ide/drjacoby |
Estimate autocorrelation
acf_data(x, lag)
acf_data(x, lag)
x |
Single chain. |
lag |
maximum lag. Must be an integer between 1 and 500. |
Simple function to check that drjacoby package has loaded successfully. Prints "drjacoby loaded successfully!" if so.
check_drjacoby_loaded()
check_drjacoby_loaded()
Create template for cpp
cpp_template(save_as)
cpp_template(save_as)
save_as |
Path of (.cpp) file to create, relative to root of active project. |
Provides a convenient way of defining parameters in the format
required by run_mcmc()
. Each parameter must have the following three
elements, defined in order:
name
- the parameter name.
min
- the minimum value of the parameter. -Inf
is
allowed.
max
- the maximum value of the parameter. Inf
is
allowed.
There following arguments are also optional:
init
- the initial value of the parameter. If running
multiple chains a vector of initial values can be used to specify distinct
values for each chain.
block
- which likelihood block(s) this parameter belongs to.
See vignettes for instructions on using likelihood blocks.
define_params(...)
define_params(...)
... |
a series of named input arguments. |
define_params(name = "mu", min = -10, max = 10, init = 0, name = "sigma", min = 0, max = 5, init = c(1, 2)) define_params(name = "mu1", min = -10, max = 10, init = 0, block = 1, name = "mu2", min = -10, max = 10, init = 0, block = 2, name = "sigma", min = 0, max = 5, init = 1, block = c(1, 2))
define_params(name = "mu", min = -10, max = 10, init = 0, name = "sigma", min = 0, max = 5, init = c(1, 2)) define_params(name = "mu1", min = -10, max = 10, init = 0, block = 1, name = "mu2", min = -10, max = 10, init = 0, block = 2, name = "sigma", min = 0, max = 5, init = 1, block = c(1, 2))
Flexible Markov chain monte carlo via reparameterization using the Jacobean matrix.
_PACKAGE
Estimate sthe Gelman-Rubin (rhat) convergence statistic for a single parameter across multiple chains. Basic method, assuming all chains are of equal length
gelman_rubin(par_matrix, chains, samples)
gelman_rubin(par_matrix, chains, samples)
par_matrix |
Matrix (interations x chains) |
chains |
number of chains |
samples |
number of samples |
Gelman-Rubin statistic
Gelman, A., and D. B. Rubin. 1992. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science 7: 457–511.
https://astrostatistics.psu.edu/RLectures/diagnosticsMCMC.pdf
Plot autocorrelation for specified parameters
plot_autocorrelation(x, lag = 20, par = NULL, chain = 1, phase = "sampling")
plot_autocorrelation(x, lag = 20, par = NULL, chain = 1, phase = "sampling")
x |
an object of class |
lag |
maximum lag. Must be an integer between 1 and 500. |
par |
vector of parameter names. If |
chain |
which chain to plot. |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
Produces a matrix showing the correlation between all parameters from posterior draws.
plot_cor_mat(x, show = NULL, phase = "sampling", param_names = NULL)
plot_cor_mat(x, show = NULL, phase = "sampling", param_names = NULL)
x |
an object of class |
show |
Vector of parameter names to plot. |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
param_names |
Optional vector of names to replace the default parameter names. |
Plots posterior 95% credible intervals over specified set of parameters (defauls to all parameters).
plot_credible(x, show = NULL, phase = "sampling", param_names = NULL)
plot_credible(x, show = NULL, phase = "sampling", param_names = NULL)
x |
an object of class |
show |
vector of parameter names to plot. |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
param_names |
optional vector of names to replace the default parameter names. |
Density plots of all parameters. Use show
and hide
to be more specific about which parameters to plot.
plot_density(x, show = NULL, hide = NULL)
plot_density(x, show = NULL, hide = NULL)
x |
an object of class |
show |
optional vector of parameter names to plot. Parameters matching show will be included. |
hide |
optional vector of parameter names to filter out. Parameters matching hide will be hidden. |
Plot Metropolis coupling acceptance rates between all rungs.
plot_mc_acceptance(x, chain = NULL, phase = "sampling", x_axis_type = 1)
plot_mc_acceptance(x, chain = NULL, phase = "sampling", x_axis_type = 1)
x |
an object of class |
chain |
which chain to plot. If |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
x_axis_type |
how to format the x-axis. 1 = integer rungs, 2 = values of the thermodynamic power. |
Uses ggpairs
function from the GGally
package to
produce scatterplots between all named parameters.
plot_pairs(x, show = NULL, hide = NULL)
plot_pairs(x, show = NULL, hide = NULL)
x |
an object of class |
show |
optional vector of parameter names to plot. Parameters matching show will be included. |
hide |
optional vector of parameter names to filter out. Parameters matching hide will be hidden. |
Plot loglikelihood 95% credible intervals.
plot_rung_loglike( x, chain = 1, phase = "sampling", x_axis_type = 1, y_axis_type = 1 )
plot_rung_loglike( x, chain = 1, phase = "sampling", x_axis_type = 1, y_axis_type = 1 )
x |
an object of class |
chain |
which chain to plot. |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
x_axis_type |
how to format the x-axis. 1 = integer rungs, 2 = values of the thermodynamic power. |
y_axis_type |
how to format the y-axis. 1 = raw values, 2 = truncated at auto-chosen lower limit. 3 = double-log scale. |
Produces scatterplot between two named parameters.
plot_scatter( x, parameter1, parameter2, downsample = TRUE, phase = "sampling", chain = NULL )
plot_scatter( x, parameter1, parameter2, downsample = TRUE, phase = "sampling", chain = NULL )
x |
an object of class |
parameter1 |
name of parameter first parameter. |
parameter2 |
name of parameter second parameter. |
downsample |
whether to downsample output to 200 values max to speed up plotting. |
phase |
which phase to plot. Must be either "burnin" or "sampling". |
chain |
which chain to plot. |
Produce a series of plots corresponding to each parameter, including the raw trace, the posterior histogram and an autocorrelation plot. Plotting objects can be cycled through interactively, or can be returned as an object allowing them to be viewed/edited by the user.
plot_trace( x, show = NULL, hide = NULL, lag = 20, downsample = TRUE, phase = "sampling", chain = NULL, display = TRUE )
plot_trace( x, show = NULL, hide = NULL, lag = 20, downsample = TRUE, phase = "sampling", chain = NULL, display = TRUE )
x |
an object of class |
show |
optional vector of parameter names to plot. Parameters matching show will be included. |
hide |
optional vector of parameter names to filter out. Parameters matching hide will be hidden. |
lag |
maximum lag. Must be an integer between 1 and 500. |
downsample |
boolean. Whether to downsample chain to make plotting more efficient. |
phase |
which phase to plot. Must be either "burnin", "sampling" or "both". |
chain |
which chain to plot. |
display |
boolean. Whether to show plots, if |
Example population growth data.
population_data
population_data
A data frame with 20 rows and 2 variables:
population size
time
...
Run MCMC either with or without parallel tempering turned on.
Minimum inputs include a data object, a data.frame of parameters, a
log-likelihood function and a log-prior function. Produces an object of
class drjacoby_output
, which contains all MCMC output along with
some diagnostics and a record of inputs.
run_mcmc( data, df_params, misc = list(), loglike, logprior, burnin = 1000, samples = 10000, rungs = 1, chains = 5, beta_manual = NULL, alpha = 1, target_acceptance = 0.44, cluster = NULL, coupling_on = TRUE, pb_markdown = FALSE, save_data = TRUE, save_hot_draws = FALSE, silent = FALSE )
run_mcmc( data, df_params, misc = list(), loglike, logprior, burnin = 1000, samples = 10000, rungs = 1, chains = 5, beta_manual = NULL, alpha = 1, target_acceptance = 0.44, cluster = NULL, coupling_on = TRUE, pb_markdown = FALSE, save_data = TRUE, save_hot_draws = FALSE, silent = FALSE )
data |
a named list or data frame or data values. |
df_params |
a data.frame of parameters (see |
misc |
optional list object passed to likelihood and prior. This can be useful for passing values that are not strictly data, for example passing a lookup table to the prior function. |
loglike , logprior
|
the log-likelihood and log-prior functions used in the MCMC. Can either be passed in as R functions (not in quotes), or as character strings naming compiled C++ functions (in quotes). |
burnin |
the number of burn-in iterations. Automatic tuning of proposal standard deviations is only active during the burn-in period. |
samples |
the number of sampling iterations. |
rungs |
the number of temperature rungs used in the parallel tempering
method. By default, |
chains |
the number of independent replicates of the MCMC to run. If a
|
beta_manual |
vector of manually defined |
alpha |
the likelihood for the ith heated chain is
raised to the power |
target_acceptance |
Target acceptance rate. Should be between 0 and 1. Default of 0.44, set as optimum for unvariate proposal distributions. |
cluster |
option to pass in a cluster environment, allowing chains to be run in parallel (see package "parallel"). |
coupling_on |
whether to implement Metropolis-coupling over temperature
rungs. The option of deactivating coupling has been retained for general
interest and debugging purposes only. If this parameter is |
pb_markdown |
whether to run progress bars in markdown mode, meaning they are only updated when they reach 100% to avoid large amounts of output being printed to markdown files. |
save_data |
if |
save_hot_draws |
if |
silent |
whether to suppress all console output. |
Note that both data
and misc
are passed into
log-likelihood/log-prior functions *by reference*. This means if you modify
these objects inside the functions then any changes will persist.
Sample posterior draws from all available chains
sample_chains(x, sample_n, keep_chain_index = FALSE)
sample_chains(x, sample_n, keep_chain_index = FALSE)
x |
an object of class |
sample_n |
An integer number of samples. |
keep_chain_index |
if |
A data.frame of posterior samples