dream_old: Differential Evolution Adaptive Metropolis (DREAM) algorithm
In tpilz/HydroBayes: Bayesian statistics for use in hydrological applications

Description Usage Arguments Details Value Author(s) References

Differential Evolution Adaptive Metropolis (DREAM) algorithm

dream_old(fun, ..., par.info = list(initial = NULL, min = NULL, max = NULL, mu
  = NULL, cov = NULL, val_ini = NULL, bound = NULL, names = NULL), nc, t, d,
  burnin = 0, adapt = 0.1, updateInterval = 10, delta = 3,
  c_val = 0.1, c_star = 1e-12, nCR = 3, p_g = 0.2, beta0 = 1,
  thin = 1, keep_sim = FALSE, checkConvergence = FALSE, verbose = TRUE,
  DEBUG = FALSE)

`fun`	`character`. Name of a function(x, ...) which is evaluated at x being a d-dimensional parameter vector.
`...`	Additional arguments for `fun`.
`par.info`	A `list` characterising the prior sampling.
`nc`	`numeric`. Number of chains evolved in parallel.
`t`	`numeric`. Number of samples from the Markov chain.
`d`	`numeric`. Number of parameters.
`burnin`	`numeric`. Length of the burn-in period as portion of t (`burnin period = burnin * t`). These samples from the Markov chain will not be included in the output. Default: 0.
`adapt`	`numeric`. Length of the adaptation period as portion of t (`adaptation period = adapt * t`). Will be used for the update of crossover probabilities and the replacement of outlier chains. Default: 0.1.
`updateInterval`	`integer`. Interval for crossover probability updates during the adaptation period.
`delta`	`integer`. Maximum number of chain pairs used to generate the jump (default: 3).
`c_val`	`numeric`. Lambda value is sampled from U[-c_val,c_val] (default: 0.1).
`c_star`	`numeric`. Zeta value sampled from N[0,c_star]. Should be small compared to target (i.e. in this case the normal) distribution. Default: 1e-12.
`nCR`	`integer`. Length of vector with crossover probabilities for parameter subspace sampling (default: 3).
`p_g`	`numeric`. Probability for gamma, the jump rate, being equal to 1. Default: 0.2.
`beta0`	`numeric`. Reduce jump distance, e.g. if the average acceptance rate is low (less than 15 %). `0 < beta0 <= 1`. Default: 1 (i.e. jump distance is not adjusted).
`thin`	`integer`. Thinning to be applied to output in case of large `t`. See below.
`keep_sim`	`logical`. Shall simulations generated with `fun` be returned in the output list? If `TRUE`, `fun` needs to return a list with elements 'lp' (the log posterior density) and 'sim' (the simulation time series). Default: FALSE.
`checkConvergence`	`logical`. Shall convergence of the MCMC chain be checked? Currently implemented: Calculating the Gelman-Rubin diagnostic. Takes a lot of time! Default: FALSE.
`verbose`	`logical`. Print progress bar to console? Default: TRUE.
`DEBUG`	`logical`. Option enables further output for error and/or more in-depth analysis. See below. Default: FALSE.
`initial`	`character`. Method for prior sampling. One of: uniform - sampling from a uniform distribution; normal - a (multivariate) normal distribution; latin - latin hypercube sampling; user - value(s) given by the user.
`min`	`numeric`. A d-dimensional vector of minimum values for each parameter to sample from if `initial` is 'uniform' or 'latin'. Also defines bounding region for the proposal, see `bound`.
`max`	`numeric`. A d-dimensional vector of maximum values for each parameter to sample from if `initial` is 'uniform' or 'latin'. Also defines bounding region for the proposal, see `bound`.
`mu`	`numeric`. d-dimensional vector of parameter means if `initial` is 'normal'.
`cov`	`numeric`. d-by-d positive-definite symmetric matrix of parameter covariances if `initial` is 'normal'.
`val_ini`	`numeric` nc-by-d-dimensional matrix of prior values if `initial` is 'user'.
`bound`	`character`. What to do if the proposal parameter is outside the defined min-max limits. One of: bound - proposal is set to min/max value if it is smaller/larger than the defined limit.
`names`	`character` vector of length d with names for the parameters. These can be used within `fun` (in this case, parameter input x of `fun` is a named vector) and will appear in the output list element 'chain'.

To understand the notation (e.g. what is lambda, nCR etc.), have a look at Sect. 3.3 of the reference paper (see below).

list with named elements:

chain: a (1-burnin)*t/thin-by-d-by-nc array of parameter realisations for each iteration and Markov chain;

density: a (1-burnin)*t/thin-by-nc matrix of log-densities computed by pdf at each iteration for each Markov chain;

runtime: time of function execution in seconds;

outlier: a list with adapt*t vectors of outlier indices in nc (value of 0 means no outliers);

AR: a (1-burnin)*t/thin-by-nCR matrix giving the acceptance rate for each sample number and crossover value (first element is NA due to computational reasons);

CR: a (1-burnin)*t/thin-by-nCR matrix giving the selection probability for each sample number and crossover value (first element is NA due to computational reasons).

IF keep_sim == TRUE:

fun_sim: a (1-burnin)*t/thin-by-k-by-nc array of simulation time series (length k) generated with fun coresponding to the parameter realisations of output element chain.

IF checkConvergence == TRUE:

R_stat: a (1-burnin)*t/thin-50+1-by-d matrix giving the Gelman-Rubin convergence diagnostic (note that at least 50 observations are used to compute R_stat).

IF DEBUG == TRUE:

DEBUG: a list with the elements:

J: a t-by-nCR matrix of cumulated Euclidean jump distances during the Markov chain progressing;

dx: a t-by-nc-by-d array of jump proposals;

dx_eff: a t-by-nc-by-d array of accepted jumps;

std: a t-by-d matrix of standard deviations of the chain for each parameter. Note: J_i = J_i-1 + sum( (dx_eff_i/std_i)^2 );

gamma: a t-by-nc matrix of jump rate values;

lambda: a t-by-nc-by-d array of lambda values;

zeta: a t-by-nc-by-d array of zeta values;

jump_diff: a t-by-nc-by-d array of jump differentials ( sum(X_a - X_b) ).

Tobias Pilz tpilz@uni-potsdam.de

Code based on 'Algorithm 5' and 'Algorithm 6' of:

Vrugt, J. A.: "Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts, and MATLAB implementation." Environmental Modelling & Software, 2016, 75, 273 – 316, http://dx.doi.org/10.1016/j.envsoft.2015.08.013.

tpilz/HydroBayes documentation built on May 6, 2019, 3:44 p.m.