dream_test: TEST Implementation of DREAM
In tpilz/HydroBayes: Bayesian statistics for use in hydrological applications

Description Usage Arguments Details Value Author(s) References

TEST Implementation of DREAM

dream_test(prior, pdf, 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, checkConvergence = FALSE,
  ncores = 1, verbose = TRUE, DEBUG = FALSE)

`prior`	A function(N,d) that draws N samples from a d-variate prior distribution. Returns an N-by-d matrix.
`pdf`	A function(prior) that calculates the log-density of the target distribution for given prior. Returns an N-variate vector.
`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.
`checkConvergence`	`logical`. Shall convergence of the MCMC chain be checked? Currently implemented: Calculating the Gelman-Rubin diagnostic. Takes a lot of time! Default: FALSE.
`ncores`	`integer` specifying the number of cores to be used for parallel `pdf` evaluation using the `doMC` package (Linux only!). A value > 1 is only useful for complex `pdf`. Default: 1.
`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.

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

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

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.