diff --git a/BayesianTools/R/MAP.R b/BayesianTools/R/MAP.R index 6877e8c..ff8476d 100644 --- a/BayesianTools/R/MAP.R +++ b/BayesianTools/R/MAP.R @@ -1,8 +1,8 @@ -#' calculates the Maxiumum APosteriori value (MAP) +#' Calculates the Maximum APosteriori value (MAP) #' @author Florian Hartig #' @param bayesianOutput an object of class BayesianOutput (mcmcSampler, smcSampler, or mcmcList) #' @param ... optional values to be passed on the the getSample function -#' @details Currently, this function simply returns the parameter combination with the highest posterior in the chain. A more refined option would be to take the MCMC sample and do additional calculations, e.g. use an optimizer, a kerne delnsity estimator, or some other tool to search / interpolate around the best value in the chain +#' @details Currently, this function simply returns the parameter combination with the highest posterior in the chain. A more refined option would be to take the MCMC sample and do additional calculations, e.g. use an optimizer, a kernel density estimator, or some other tool to search / interpolate around the best value in the chain. #' @seealso \code{\link{WAIC}}, \code{\link{DIC}}, \code{\link{marginalLikelihood}} #' @export MAP <- function(bayesianOutput, ...){ diff --git a/BayesianTools/R/VSEM.R b/BayesianTools/R/VSEM.R index 9dbc09f..4a67b82 100644 --- a/BayesianTools/R/VSEM.R +++ b/BayesianTools/R/VSEM.R @@ -141,7 +141,7 @@ VSEMgetDefaults <- function(){ #' Allows to mix a given parameter vector with a default parameter vector #' @param pars vector with new parameter values -#' @param defaults vector with defaukt parameter values +#' @param defaults vector with default parameter values #' @param locations indices of the new parameter values #' @rdname package-deprecated #' @description This function is deprecated and will be removed by v0.2. @@ -166,7 +166,7 @@ VSEMcreatePAR <- function(days = 1:(3*365)){ #' Create an example dataset, and from that a likelihood or posterior for the VSEM model #' @author Florian Hartig -#' @param likelihoodOnly switch to devide whether to create only a likelihood, or a full bayesianSetup with uniform priors. +#' @param likelihoodOnly switch to decide whether to create only a likelihood, or a full bayesianSetup with uniform priors. #' @param plot switch to decide whether data should be plotted #' @param selection vector containing the indices of the selected parameters #' @details The purpose of this function is to be able to conveniently create a likelihood for the VSEM model for demonstration purposes. The function creates example data --> likelihood --> BayesianSetup, where the latter is the diff --git a/BayesianTools/R/blockUpdate.R b/BayesianTools/R/blockUpdate.R index e63fbc4..9bd796f 100644 --- a/BayesianTools/R/blockUpdate.R +++ b/BayesianTools/R/blockUpdate.R @@ -1,8 +1,7 @@ - #' Determine the groups of correlated parameters #' @author Stefan Paul -#' @param chain MCMC chain including only the parameters (not logP,ll, logP) -#' @param blockSettings list with settings +#' @param chain MCMC chain including only the parameters (not logP, ll, logP) +#' @param blockSettings a list with settings #' @return groups #' @keywords internal updateGroups <- function(chain,blockSettings){ @@ -60,7 +59,7 @@ getBlock <- function(blockSettings){ #' getblockSettings -#' @description Transforms the original settings in settings used in the model runs +#' @description Transforms the original settings to settings used in the model runs #' @param blockUpdate input settings #' @return list with block settings #' @keywords internal diff --git a/BayesianTools/R/classBayesianOutput.R b/BayesianTools/R/classBayesianOutput.R index 2bc278c..51d5011 100644 --- a/BayesianTools/R/classBayesianOutput.R +++ b/BayesianTools/R/classBayesianOutput.R @@ -1,21 +1,20 @@ # NOTE: The functions in this class are just templates that are to be implemented for all subclasses of BayesianOutput. They are not functional. - #' Extracts the sample from a bayesianOutput #' @author Florian Hartig #' @param sampler an object of class mcmcSampler, mcmcSamplerList, smcSampler, smcSamplerList, mcmc, mcmc.list, double, numeric #' @param parametersOnly for a BT output, if F, likelihood, posterior and prior values are also provided in the output -#' @param coda works only for mcmc classes - provides output as a coda object. Note: if mcmcSamplerList contains mcmc samplers such as DE that have several chains, the internal chains will be collapsed. This may not be the desired behavior for all applications. -#' @param start for mcmc samplers start value in the chain. For SMC samplers, start particle +#' @param coda works only for mcmc classes - returns output as a coda object. Note: if mcmcSamplerList contains mcmc samplers such as DE that have several chains, the internal chains will be collapsed. This may not be desired for all applications. +#' @param start for mcmc samplers, start value in the chain. For SMC samplers, start particle #' @param end for mcmc samplers end value in the chain. For SMC samplers, end particle -#' @param thin thinning parameter. Either an integer determining the thinning intervall (default is 1) or "auto" for automatic thinning. -#' @param numSamples sample size (only used if thin = 1). If you want to use numSamples set thin to 1. +#' @param thin thinning parameter. Either an integer determining the thinning interval (default is 1) or "auto" for automatic thinning. +#' @param numSamples sample size (only used if thin = 1). If you want to use numSamples, set thin to 1. #' @param whichParameters possibility to select parameters by index #' @param reportDiagnostics logical, determines whether settings should be included in the output #' @param ... further arguments #' @example /inst/examples/getSampleHelp.R -#' @details If thin is greater than the total number of samples in the sampler object the first and the last element (of each chain if a sampler with multiples chains is used) are sampled. If numSamples is greater than the total number of samples all samples are selected. In both cases a warning is displayed. -#' @details If thin and numSamples is passed, the function will use the thin argument if it is valid and greater than 1, else numSamples will be used. +#' @details If thin is greater than the total number of samples in the sampler object, the first and the last element (of each chain if a sampler with multiples chains is used) are sampled. If numSamples is greater than the total number of samples all samples are selected. A warning will be displayed in both cases. +#' @details If both thin and numSamples are provided, the function will use thin only if it is valid and greater than 1; otherwise, numSamples will be used. #' @export getSample <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics = FALSE, ...) UseMethod("getSample") diff --git a/BayesianTools/R/classBayesianSetup.R b/BayesianTools/R/classBayesianSetup.R index f7dbc8b..0fae73b 100644 --- a/BayesianTools/R/classBayesianSetup.R +++ b/BayesianTools/R/classBayesianSetup.R @@ -9,20 +9,20 @@ #' @param best vector with best prior values #' @param names optional vector with parameter names #' @param parallel parallelization option. Default is F. Other options include T, or "external". See details. -#' @param parallelOptions list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global environment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details and Examples). -#' @param catchDuplicates Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns. +#' @param parallelOptions list containing three lists.\itemize{ \item First, "packages" determines the R packages necessary to run the likelihood function.\item Second, "variables" - the objects in the global environment needed to run the likelihood function and \item Third, "dlls" is needed to run the likelihood function (see Details and Examples). } +#' @param catchDuplicates logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns. #' @param plotLower vector with lower limits for plotting #' @param plotUpper vector with upper limits for plotting #' @param plotBest vector with best values for plotting #' @details If prior is of class prior (e.g. create with \code{\link{createPrior}}), priorSampler, lower, upper and best will be ignored.\cr If prior is a function (log prior density), priorSampler (custom sampler), or lower/upper (uniform sampler) is required.\cr If prior is NULL, and lower and upper are passed, a uniform prior (see \code{\link{createUniformPrior}}) will be created with boundaries lower and upper. #' -#' For parallelization, Bayesiantools requies that the likelihood can evaluate several parameter vectors (supplied as a matrix) in parallel. +#' For parallelization, Bayesiantools requires that the likelihood can evaluate multiple parameter vectors (supplied as a matrix) in parallel. #' -#' * parallel = T means that an automatic parallelization of the likelihood via a standard R socket cluster is attempted, using the function \code{\link{generateParallelExecuter}}. By default, of the N cores detected on the computer, N-1 cores are requested. Alternatively, you can provide a integer number to parallel, specifying the cores reserved for the cluster. When the cluster is cluster is created, a copy of your workspace, including DLLs and objects are exported to the cluster workers. Because this can be very inefficient, you can explicitly specify the packages, objects and DLLs that are to be exported via parallelOptions. Using parallel = T requires that the function to be parallelized is well encapsulate, i.e. can run on a shared memory / shared hard disk machine in parallel without interfering with each other. +#' * parallel = T attempts to parallelize likelihood via a standard R socket cluster using the \code{\link{generateParallelExecuter}} function. By default, of the N cores detected on the computer, N-1 cores are requested. Alternatively, you can provide a integer number to parallel, specifying the cores reserved for the cluster. When the cluster is created, a copy of your workspace, including DLLs and objects are exported to the cluster workers. As this approach can be highly inefficient, it is recommended to explicitly specify the packages, objects and DLLs to export using parallelOptions. Using parallel = T requires that the function to be parallelized is well encapsulated, i.e. can run in parallel on a shared memory / shared hard disk machine in parallel without interfering with each other. #' -#' If automatic parallelization cannot be done (e.g. because dlls are not thread-safe or write to shared disk), and only in this case, you should specify parallel = "external". In this case, it is assumed that the likelihood is programmed such that it accepts a matrix with parameters as columns and the different model runs as rows. It is then up to the user if and how to parallelize this function. This option gives most flexibility to the user, in particular for complicated parallel architecture or shared memory problems. +#' If automatic parallelization is not possible (e.g., because dlls are not thread-safe or write to shared disk), and only in this case, you should specify parallel = "external". In this case, it is assumed that the likelihood is programmed to accept a matrix with parameters as columns and the different model runs as rows. The user can then choose whether and how to parallelize this function. This option provides optimal flexibility for the user, especially regarding complicated parallel architectures or shared memory issues. #' -#' For more details on parallelization, make sure to read both vignettes, in particular the section on the likelihood in the main vignette, and the section on parallelization in the vignette on interfacing models. +#' For more details on parallelization, make sure to read both vignettes, especially the section on likelihood in the main vignette and the section on parallelization in the vignette on interfacing models. #' #' @export #' @seealso \code{\link{checkBayesianSetup}} \cr @@ -168,7 +168,7 @@ checkBayesianSetup <- function(bayesianSetup, parallel = F){ #' Function to close cluster in BayesianSetup #' @author Stefan Paul #' @description Function closes -#' the parallel executer (if available) +#' the parallel executor (if available) #' @param bayesianSetup object of class BayesianSetup #' @export stopParallel <- function(bayesianSetup){ diff --git a/BayesianTools/R/classLikelihood.R b/BayesianTools/R/classLikelihood.R index 905d9a5..b6bdfb6 100644 --- a/BayesianTools/R/classLikelihood.R +++ b/BayesianTools/R/classLikelihood.R @@ -1,13 +1,14 @@ #' Creates a standardized likelihood class#' #' @author Florian Hartig -#' @param likelihood Log likelihood density -#' @param names Parameter names (optional) -#' @param parallel parallelization , either i) no parallelization --> F, ii) native R parallelization --> T / "auto" will select n-1 of your available cores, or provide a number for how many cores to use, or iii) external parallelization --> "external". External means that the likelihood is already able to execute parallel runs in form of a matrix with -#' @param catchDuplicates Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns. -#' @param parallelOptions list containing two lists. First "packages" determines the R packages necessary to run the likelihood function. Second "objects" the objects in the global envirnment needed to run the likelihood function (for details see \code{\link{createBayesianSetup}}). +#' @param likelihood log likelihood density +#' @param names parameter names (optional) +#' @param parallel parallelization , either i) no parallelization --> F, ii) native R parallelization --> T / "auto" will select n-1 of your available cores, or provide a number for how many cores to use, or iii) external parallelization --> "external". External means that the likelihood is already able to execute parallel runs in the form of a matrix. +#' @param catchDuplicates logical, determines whether unique parameter combinations should only be evaluated once. This is only applicable when the likelihood accepts a matrix with parameters as columns. +#' @param parallelOptions a list containing two lists. First, "packages" specifies the R packages necessary to run the likelihood function. Second, "objects" contains the objects in the global environment needed to run the likelihood function (for details see \code{\link{createBayesianSetup}}). #' @param sampler sampler #' @seealso \code{\link{likelihoodIidNormal}} \cr #' \code{\link{likelihoodAR1}} \cr +#' @example /inst/examples/createLikelihoodHelp.R #' @export createLikelihood <- function(likelihood, names = NULL, parallel = F, catchDuplicates=T, sampler = NULL, parallelOptions = NULL){ @@ -105,7 +106,6 @@ createLikelihood <- function(likelihood, names = NULL, parallel = F, catchDuplic #library(mvtnorm) #library(sparseMVN) - #' Normal / Gaussian Likelihood function #' @author Florian Hartig #' @param predicted vector of predicted values @@ -137,7 +137,7 @@ likelihoodAR1 <- function(predicted, observed, sd, a){ res = predicted - observed - # this calculates the unconditiona LL for this data, see e.g. http://stat.unicas.it/downloadStatUnicas/seminari/2008/Julliard0708_1.pdf + # this calculates the unconditional LL for this data, see e.g. http://stat.unicas.it/downloadStatUnicas/seminari/2008/Julliard0708_1.pdf ll = 0.5 * ( - n * log(2*pi) - n * log(sd^2) diff --git a/BayesianTools/R/classMcmcSamplerList.R b/BayesianTools/R/classMcmcSamplerList.R index 67fad61..d894696 100644 --- a/BayesianTools/R/classMcmcSamplerList.R +++ b/BayesianTools/R/classMcmcSamplerList.R @@ -1,7 +1,7 @@ #' Convenience function to create an object of class mcmcSamplerList from a list of mcmc samplers #' @author Florian Hartig -#' @param mcmcList a list with each object being an mcmcSampler -#' @return Object of class "mcmcSamplerList" +#' @param mcmcList list of objects, each of which is an mcmcSampler +#' @return object of class "mcmcSamplerList" #' @export createMcmcSamplerList <- function(mcmcList){ # mcmcList <- list(mcmcList) -> This line didn't make any sense at all. Better would be to allow the user to simply provide several inputs without a list, but I guess the list option should be maintained, as this is convenient when scripting. @@ -12,6 +12,7 @@ createMcmcSamplerList <- function(mcmcList){ return(mcmcList) } + #' @author Stefan Paul #' @method summary mcmcSamplerList #' @export diff --git a/BayesianTools/R/classPosterior.R b/BayesianTools/R/classPosterior.R index 9c4b3c5..e083b3c 100644 --- a/BayesianTools/R/classPosterior.R +++ b/BayesianTools/R/classPosterior.R @@ -1,8 +1,8 @@ #' Creates a standardized posterior class #' @author Florian Hartig #' @param prior prior class -#' @param likelihood Log likelihood density -#' @details Function is internally used in \code{\link{createBayesianSetup}} to create a standarized posterior class. +#' @param likelihood log likelihood density +#' @details Function is internally used in \code{\link{createBayesianSetup}} to create a standardized posterior class. #' @export createPosterior <- function(prior, likelihood){ @@ -53,5 +53,4 @@ createPosterior <- function(prior, likelihood){ # x$density(c(0.2,0.2)) # prior$density(c(2,2)) # -# # x = c(0.2,0.2) diff --git a/BayesianTools/R/classPrior.R b/BayesianTools/R/classPrior.R index efb0215..57ca064 100644 --- a/BayesianTools/R/classPrior.R +++ b/BayesianTools/R/classPrior.R @@ -1,11 +1,11 @@ #' Creates a standardized prior class #' @author Florian Hartig -#' @param density Prior density +#' @param density prior density #' @param sampler Sampling function for density (optional) #' @param lower vector with lower bounds of parameters #' @param upper vector with upper bounds of parameter #' @param best vector with "best" parameter values -#' @details This is the general prior generator. It is highly recommended to not only implement the density, but also the sampler function. If this is not done, the user will have to provide explicit starting values for many of the MCMC samplers. Note the existing, more specialized prior function. If your prior can be created by those, they are preferred. Note also that priors can be created from an existing MCMC output from BT, or another MCMC sample, via \code{\link{createPriorDensity}}. +#' @details This is the general prior generator. It is highly recommended to implement both the density and sampler function. If not, the user will have to provide explicit starting values for many of the MCMC samplers. Note the existing, more specialized prior functions. It is recommended to use those specialized prior functions, if possible. Also note that priors can be created from an existing MCMC output from BT, or another MCMC sample, via \code{\link{createPriorDensity}}. #' @note min and max truncate, but not re-normalize the prior density (so, if a pdf that integrated to one is truncated, the integral will in general be smaller than one). For MCMC sampling, this doesn't make a difference, but if absolute values of the prior density are a concern, one should provide a truncated density function for the prior. #' @export #' @seealso \code{\link{createPriorDensity}} \cr @@ -247,7 +247,6 @@ createPriorDensity <- function(sampler, method = "multivariate", eps = 1e-10, lo #' @author Maximilian Pichler - #' @export print.prior <- function(x, ...){ diff --git a/BayesianTools/R/classSMCSamplerList.R b/BayesianTools/R/classSMCSamplerList.R index 529e0e0..d4e462b 100644 --- a/BayesianTools/R/classSMCSamplerList.R +++ b/BayesianTools/R/classSMCSamplerList.R @@ -1,7 +1,7 @@ #' Convenience function to create an object of class SMCSamplerList from a list of mcmc samplers #' @author Florian Hartig #' @param ... a list of MCMC samplers -#' @return a list of class smcSamplerList with each object being an smcSampler +#' @return a list of class smcSamplerList with objects of smcSampler #' @export createSmcSamplerList <- function(...){ smcList <- list(...) @@ -35,6 +35,7 @@ plot.smcSamplerList <- function(x, ...){ marginalPlot(x, ...) } + #' @export getSample.smcSamplerList <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics = FALSE, ...){ diff --git a/BayesianTools/R/classSmcSampler.R b/BayesianTools/R/classSmcSampler.R index 2eb931a..051d2b1 100644 --- a/BayesianTools/R/classSmcSampler.R +++ b/BayesianTools/R/classSmcSampler.R @@ -52,6 +52,7 @@ summary.smcSampler<- function(object, ...){ summary(getSample(sampler, ...)) } + #' @method plot smcSampler #' @export plot.smcSampler<- function(x, ...){ diff --git a/BayesianTools/R/codaFunctions.R b/BayesianTools/R/codaFunctions.R index 715d27d..e97afb9 100644 --- a/BayesianTools/R/codaFunctions.R +++ b/BayesianTools/R/codaFunctions.R @@ -1,7 +1,7 @@ #' Function to combine chains #' #' @param x a list of MCMC chains -#' @param merge logical determines whether chains should be merged +#' @param merge logical, should chains be merged? (T or F) #' @return combined chains #' #' @note to combine several chains to a single McmcSamplerList, see \code{\link{createMcmcSamplerList}} @@ -44,10 +44,10 @@ combineChains <- function(x, merge = T){ #' Helper function to change an object to a coda mcmc class, #' #' @param chain mcmc Chain -#' @param start for mcmc samplers start value in the chain. For SMC samplers, start particle -#' @param end for mcmc samplers end value in the chain. For SMC samplers, end particle +#' @param start For MCMC samplers, the initial value in the chain. For SMC samplers, initial particle +#' @param end For MCMC samplers, the end value in the chain. For SMC samplers, end particle. #' @param thin thinning parameter -#' @return object of class coda::mcmc +#' @return object an object of class coda::mcmc #' @details Very similar to coda::mcmc but with less overhead #' @keywords internal makeObjectClassCodaMCMC <- function (chain, start = 1, end = numeric(0), thin = 1){ diff --git a/BayesianTools/R/convertCoda.R b/BayesianTools/R/convertCoda.R index b00739a..c269ac0 100644 --- a/BayesianTools/R/convertCoda.R +++ b/BayesianTools/R/convertCoda.R @@ -1,16 +1,12 @@ - #' Convert coda::mcmc objects to BayesianTools::mcmcSampler -#' @description Function is used to make the plot and diagnostic functions -#' available for coda::mcmc objects -#' @param sampler An object of class mcmc or mcmc.list -#' @param names vector giving the parameter names (optional) -#' @param info matrix (or list with matrices for mcmc.list objects) with three coloumns containing log posterior, log likelihood and log prior of the sampler for each time step (optional; but see Details) -#' @param likelihood likelihood function used in the sampling (see Details) -#' @details The parameter 'likelihood' is optional for most functions but can be needed e.g for -#' using the \code{\link{DIC}} function. +#' @description Function to support plotting and diagnostic functions for coda::mcmc objects. +#' @param sampler an object of class mcmc or mcmc.list +#' @param names a vector with parameter names (optional) +#' @param info a matrix (or list with matrices for mcmc.list objects) with three columns containing log posterior, log likelihood and log prior of the sampler for each time step (optional; but see Details) +#' @param likelihood likelihood function used for sampling (see Details) +#' @details The parameter 'likelihood' is optional for most functions but can be needed e.g for \code{\link{DIC}} function. #' -#' Also the parameter info is optional for most uses. However for some functions (e.g. \code{\link{MAP}}) -#' the matrix or single coloumns (e.g. log posterior) are necessary for the diagnostics. +#' Also, the parameter information is typically optional for most uses. However, for certain functions (e.g. \code{\link{MAP}}), the matrix or single columns (e.g. log posterior) are necessary for diagnostics. #' @export convertCoda <- function(sampler, names = NULL, info = NULL, likelihood = NULL){ diff --git a/BayesianTools/R/getVolume.R b/BayesianTools/R/getVolume.R index 3eef157..d9a7ede 100644 --- a/BayesianTools/R/getVolume.R +++ b/BayesianTools/R/getVolume.R @@ -1,10 +1,10 @@ #' Calculate posterior volume #' @author Florian Hartig -#' @param sampler an object of superclass bayesianOutput or any other class that has the getSample function implemented (e.g. Matrix) -#' @param prior schould also prior volume be calculated +#' @param sampler an object of superclass bayesianOutput or any other class that has implemented the getSample function (e.g. Matrix) +#' @param prior logical, should prior volume be calculated? #' @param method method for volume estimation. Currently, the only option is "MVN" #' @param ... additional parameters to pass on to the \code{\link{getSample}} -#' @details The idea of this function is to provide an estimate of the "posterior volume", i.e. how "broad" the posterior is. One potential application is to the overall reduction of parametric uncertainty between different data types, or between prior and posterior. +#' @details The idea of this function is to provide an estimate of the "posterior volume", i.e. how "broad" the posterior is. One potential application is the overall reduction of parametric uncertainty between different data types, or between prior and posterior. #' #' Implemented methods for volume estimation: #' @@ -30,6 +30,4 @@ getVolume <- function(sampler, prior = F, method = "MVN", ...){ }else stop("BayesianTools: unknown method argument in getVolume") return(list(priorVol = priorVol, postVol = postVol)) }else return(postVol) -} - - +} \ No newline at end of file diff --git a/BayesianTools/R/marginalLikelihood.R b/BayesianTools/R/marginalLikelihood.R index 764c4e3..15be3da 100644 --- a/BayesianTools/R/marginalLikelihood.R +++ b/BayesianTools/R/marginalLikelihood.R @@ -1,9 +1,7 @@ - # Motivation for this functions from # https://radfordneal.wordpress.com/2008/08/17/the-harmonic-mean-of-the-likelihood-worst-monte-carlo-method-ever/ # https://gist.github.com/gaberoo/4619102 - # ' @export #marginalLikelihood <- function(x,lik,V,sampler$setup$likelihood$density,sampler$setup$prior$density,..., num.samples=1000,log=TRUE) UseMethod("marginalLikelihood") @@ -24,17 +22,17 @@ #' #' In BT, we return the log ML, so you will have to exp all values for this formula. #' -#' It is well-known that the ML is VERY dependent on the prior, and in particular the choice of the width of uninformative priors may have major impacts on the relative weights of the models. It has therefore been suggested to not use the ML for model averaging / selection on uninformative priors. If you have no informative priors, and option is to split the data into two parts, use one part to generate informative priors for the model, and the second part for the model selection. See help for an example. +#' It is well-known that the ML is strongly dependent on the prior, and in particular the choice of the width of uninformative priors may have major impacts on the relative weights of the models. It has therefore been suggested to not use the ML for model averaging / selection on uninformative priors. If you have no informative priors, and option is to split the data into two parts, use one part to generate informative priors for the model, and the second part for the model selection. See help for an example. #' #' The marginalLikelihood function currently implements four ways to calculate the marginal likelihood. Be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that any of the implemented algorithms will converge reasonably fast. The recommended (and default) method is the method "Chib" (Chib and Jeliazkov, 2001), which is based on MCMC samples, with a limited number of additional calculations. Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions. #' -#' The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unrealiable and usually should not be used. +#' The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unreliable and usually should not be used. #' #' The third method is simply sampling from the prior. While in principle unbiased, it will only converge for a large number of samples, and is therefore numerically inefficient. #' #' The Bridge method uses bridge sampling as implemented in the R package "bridgesampling". It is potentially more exact than the Chib method, but might require more computation time. However, this may be very dependent on the sampler. #' -#' @return A list with log of the marginal likelihood, as well as other diagnostics depending on the chose method +#' @return A list with log of the marginal likelihood, as well as other diagnostics depending on the chosen method #' #' @example /inst/examples/marginalLikelihoodHelp.R #' @references diff --git a/BayesianTools/R/mcmcDE.R b/BayesianTools/R/mcmcDE.R index 734f080..a338947 100644 --- a/BayesianTools/R/mcmcDE.R +++ b/BayesianTools/R/mcmcDE.R @@ -1,17 +1,15 @@ - - #' Differential-Evolution MCMC #' @author Francesco Minunno and Stefan Paul #' @param bayesianSetup a BayesianSetup with the posterior density function to be sampled from #' @param settings list with parameter settings -#' @param startValue (optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population. +#' @param startValue (optional) either a matrix with start population, a number defining the number of chains to be run or a function that samples a starting population. #' @param iterations number of function evaluations. #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. #' @param thin thinning parameter. Determines the interval in which values are recorded. #' @param f scaling factor gamma #' @param eps small number to avoid singularity #' @param blockUpdate list determining whether parameters should be updated in blocks. For possible settings see Details. -#' @param message logical determines whether the sampler's progress should be printed +#' @param message logical, specifies whether to print the progress of the sampler. #' @references Braak, Cajo JF Ter. "A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces." Statistics and Computing 16.3 (2006): 239-249. #' @export #' @example /inst/examples/DEfamilyHelp.R @@ -24,19 +22,7 @@ #' \item{"random"} { random blocking. Using k (see below)} #' \item{"user"} { user defined groups. Using groups (see below)} #' } -#' Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength -#' of the correlation used to group parameter and "groups" is used for user defined groups. -#' "groups" is a vector containing the group number for each parameter. E.g. for three parameters -#' with the first two in one group, "groups" would be c(1,1,2). -#' Further pSel and pGroup can be used to influence the choice of groups. In the sampling process -#' a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities -#' for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. -#' For updating one or two groups with the same probability pSel = c(1,1). By default all numbers -#' have the same probability. -#' The same principle is used in pGroup. Here the user can influence the probability of each group -#' to be updated. By default all groups have the same probability. -#' Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall -#' in which the groups are evaluated. +#' Further, seven parameters can be specified. "k" defines the number of groups, "h" the strength of the correlation used to group the parameters and "groups" is used for user defined groups. "groups" is a vector containing the group number for each parameter. E.g. for three parameters with the first two in one group, "groups" would be c(1,1,2). Moreover, pSel and pGroup can be used to influence the choice of groups. In the sampling process a number of groups are drawn at random and updated. pSel is a vector containing relative probabilities for updating the respective number of groups. E.g. To update one group at a time pSel = 1. For updating one or two groups with the same probability pSel = c(1,1). By default all numbers have the same probability. The same principle is used for pGroup. Here the user can influence the probability of each group to be updated. By default all groups have the same probability. Finally, "groupStart" defines the starting point of the groupUpdate and "groupIntervall" - the interval in which the groups are evaluated. DE <- function(bayesianSetup, settings = list( diff --git a/BayesianTools/R/mcmcDEzs.R b/BayesianTools/R/mcmcDEzs.R index ed99255..f774508 100644 --- a/BayesianTools/R/mcmcDEzs.R +++ b/BayesianTools/R/mcmcDEzs.R @@ -4,9 +4,9 @@ #' @author Francesco Minunno and Stefan Paul #' @param bayesianSetup a BayesianSetup with the posterior density function to be sampled from #' @param settings list with parameter settings -#' @param startValue (optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population. +#' @param startValue (optional) either a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population. #' @param Z starting Z population -#' @param iterations iterations to run +#' @param iterations number of iterations to run #' @param pSnooker probability of Snooker update #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. #' @param thin thinning parameter. Determines the interval in which values are recorded. @@ -17,12 +17,12 @@ #' @param eps.mult random term (multiplicative error) #' @param eps.add random term #' @param blockUpdate list determining whether parameters should be updated in blocks. For possible settings see Details. -#' @param message logical determines whether the sampler's progress should be printed +#' @param message logical, specifies whether to print the progress of the sampler. #' @references ter Braak C. J. F., and Vrugt J. A. (2008). Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing http://dx.doi.org/10.1007/s11222-008-9104-9 #' @export #' @example /inst/examples/DEfamilyHelp.R #' @seealso \code{\link{DE}} -#' @details For parallel computing, the likelihood density in the bayesianSetup needs to be parallelized, i.e. needs to be able to operate on a matrix of proposals +#' @details For parallel computing, the likelihood density in the bayesianSetup needs to be parallelized, i.e., it needs to be able to operate on a matrix of proposals #' #' For blockUpdate the first element in the list determines the type of blocking. #' Possible choices are @@ -32,18 +32,18 @@ #' \item{"random"} { random blocking. Using k (see below)} #' \item{"user"} { user defined groups. Using groups (see below)} #' } -#' Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength +#' Further, seven parameters can be specified. "k" defines the number of groups, "h" the strength #' of the correlation used to group parameter and "groups" is used for user defined groups. #' "groups" is a vector containing the group number for each parameter. E.g. for three parameters #' with the first two in one group, "groups" would be c(1,1,2). -#' Further pSel and pGroup can be used to influence the choice of groups. In the sampling process -#' a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities -#' for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. +#' Moreover, pSel and pGroup can be used to influence the choice of groups. In the sampling process +#' a number of groups is drawn at random and updated. pSel is a vector containing relative probabilities +#' for updating the respective number of groups. E.g. To update one group at a time pSel = 1. #' For updating one or two groups with the same probability pSel = c(1,1). By default all numbers #' have the same probability. -#' The same principle is used in pGroup. Here the user can influence the probability of each group +#' The same principle is used in pGroup. Here, the user can influence the probability of each group #' to be updated. By default all groups have the same probability. -#' Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall +#' Finally, "groupStart" defines the starting point of the groupUpdate and "groupIntervall" - the interval #' in which the groups are evaluated. DEzs <- function(bayesianSetup, settings = list(iterations=10000, diff --git a/BayesianTools/R/mcmcDREAM.R b/BayesianTools/R/mcmcDREAM.R index ead86c5..969222e 100644 --- a/BayesianTools/R/mcmcDREAM.R +++ b/BayesianTools/R/mcmcDREAM.R @@ -2,43 +2,42 @@ #' DREAM #' @author Stefan Paul -#' @param bayesianSetup Object of class 'bayesianSetup' or 'bayesianOuput'. +#' @param bayesianSetup object of class 'bayesianSetup' or 'bayesianOuput'. #' @param settings list with parameter values -#' @param iterations Number of model evaluations +#' @param iterations number of model evaluations #' @param nCR parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly. -#' @param updateInterval determining the intervall for the pCR update +#' @param updateInterval determines the interval for the pCR update #' @param gamma Kurtosis parameter Bayesian Inference Scheme #' @param eps Ergodicity term #' @param e Ergodicity term -#' @param pCRupdate If T, crossover probabilities will be updated +#' @param pCRupdate logical, if T, crossover probabilities will be updated #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. -#' @param thin thin thinning parameter. Determines the interval in which values are recorded. -#' @param adaptation Number or percentage of samples that are used for the adaptation in DREAM (see Details). -#' @param DEpairs Number of pairs used to generate proposal -#' @param startValue eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior. -#' @param consoleUpdates Intervall in which the sampling progress is printed to the console -#' @param message logical determines whether the sampler's progress should be printed +#' @param thin thinning parameter. Determines the interval in which values are recorded. +#' @param adaptation number or percentage of samples that are used for the adaptation in DREAM (see Details). +#' @param DEpairs number of pairs used to generate proposal +#' @param startValue either a matrix containing the start values (see details), an integer to define the number of chains to be run, a function to sample the start values or NUll - in which case the values are sampled from the prior. +#' @param consoleUpdates interval at which the sampling progress is printed to the console +#' @param message logical, determines whether the sampler's progress should be printed #' @return mcmc.object containing the following elements: chains, X, pCR #' @references Vrugt, Jasper A., et al. "Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling." International Journal of Nonlinear Sciences and Numerical Simulation 10.3 (2009): 273-290. -#' @details Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler -#' from the last iteration. Due to the sampler's internal structure you can only use the output -#' of DREAM. -#' If you provide a matrix with start values the number of rows determines the number of chains that are run. +#' @details Instead of a bayesianSetup, the function can take the output of a previous run to restart the sampler +#' from the last iteration. Due to the sampler's internal structure you can only use the output of DREAM. +#' If you provide a matrix with start values, the number of rows determines the number of chains that will be run. #' The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr #' There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly -#' the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +#' The algorithm implemented here does not have an automatic stopping criterion. Hence, it will #' always run the number of iterations specified by the user. Also, convergence is not #' monitored and left to the user. This can easily be done with coda::gelman.diag(chain). -#' Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here.\cr\cr +#' Furthermore, the delayed rejection step proposed in Vrugt et al. (2009) is not implemented here.\cr\cr #' #' During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. -#' First the disribution of crossover values is tuned to favor large jumps in the parameter space. +#' First, the disribution of crossover values is tuned to favor large jumps in the parameter space. #' The crossover probabilities determine how many parameters are updated simultaneously. -#' Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +#' Second, outlier chains are replaced as they can largely deteriorate the sampler's performance. #' However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain #' should be discarded when summarizing posterior moments. This can be done automatically during the -#' sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between -#' the burnin and adaptation phase to allow the user more flexibility in the sampler's settings. +#' sampling process (i.e. burn-in > adaptation) or subsequently by the user. We chose to distinguish between +#' the burn-in and adaptation phase to allow the user more flexibility in the sampler's settings. #' #' #' @example /inst/examples/DEfamilyHelp.R diff --git a/BayesianTools/R/mcmcDREAM_helperFunctions.R b/BayesianTools/R/mcmcDREAM_helperFunctions.R index a798bef..4781ccc 100644 --- a/BayesianTools/R/mcmcDREAM_helperFunctions.R +++ b/BayesianTools/R/mcmcDREAM_helperFunctions.R @@ -1,8 +1,6 @@ - - ##' Generates matrix of CR values based on pCR -##' @param pCR Vector of crossover probabilities. Needs to be of length nCR. -##' @param settings settings list +##' @param pCR vector of crossover probabilities. Needs to be of length nCR. +##' @param settings list of settings ##' @param Npop number of chains ##' @return Matrix with CR values #' @keywords internal @@ -34,11 +32,9 @@ generateCRvalues <- function(pCR,settings, Npop){ } - - #' Adapts pCR values -#' @param CR Vector of crossover probabilities. Needs to be of length nCR. -#' @param settings settings list +#' @param CR vector of crossover probabilities. Needs to be of length nCR. +#' @param settings list of settings #' @param delta vector with differences #' @param lCR values to weight delta #' @param Npop number of chains. diff --git a/BayesianTools/R/mcmcDREAMzs.R b/BayesianTools/R/mcmcDREAMzs.R index 0467da3..634d1e2 100644 --- a/BayesianTools/R/mcmcDREAMzs.R +++ b/BayesianTools/R/mcmcDREAMzs.R @@ -2,46 +2,46 @@ #' DREAMzs #' @author Stefan Paul -#' @param bayesianSetup Object of class 'bayesianSetup' or 'bayesianOuput'. +#' @param bayesianSetup object of class 'bayesianSetup' or 'bayesianOuput'. #' @param settings list with parameter values -#' @param iterations Number of model evaluations -#' @param nCR parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly. -#' @param updateInterval determining the intervall for the pCR (crossover probabilities) update -#' @param gamma Kurtosis parameter Bayesian Inference Scheme. +#' @param iterations number of model evaluations +#' @param nCR parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly. +#' @param updateInterval determines the interval for the pCR (crossover probabilities) update +#' @param gamma kurtosis parameter for Bayesian Inference Scheme. #' @param eps Ergodicity term #' @param e Ergodicity term -#' @param pCRupdate Update of crossover probabilities +#' @param pCRupdate update of crossover probabilities #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. -#' @param thin thin thinning parameter. Determines the interval in which values are recorded. -#' @param adaptation Number or percentage of samples that are used for the adaptation in DREAM (see Details) -#' @param DEpairs Number of pairs used to generate proposal +#' @param thin thinning parameter. Determines the interval in which values are recorded. +#' @param adaptation number or percentage of samples that are used for the adaptation in DREAM (see Details) +#' @param DEpairs number of pairs used to generate proposal #' @param ZupdateFrequency frequency to update Z matrix #' @param pSnooker probability of snooker update #' @param Z starting matrix for Z -#' @param startValue eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior. -#' @param consoleUpdates Intervall in which the sampling progress is printed to the console -#' @param message logical determines whether the sampler's progress should be printed +#' @param startValue either a matrix containing the start values (see details), an integer to define the number of chains to be run, a function to sample the start values or NUll - in which case the values are sampled from the prior. +#' @param consoleUpdates interval in which the sampling progress is printed to the console +#' @param message logical, determines whether the sampler's progress should be printed #' @return mcmc.object containing the following elements: chains, X, pCR, Z #' @references Vrugt, Jasper A., et al. "Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling." International Journal of Nonlinear Sciences and Numerical Simulation 10.3 (2009): 273-290. #' @references ter Braak C. J. F., and Vrugt J. A. (2008). Differential Evolution Markov Chain with snooker updater and fewer chains. Statistics and Computing http://dx.doi.org/10.1007/s11222-008-9104-9 -#' @details Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler +#' @details Instead of a bayesianSetup, the function can take the output of a previous run to restart the sampler #' from the last iteration. Due to the sampler's internal structure you can only use the output #' of DREAMzs. -#' If you provide a matrix with start values the number of rows detemines the number of chains that are run. -#' The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr +#' If you provide a matrix with start values, the number of rows determines the number of chains that will be run. +#' The number of columns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr #' There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly -#' the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +#' The algorithm implemented here does not have an automatic stopping criterion. Hence, it will #' always run the number of iterations specified by the user. Also, convergence is not #' monitored and left to the user. This can easily be done with coda::gelman.diag(chain). -#' Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here.\cr\cr +#' Furthermore, the delayed rejection step proposed in Vrugt et al. (2009) is not implemented here.\cr\cr #' During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. -#' First the disribution of crossover values is tuned to favor large jumps in the parameter space. +#' First, the distribution of crossover values is tuned to favor large jumps in the parameter space. #' The crossover probabilities determine how many parameters are updated simultaneously. -#' Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +#' Second, outlier chains are replaced as they can largely deteriorate the sampler's performance. #' However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain #' should be discarded when summarizing posterior moments. This can be done automatically during the -#' sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between -#' the burnin and adaptation phase to allow the user more flexibility in the sampler's settings. +#' sampling process (i.e. burn-in > adaptation) or subsequently by the user. We chose to distinguish between +#' the burn-in and adaptation phase to allow the user more flexibility in the sampler's settings. #' @example /inst/examples/DEfamilyHelp.R #' @seealso \code{\link{DREAM}} #' @export diff --git a/BayesianTools/R/mcmcFrancesco.R b/BayesianTools/R/mcmcFrancesco.R index be874ea..635a4cc 100644 --- a/BayesianTools/R/mcmcFrancesco.R +++ b/BayesianTools/R/mcmcFrancesco.R @@ -1,14 +1,14 @@ #' The Metropolis Algorithm #' @author Francesco Minunno #' @description The Metropolis Algorithm (Metropolis et al. 1953) -#' @param startValue vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior. -#' @param iterations iterations to run -#' @param nBI number of burnin +#' @param startValue vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case, startValues are sampled from the prior. +#' @param iterations number of iterations to run +#' @param nBI number of burn-in #' @param parmin minimum values for the parameter vector or NULL if FUN is of class BayesianSetup #' @param parmax maximum values for the parameter vector or NULL if FUN is of class BayesianSetup #' @param f scaling factor -#' @param FUN function to be sampled from or object of class bayesianSetup -#' @param consoleUpdates interger, determines the frequency with which sampler progress is printed to the console +#' @param FUN function to be sampled from or object of class BayesianSetup +#' @param consoleUpdates integer, determines the frequency with which sampler progress is printed to the console #' @references Metropolis, Nicholas, et al. "Equation of state calculations by fast computing machines." The journal of chemical physics 21.6 (1953): 1087-1092. #' @keywords internal # #' @export diff --git a/BayesianTools/R/mcmcMetropolis.R b/BayesianTools/R/mcmcMetropolis.R index 3a6a4bc..3d0677e 100644 --- a/BayesianTools/R/mcmcMetropolis.R +++ b/BayesianTools/R/mcmcMetropolis.R @@ -1,21 +1,21 @@ -#' Creates a Metropolis-type MCMC with options for covariance adaptatin, delayed rejection, Metropolis-within-Gibbs, and tempering +#' Creates a Metropolis-type MCMC with options for covariance adaptation, delayed rejection, Metropolis-within-Gibbs, and tempering #' @author Florian Hartig #' @param bayesianSetup either an object of class bayesianSetup created by \code{\link{createBayesianSetup}} (recommended), or a log target function -#' @param settings a list of settings - possible options follow below +#' @param settings a list of settings - possible options follow #' @param startValue startValue for the MCMC and optimization (if optimize = T). If not provided, the sampler will attempt to obtain the startValue from the bayesianSetup #' @param optimize logical, determines whether an optimization for start values and proposal function should be run before starting the sampling -#' @param proposalGenerator optional proposalgenerator object (see \code{\link{createProposalGenerator}}) +#' @param proposalGenerator optional, proposalgenerator object (see \code{\link{createProposalGenerator}}) #' @param proposalScaling additional scaling parameter for the proposals that controls the different scales of the proposals after delayed rejection (typical, after a rejection, one would want to try a smaller scale). Needs to be as long as DRlevels. Defaults to 0.5^(- 0:(mcmcSampler$settings$DRlevels -1) #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. #' @param thin thinning parameter. Determines the interval in which values are recorded. #' @param consoleUpdates integer, determines the frequency with which sampler progress is printed to the console -#' @param adapt logical, determines wheter an adaptive algorithm should be implemented. Default is TRUE. +#' @param adapt logical, determines whether an adaptive algorithm should be implemented. Default is TRUE. #' @param adaptationInterval integer, determines the interval of the adaption if adapt = TRUE. #' @param adaptationNotBefore integer, determines the start value for the adaption if adapt = TRUE. #' @param DRlevels integer, determines the number of levels for a delayed rejection sampler. Default is 1, which means no delayed rejection is used. #' @param temperingFunction function to implement simulated tempering in the algorithm. The function describes how the acceptance rate will be influenced in the course of the iterations. -#' @param gibbsProbabilities vector that defines the relative probabilities of the number of parameters to be changes simultaniously. -#' @param message logical determines whether the sampler's progress should be printed +#' @param gibbsProbabilities vector that defines the relative probabilities of the number of parameters to be changed simultaneously. +#' @param message logical, determines whether the sampler's progress should be printed #' @details The 'Metropolis' function is the main function for all Metropolis based samplers in this package. To call the derivatives from the basic Metropolis-Hastings MCMC, you can either use the corresponding function (e.g. \code{\link{AM}} for an adaptive Metropolis sampler) or use the parameters to adapt the basic Metropolis-Hastings. The advantage of the latter case is that you can easily combine different properties (e.g. adapive sampling and delayed rejection sampling) without changing the function. #' @import coda #' @example /inst/examples/MetropolisHelp.R diff --git a/BayesianTools/R/mcmcMultipleChains.R b/BayesianTools/R/mcmcMultipleChains.R index 84fa49c..0e9fab3 100644 --- a/BayesianTools/R/mcmcMultipleChains.R +++ b/BayesianTools/R/mcmcMultipleChains.R @@ -1,6 +1,6 @@ #' Run multiple chains -#' @param bayesianSetup Object of class "BayesianSetup" -#' @param settings list with settings for sampler +#' @param bayesianSetup object of class "BayesianSetup" +#' @param settings list containing settings for sampler #' @param sampler character, either "Metropolis" or "DE" #' @return list containing the single runs ($sampler) and the chains in a coda::mcmc.list ($mcmc.list) #' @keywords internal diff --git a/BayesianTools/R/mcmcRun.R b/BayesianTools/R/mcmcRun.R index b273a90..5258cef 100644 --- a/BayesianTools/R/mcmcRun.R +++ b/BayesianTools/R/mcmcRun.R @@ -24,7 +24,7 @@ #' #' Note that even if you specify parallel = T, this will only turn on internal parallelization of the samplers. The independent samplers controlled by nrChains are not evaluated in parallel, so if time is an issue it will be better to run the MCMCs individually and then combine them via \code{\link{createMcmcSamplerList}} into one joint object. #' -#' Note that DE and DREAM variants as well as SMC and T-walk require a population to start, which should be provided as a matrix. Default (NULL) sets the population size for DE to 3 x dimensions of parameters, for DREAM to 2 x dimensions of parameters and for DEzs and DREAMzs to three, sampled from the prior. Note also that the zs variants of DE and DREAM require two populations, the current population and the z matrix (a kind of memory) - if you want to set both, provide a list with startvalue$X and startvalue$Z. +#' Note that, DE and DREAM variants as well as SMC and T-walk require a population to start, which should be provided as a matrix. Default (NULL) sets the population size for DE to 3 x dimensions of parameters, for DREAM to 2 x dimensions of parameters and for DEzs and DREAMzs to three, sampled from the prior. Note also that the zs variants of DE and DREAM require two populations, the current population and the z matrix (a kind of memory) - if you want to set both, provide a list with startvalue$X and startvalue$Z. #' #' setting startValue for sampling with nrChains > 1 : if you want to provide different start values for the different chains, provide them as a list #' @@ -39,7 +39,7 @@ runMCMC <- function(bayesianSetup , sampler = "DEzs", settings = NULL){ ptm <- proc.time() - ####### RESTART ########## + ######## RESTART ########## if("bayesianOutput" %in% class(bayesianSetup)){ diff --git a/BayesianTools/R/mcmcTwalk.R b/BayesianTools/R/mcmcTwalk.R index ffceb46..2d311e3 100644 --- a/BayesianTools/R/mcmcTwalk.R +++ b/BayesianTools/R/mcmcTwalk.R @@ -1,19 +1,19 @@ #' T-walk MCMC #' @author Stefan Paul -#' @param bayesianSetup Object of class 'bayesianSetup' or 'bayesianOuput'. +#' @param bayesianSetup object of class 'BayesianSetup' or 'BayesianOuput'. #' @param settings list with parameter values. -#' @param iterations Number of model evaluations +#' @param iterations number of model evaluations #' @param at "traverse" move proposal parameter. Default to 6 #' @param aw "walk" move proposal parameter. Default to 1.5 -#' @param pn1 Probability determining the number of parameters that are changed -#' @param Ptrav Move probability of "traverse" moves, default to 0.4918 -#' @param Pwalk Move probability of "walk" moves, default to 0.4918 -#' @param Pblow Move probability of "traverse" moves, default to 0.0082 +#' @param pn1 probability determining the number of parameters that are changed +#' @param Ptrav move probability of "traverse" moves, default to 0.4918 +#' @param Pwalk move probability of "walk" moves, default to 0.4918 +#' @param Pblow move probability of "traverse" moves, default to 0.0082 #' @param burnin number of iterations treated as burn-in. These iterations are not recorded in the chain. #' @param thin thinning parameter. Determines the interval in which values are recorded. -#' @param startValue Matrix with start values -#' @param consoleUpdates Intervall in which the sampling progress is printed to the console -#' @param message logical determines whether the sampler's progress should be printed +#' @param startValue matrix with start values +#' @param consoleUpdates intervall in which the sampling progress is printed to the console +#' @param message logical, determines whether the sampler's progress should be printed #' @details ##' The probability of "hop" moves is 1 minus the sum of all other probabilities. #' @return Object of class bayesianOutput. diff --git a/BayesianTools/R/mcmcTwalk_helperFunctions.R b/BayesianTools/R/mcmcTwalk_helperFunctions.R index 8aaaf41..04dd9ef 100644 --- a/BayesianTools/R/mcmcTwalk_helperFunctions.R +++ b/BayesianTools/R/mcmcTwalk_helperFunctions.R @@ -1,10 +1,10 @@ ###### -# Twalk helper functions +# T-walk helper functions ###### #' Wrapper for step function -#' @param Npar Number of parameters -#' @param FUN Log posterior density +#' @param Npar number of parameters +#' @param FUN log posterior density #' @param x parameter vector of chain 1 #' @param Eval last evaluation of x #' @param x2 parameter vector of chain 2 diff --git a/BayesianTools/inst/examples/VSEMHelp.R b/BayesianTools/inst/examples/VSEMHelp.R index 83157fd..e1b5ff7 100644 --- a/BayesianTools/inst/examples/VSEMHelp.R +++ b/BayesianTools/inst/examples/VSEMHelp.R @@ -1,5 +1,5 @@ - + ## This example shows how to run and calibrate the VSEM model library(BayesianTools) diff --git a/BayesianTools/inst/examples/marginalLikelihoodHelp.R b/BayesianTools/inst/examples/marginalLikelihoodHelp.R index 171b693..4b6852d 100644 --- a/BayesianTools/inst/examples/marginalLikelihoodHelp.R +++ b/BayesianTools/inst/examples/marginalLikelihoodHelp.R @@ -76,7 +76,7 @@ exp(M1$ln.ML - M2$ln.ML) # it has therefore been suggested that ML should not be calculated on uninformative priors. But # what to do if there are no informative priors? -# one option is to calculate the fractional BF, which means that one splites the data in half, +# one option is to calculate the fractional BF, which means that one splits the data in half, # uses the first half to fit the model, and then use the posterior as a new (now informative) # prior for the ML - let's do this for the previous case @@ -145,7 +145,7 @@ exp(M1$ln.ML - M2$ln.ML) # Low dimensional case with narrow priors - all methods have low error # we use a truncated normal for the likelihood to make sure that the density -# integrates to 1 - makes it easier to calcuate the theoretical ML +# integrates to 1 - makes it easier to calculate the theoretical ML likelihood <- function(x) sum(msm::dtnorm(x, log = TRUE, lower = -1, upper = 1)) prior = createUniformPrior(lower = rep(-1,2), upper = rep(1,2)) bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior) @@ -180,4 +180,3 @@ marginalLikelihood(out, method = "HM", numSamples = 500)$ln.ML - theory marginalLikelihood(out, method = "Bridge", numSamples = 500)$ln.ML - theory - diff --git a/BayesianTools/inst/examples/plotSensitivityHelp.R b/BayesianTools/inst/examples/plotSensitivityHelp.R index 7afbdb5..523af53 100644 --- a/BayesianTools/inst/examples/plotSensitivityHelp.R +++ b/BayesianTools/inst/examples/plotSensitivityHelp.R @@ -1,5 +1,4 @@ - ll <- testDensityBanana bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 2), upper = rep(10, 2)) - plotSensitivity(bayesianSetup) + diff --git a/BayesianTools/inst/examples/plotTimeSeriesResultsHelp.R b/BayesianTools/inst/examples/plotTimeSeriesResultsHelp.R index 3991e0f..8109d31 100644 --- a/BayesianTools/inst/examples/plotTimeSeriesResultsHelp.R +++ b/BayesianTools/inst/examples/plotTimeSeriesResultsHelp.R @@ -25,7 +25,7 @@ likelihood <- function(par, sum = TRUE){ x[parSel] = par predicted <- VSEM(x[1:11], PAR) # replace here VSEM with your model predicted[,1] = 1000 * predicted[,1] # this is just rescaling - diff <- c(predicted[,1:4] - obs[,1:4]) # difference betweeno observed and predicted + diff <- c(predicted[,1:4] - obs[,1:4]) # difference between observed and predicted # univariate normal likelihood. Note that there is a parameter involved here that is fit llValues <- dnorm(diff, sd = x[12], log = TRUE) if (sum == FALSE) return(llValues) @@ -67,3 +67,4 @@ plotTimeSeriesResults(sampler = out, model = createPredictions, observed = obs[, error = createError, prior = TRUE, main = "Prior predictive") } + diff --git a/BayesianTools/inst/examples/proposalGeneratorHelp.R b/BayesianTools/inst/examples/proposalGeneratorHelp.R index 46dfebf..ffd94d6 100644 --- a/BayesianTools/inst/examples/proposalGeneratorHelp.R +++ b/BayesianTools/inst/examples/proposalGeneratorHelp.R @@ -59,3 +59,4 @@ testGenerator$returnProposal(testVector) x <- testGenerator$returnProposalMatrix(testMatrix) boxplot(x) table(x[,4]) + diff --git a/BayesianTools/inst/examples/testLinearModel.R b/BayesianTools/inst/examples/testLinearModel.R index 048c7d9..910dec5 100644 --- a/BayesianTools/inst/examples/testLinearModel.R +++ b/BayesianTools/inst/examples/testLinearModel.R @@ -1,3 +1,4 @@ + x = c(1,2) y = testLinearModel(x) plot(y) diff --git a/BayesianTools/inst/examples/tracePlotHelp.R b/BayesianTools/inst/examples/tracePlotHelp.R index d8f10ee..4de0405 100644 --- a/BayesianTools/inst/examples/tracePlotHelp.R +++ b/BayesianTools/inst/examples/tracePlotHelp.R @@ -1,3 +1,5 @@ + + # set up and run the MCMC ll <- function(x) sum(dnorm(x, log = TRUE)) setup <- createBayesianSetup(likelihood = ll, lower = c(-10, -10), upper = c(10,10)) diff --git a/BayesianTools/man/AdaptpCR.Rd b/BayesianTools/man/AdaptpCR.Rd index 3a7bfb4..1801699 100644 --- a/BayesianTools/man/AdaptpCR.Rd +++ b/BayesianTools/man/AdaptpCR.Rd @@ -7,13 +7,13 @@ AdaptpCR(CR, delta, lCR, settings, Npop) } \arguments{ -\item{CR}{Vector of crossover probabilities. Needs to be of length nCR.} +\item{CR}{vector of crossover probabilities. Needs to be of length nCR.} \item{delta}{vector with differences} \item{lCR}{values to weight delta} -\item{settings}{settings list} +\item{settings}{list of settings} \item{Npop}{number of chains.} } diff --git a/BayesianTools/man/DE.Rd b/BayesianTools/man/DE.Rd index a55b88f..9f35d67 100644 --- a/BayesianTools/man/DE.Rd +++ b/BayesianTools/man/DE.Rd @@ -17,7 +17,7 @@ DE( \item{settings}{list with parameter settings} -\item{startValue}{(optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population.} +\item{startValue}{(optional) either a matrix with start population, a number defining the number of chains to be run or a function that samples a starting population.} \item{iterations}{number of function evaluations.} @@ -31,7 +31,7 @@ DE( \item{blockUpdate}{list determining whether parameters should be updated in blocks. For possible settings see Details.} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, specifies whether to print the progress of the sampler.} } \description{ Differential-Evolution MCMC @@ -45,19 +45,7 @@ Possible choices are \item{"random"} { random blocking. Using k (see below)} \item{"user"} { user defined groups. Using groups (see below)} } -Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength -of the correlation used to group parameter and "groups" is used for user defined groups. -"groups" is a vector containing the group number for each parameter. E.g. for three parameters -with the first two in one group, "groups" would be c(1,1,2). -Further pSel and pGroup can be used to influence the choice of groups. In the sampling process -a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities -for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. -For updating one or two groups with the same probability pSel = c(1,1). By default all numbers -have the same probability. -The same principle is used in pGroup. Here the user can influence the probability of each group -to be updated. By default all groups have the same probability. -Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall -in which the groups are evaluated. +Further, seven parameters can be specified. "k" defines the number of groups, "h" the strength of the correlation used to group the parameters and "groups" is used for user defined groups. "groups" is a vector containing the group number for each parameter. E.g. for three parameters with the first two in one group, "groups" would be c(1,1,2). Moreover, pSel and pGroup can be used to influence the choice of groups. In the sampling process a number of groups are drawn at random and updated. pSel is a vector containing relative probabilities for updating the respective number of groups. E.g. To update one group at a time pSel = 1. For updating one or two groups with the same probability pSel = c(1,1). By default all numbers have the same probability. The same principle is used for pGroup. Here the user can influence the probability of each group to be updated. By default all groups have the same probability. Finally, "groupStart" defines the starting point of the groupUpdate and "groupIntervall" - the interval in which the groups are evaluated. } \examples{ library(BayesianTools) diff --git a/BayesianTools/man/DEzs.Rd b/BayesianTools/man/DEzs.Rd index df53820..37bc139 100644 --- a/BayesianTools/man/DEzs.Rd +++ b/BayesianTools/man/DEzs.Rd @@ -18,11 +18,11 @@ DEzs( \item{settings}{list with parameter settings} -\item{startValue}{(optional) eiter a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population.} +\item{startValue}{(optional) either a matrix with start population, a number to define the number of chains that are run or a function that samples a starting population.} \item{Z}{starting Z population} -\item{iterations}{iterations to run} +\item{iterations}{number of iterations to run} \item{pSnooker}{probability of Snooker update} @@ -44,13 +44,13 @@ DEzs( \item{blockUpdate}{list determining whether parameters should be updated in blocks. For possible settings see Details.} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, specifies whether to print the progress of the sampler.} } \description{ Differential-Evolution MCMC zs } \details{ -For parallel computing, the likelihood density in the bayesianSetup needs to be parallelized, i.e. needs to be able to operate on a matrix of proposals +For parallel computing, the likelihood density in the bayesianSetup needs to be parallelized, i.e., it needs to be able to operate on a matrix of proposals For blockUpdate the first element in the list determines the type of blocking. Possible choices are @@ -60,18 +60,18 @@ Possible choices are \item{"random"} { random blocking. Using k (see below)} \item{"user"} { user defined groups. Using groups (see below)} } -Further seven parameters can be specified. "k" determnined the number of groups, "h" the strength +Further, seven parameters can be specified. "k" defines the number of groups, "h" the strength of the correlation used to group parameter and "groups" is used for user defined groups. "groups" is a vector containing the group number for each parameter. E.g. for three parameters with the first two in one group, "groups" would be c(1,1,2). -Further pSel and pGroup can be used to influence the choice of groups. In the sampling process -a number of groups is randomly drawn and updated. pSel is a vector containing relative probabilities -for an update of the respective number of groups. E.g. for always updating only one group pSel = 1. +Moreover, pSel and pGroup can be used to influence the choice of groups. In the sampling process +a number of groups is drawn at random and updated. pSel is a vector containing relative probabilities +for updating the respective number of groups. E.g. To update one group at a time pSel = 1. For updating one or two groups with the same probability pSel = c(1,1). By default all numbers have the same probability. -The same principle is used in pGroup. Here the user can influence the probability of each group +The same principle is used in pGroup. Here, the user can influence the probability of each group to be updated. By default all groups have the same probability. -Finally "groupStart" defines the starting point of the groupUpdate and "groupIntervall" the intervall +Finally, "groupStart" defines the starting point of the groupUpdate and "groupIntervall" - the interval in which the groups are evaluated. } \examples{ diff --git a/BayesianTools/man/DREAM.Rd b/BayesianTools/man/DREAM.Rd index e6da437..c509399 100644 --- a/BayesianTools/man/DREAM.Rd +++ b/BayesianTools/man/DREAM.Rd @@ -13,15 +13,15 @@ DREAM( ) } \arguments{ -\item{bayesianSetup}{Object of class 'bayesianSetup' or 'bayesianOuput'.} +\item{bayesianSetup}{object of class 'bayesianSetup' or 'bayesianOuput'.} \item{settings}{list with parameter values} -\item{iterations}{Number of model evaluations} +\item{iterations}{number of model evaluations} \item{nCR}{parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly.} -\item{updateInterval}{determining the intervall for the pCR update} +\item{updateInterval}{determines the interval for the pCR update} \item{gamma}{Kurtosis parameter Bayesian Inference Scheme} @@ -29,21 +29,21 @@ DREAM( \item{e}{Ergodicity term} -\item{pCRupdate}{If T, crossover probabilities will be updated} +\item{pCRupdate}{logical, if T, crossover probabilities will be updated} \item{burnin}{number of iterations treated as burn-in. These iterations are not recorded in the chain.} -\item{thin}{thin thinning parameter. Determines the interval in which values are recorded.} +\item{thin}{thinning parameter. Determines the interval in which values are recorded.} -\item{adaptation}{Number or percentage of samples that are used for the adaptation in DREAM (see Details).} +\item{adaptation}{number or percentage of samples that are used for the adaptation in DREAM (see Details).} -\item{DEpairs}{Number of pairs used to generate proposal} +\item{DEpairs}{number of pairs used to generate proposal} -\item{startValue}{eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior.} +\item{startValue}{either a matrix containing the start values (see details), an integer to define the number of chains to be run, a function to sample the start values or NUll - in which case the values are sampled from the prior.} -\item{consoleUpdates}{Intervall in which the sampling progress is printed to the console} +\item{consoleUpdates}{interval at which the sampling progress is printed to the console} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, determines whether the sampler's progress should be printed} } \value{ mcmc.object containing the following elements: chains, X, pCR @@ -52,25 +52,24 @@ mcmc.object containing the following elements: chains, X, pCR DREAM } \details{ -Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler -from the last iteration. Due to the sampler's internal structure you can only use the output -of DREAM. -If you provide a matrix with start values the number of rows determines the number of chains that are run. +Instead of a bayesianSetup, the function can take the output of a previous run to restart the sampler +from the last iteration. Due to the sampler's internal structure you can only use the output of DREAM. +If you provide a matrix with start values, the number of rows determines the number of chains that will be run. The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly -the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +The algorithm implemented here does not have an automatic stopping criterion. Hence, it will always run the number of iterations specified by the user. Also, convergence is not monitored and left to the user. This can easily be done with coda::gelman.diag(chain). -Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here.\cr\cr +Furthermore, the delayed rejection step proposed in Vrugt et al. (2009) is not implemented here.\cr\cr During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. -First the disribution of crossover values is tuned to favor large jumps in the parameter space. +First, the disribution of crossover values is tuned to favor large jumps in the parameter space. The crossover probabilities determine how many parameters are updated simultaneously. -Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +Second, outlier chains are replaced as they can largely deteriorate the sampler's performance. However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain should be discarded when summarizing posterior moments. This can be done automatically during the -sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between -the burnin and adaptation phase to allow the user more flexibility in the sampler's settings. +sampling process (i.e. burn-in > adaptation) or subsequently by the user. We chose to distinguish between +the burn-in and adaptation phase to allow the user more flexibility in the sampler's settings. } \examples{ library(BayesianTools) diff --git a/BayesianTools/man/DREAMzs.Rd b/BayesianTools/man/DREAMzs.Rd index 80f1805..87fe481 100644 --- a/BayesianTools/man/DREAMzs.Rd +++ b/BayesianTools/man/DREAMzs.Rd @@ -13,31 +13,31 @@ DREAMzs( ) } \arguments{ -\item{bayesianSetup}{Object of class 'bayesianSetup' or 'bayesianOuput'.} +\item{bayesianSetup}{object of class 'bayesianSetup' or 'bayesianOuput'.} \item{settings}{list with parameter values} -\item{iterations}{Number of model evaluations} +\item{iterations}{number of model evaluations} -\item{nCR}{parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly.} +\item{nCR}{parameter determining the number of cross-over proposals. If nCR = 1 all parameters are updated jointly.} -\item{updateInterval}{determining the intervall for the pCR (crossover probabilities) update} +\item{updateInterval}{determines the interval for the pCR (crossover probabilities) update} -\item{gamma}{Kurtosis parameter Bayesian Inference Scheme.} +\item{gamma}{kurtosis parameter for Bayesian Inference Scheme.} \item{eps}{Ergodicity term} \item{e}{Ergodicity term} -\item{pCRupdate}{Update of crossover probabilities} +\item{pCRupdate}{update of crossover probabilities} \item{burnin}{number of iterations treated as burn-in. These iterations are not recorded in the chain.} -\item{thin}{thin thinning parameter. Determines the interval in which values are recorded.} +\item{thin}{thinning parameter. Determines the interval in which values are recorded.} -\item{adaptation}{Number or percentage of samples that are used for the adaptation in DREAM (see Details)} +\item{adaptation}{number or percentage of samples that are used for the adaptation in DREAM (see Details)} -\item{DEpairs}{Number of pairs used to generate proposal} +\item{DEpairs}{number of pairs used to generate proposal} \item{ZupdateFrequency}{frequency to update Z matrix} @@ -45,11 +45,11 @@ DREAMzs( \item{Z}{starting matrix for Z} -\item{startValue}{eiter a matrix containing the start values (see details), an integer to define the number of chains that are run, a function to sample the start values or NUll, in which case the values are sampled from the prior.} +\item{startValue}{either a matrix containing the start values (see details), an integer to define the number of chains to be run, a function to sample the start values or NUll - in which case the values are sampled from the prior.} -\item{consoleUpdates}{Intervall in which the sampling progress is printed to the console} +\item{consoleUpdates}{interval in which the sampling progress is printed to the console} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, determines whether the sampler's progress should be printed} } \value{ mcmc.object containing the following elements: chains, X, pCR, Z @@ -58,24 +58,24 @@ mcmc.object containing the following elements: chains, X, pCR, Z DREAMzs } \details{ -Insted of a bayesianSetup, the function can take the output of a previous run to restart the sampler +Instead of a bayesianSetup, the function can take the output of a previous run to restart the sampler from the last iteration. Due to the sampler's internal structure you can only use the output of DREAMzs. -If you provide a matrix with start values the number of rows detemines the number of chains that are run. -The number of coloumns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr +If you provide a matrix with start values, the number of rows determines the number of chains that will be run. +The number of columns must be equivalent to the number of parameters in your bayesianSetup. \cr\cr There are several small differences in the algorithm presented here compared to the original paper by Vrugt et al. (2009). Mainly -the algorithm implemented here does not have an automatic stopping criterion. Hence, it will +The algorithm implemented here does not have an automatic stopping criterion. Hence, it will always run the number of iterations specified by the user. Also, convergence is not monitored and left to the user. This can easily be done with coda::gelman.diag(chain). -Further the proposed delayed rejectio step in Vrugt et al. (2009) is not implemented here.\cr\cr +Furthermore, the delayed rejection step proposed in Vrugt et al. (2009) is not implemented here.\cr\cr During the adaptation phase DREAM is running two mechanisms to enhance the sampler's efficiency. -First the disribution of crossover values is tuned to favor large jumps in the parameter space. +First, the distribution of crossover values is tuned to favor large jumps in the parameter space. The crossover probabilities determine how many parameters are updated simultaneously. -Second outlier chains are replanced as they can largely deteriorate the sampler's performance. +Second, outlier chains are replaced as they can largely deteriorate the sampler's performance. However, these steps destroy the detailed balance of the chain. Consequently these parts of the chain should be discarded when summarizing posterior moments. This can be done automatically during the -sampling process (i.e. burnin > adaptation) or subsequently by the user. We chose to distinguish between -the burnin and adaptation phase to allow the user more flexibility in the sampler's settings. +sampling process (i.e. burn-in > adaptation) or subsequently by the user. We chose to distinguish between +the burn-in and adaptation phase to allow the user more flexibility in the sampler's settings. } \examples{ library(BayesianTools) diff --git a/BayesianTools/man/M.Rd b/BayesianTools/man/M.Rd index d365647..3cfd803 100644 --- a/BayesianTools/man/M.Rd +++ b/BayesianTools/man/M.Rd @@ -16,11 +16,11 @@ M( ) } \arguments{ -\item{startValue}{vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case startValues are sampled from the prior.} +\item{startValue}{vector with the start values for the algorithm. Can be NULL if FUN is of class BayesianSetup. In this case, startValues are sampled from the prior.} -\item{iterations}{iterations to run} +\item{iterations}{number of iterations to run} -\item{nBI}{number of burnin} +\item{nBI}{number of burn-in} \item{parmin}{minimum values for the parameter vector or NULL if FUN is of class BayesianSetup} @@ -28,9 +28,9 @@ M( \item{f}{scaling factor} -\item{FUN}{function to be sampled from or object of class bayesianSetup} +\item{FUN}{function to be sampled from or object of class BayesianSetup} -\item{consoleUpdates}{interger, determines the frequency with which sampler progress is printed to the console} +\item{consoleUpdates}{integer, determines the frequency with which sampler progress is printed to the console} } \description{ The Metropolis Algorithm (Metropolis et al. 1953) diff --git a/BayesianTools/man/MAP.Rd b/BayesianTools/man/MAP.Rd index ce4490a..f79f6d3 100644 --- a/BayesianTools/man/MAP.Rd +++ b/BayesianTools/man/MAP.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/MAP.R \name{MAP} \alias{MAP} -\title{calculates the Maxiumum APosteriori value (MAP)} +\title{Calculates the Maximum APosteriori value (MAP)} \usage{ MAP(bayesianOutput, ...) } @@ -12,10 +12,10 @@ MAP(bayesianOutput, ...) \item{...}{optional values to be passed on the the getSample function} } \description{ -calculates the Maxiumum APosteriori value (MAP) +Calculates the Maximum APosteriori value (MAP) } \details{ -Currently, this function simply returns the parameter combination with the highest posterior in the chain. A more refined option would be to take the MCMC sample and do additional calculations, e.g. use an optimizer, a kerne delnsity estimator, or some other tool to search / interpolate around the best value in the chain +Currently, this function simply returns the parameter combination with the highest posterior in the chain. A more refined option would be to take the MCMC sample and do additional calculations, e.g. use an optimizer, a kernel density estimator, or some other tool to search / interpolate around the best value in the chain. } \seealso{ \code{\link{WAIC}}, \code{\link{DIC}}, \code{\link{marginalLikelihood}} diff --git a/BayesianTools/man/Metropolis.Rd b/BayesianTools/man/Metropolis.Rd index d17408d..8a032a1 100644 --- a/BayesianTools/man/Metropolis.Rd +++ b/BayesianTools/man/Metropolis.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/mcmcMetropolis.R \name{Metropolis} \alias{Metropolis} -\title{Creates a Metropolis-type MCMC with options for covariance adaptatin, delayed rejection, Metropolis-within-Gibbs, and tempering} +\title{Creates a Metropolis-type MCMC with options for covariance adaptation, delayed rejection, Metropolis-within-Gibbs, and tempering} \usage{ Metropolis( bayesianSetup, @@ -16,13 +16,13 @@ Metropolis( \arguments{ \item{bayesianSetup}{either an object of class bayesianSetup created by \code{\link{createBayesianSetup}} (recommended), or a log target function} -\item{settings}{a list of settings - possible options follow below} +\item{settings}{a list of settings - possible options follow} \item{startValue}{startValue for the MCMC and optimization (if optimize = T). If not provided, the sampler will attempt to obtain the startValue from the bayesianSetup} \item{optimize}{logical, determines whether an optimization for start values and proposal function should be run before starting the sampling} -\item{proposalGenerator}{optional proposalgenerator object (see \code{\link{createProposalGenerator}})} +\item{proposalGenerator}{optional, proposalgenerator object (see \code{\link{createProposalGenerator}})} \item{proposalScaling}{additional scaling parameter for the proposals that controls the different scales of the proposals after delayed rejection (typical, after a rejection, one would want to try a smaller scale). Needs to be as long as DRlevels. Defaults to 0.5^(- 0:(mcmcSampler$settings$DRlevels -1)} @@ -32,7 +32,7 @@ Metropolis( \item{consoleUpdates}{integer, determines the frequency with which sampler progress is printed to the console} -\item{adapt}{logical, determines wheter an adaptive algorithm should be implemented. Default is TRUE.} +\item{adapt}{logical, determines whether an adaptive algorithm should be implemented. Default is TRUE.} \item{adaptationInterval}{integer, determines the interval of the adaption if adapt = TRUE.} @@ -42,12 +42,12 @@ Metropolis( \item{temperingFunction}{function to implement simulated tempering in the algorithm. The function describes how the acceptance rate will be influenced in the course of the iterations.} -\item{gibbsProbabilities}{vector that defines the relative probabilities of the number of parameters to be changes simultaniously.} +\item{gibbsProbabilities}{vector that defines the relative probabilities of the number of parameters to be changed simultaneously.} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, determines whether the sampler's progress should be printed} } \description{ -Creates a Metropolis-type MCMC with options for covariance adaptatin, delayed rejection, Metropolis-within-Gibbs, and tempering +Creates a Metropolis-type MCMC with options for covariance adaptation, delayed rejection, Metropolis-within-Gibbs, and tempering } \details{ The 'Metropolis' function is the main function for all Metropolis based samplers in this package. To call the derivatives from the basic Metropolis-Hastings MCMC, you can either use the corresponding function (e.g. \code{\link{AM}} for an adaptive Metropolis sampler) or use the parameters to adapt the basic Metropolis-Hastings. The advantage of the latter case is that you can easily combine different properties (e.g. adapive sampling and delayed rejection sampling) without changing the function. diff --git a/BayesianTools/man/Twalk.Rd b/BayesianTools/man/Twalk.Rd index c32caf6..b94d361 100644 --- a/BayesianTools/man/Twalk.Rd +++ b/BayesianTools/man/Twalk.Rd @@ -12,33 +12,33 @@ Twalk( ) } \arguments{ -\item{bayesianSetup}{Object of class 'bayesianSetup' or 'bayesianOuput'.} +\item{bayesianSetup}{object of class 'BayesianSetup' or 'BayesianOuput'.} \item{settings}{list with parameter values.} -\item{iterations}{Number of model evaluations} +\item{iterations}{number of model evaluations} \item{at}{"traverse" move proposal parameter. Default to 6} \item{aw}{"walk" move proposal parameter. Default to 1.5} -\item{pn1}{Probability determining the number of parameters that are changed} +\item{pn1}{probability determining the number of parameters that are changed} -\item{Ptrav}{Move probability of "traverse" moves, default to 0.4918} +\item{Ptrav}{move probability of "traverse" moves, default to 0.4918} -\item{Pwalk}{Move probability of "walk" moves, default to 0.4918} +\item{Pwalk}{move probability of "walk" moves, default to 0.4918} -\item{Pblow}{Move probability of "traverse" moves, default to 0.0082} +\item{Pblow}{move probability of "traverse" moves, default to 0.0082} \item{burnin}{number of iterations treated as burn-in. These iterations are not recorded in the chain.} \item{thin}{thinning parameter. Determines the interval in which values are recorded.} -\item{startValue}{Matrix with start values} +\item{startValue}{matrix with start values} -\item{consoleUpdates}{Intervall in which the sampling progress is printed to the console} +\item{consoleUpdates}{intervall in which the sampling progress is printed to the console} -\item{message}{logical determines whether the sampler's progress should be printed} +\item{message}{logical, determines whether the sampler's progress should be printed} } \value{ Object of class bayesianOutput. diff --git a/BayesianTools/man/TwalkMove.Rd b/BayesianTools/man/TwalkMove.Rd index 8639ad0..06ffc78 100644 --- a/BayesianTools/man/TwalkMove.Rd +++ b/BayesianTools/man/TwalkMove.Rd @@ -21,9 +21,9 @@ TwalkMove( ) } \arguments{ -\item{Npar}{Number of parameters} +\item{Npar}{number of parameters} -\item{FUN}{Log posterior density} +\item{FUN}{log posterior density} \item{x}{parameter vector of chain 1} diff --git a/BayesianTools/man/VSEM.Rd b/BayesianTools/man/VSEM.Rd index 84571d0..cb73dfe 100644 --- a/BayesianTools/man/VSEM.Rd +++ b/BayesianTools/man/VSEM.Rd @@ -86,7 +86,7 @@ NEE Net Ecosystem Exchange kg C /m2 /day } \examples{ - + ## This example shows how to run and calibrate the VSEM model library(BayesianTools) diff --git a/BayesianTools/man/VSEMcreateLikelihood.Rd b/BayesianTools/man/VSEMcreateLikelihood.Rd index 0f29233..88e6a9b 100644 --- a/BayesianTools/man/VSEMcreateLikelihood.Rd +++ b/BayesianTools/man/VSEMcreateLikelihood.Rd @@ -7,7 +7,7 @@ VSEMcreateLikelihood(likelihoodOnly = F, plot = F, selection = c(1:6, 12)) } \arguments{ -\item{likelihoodOnly}{switch to devide whether to create only a likelihood, or a full bayesianSetup with uniform priors.} +\item{likelihoodOnly}{switch to decide whether to create only a likelihood, or a full bayesianSetup with uniform priors.} \item{plot}{switch to decide whether data should be plotted} diff --git a/BayesianTools/man/combineChains.Rd b/BayesianTools/man/combineChains.Rd index 85fcedb..814473f 100644 --- a/BayesianTools/man/combineChains.Rd +++ b/BayesianTools/man/combineChains.Rd @@ -9,7 +9,7 @@ combineChains(x, merge = T) \arguments{ \item{x}{a list of MCMC chains} -\item{merge}{logical determines whether chains should be merged} +\item{merge}{logical, should chains be merged? (T or F)} } \value{ combined chains diff --git a/BayesianTools/man/convertCoda.Rd b/BayesianTools/man/convertCoda.Rd index 6009f3b..e4213b0 100644 --- a/BayesianTools/man/convertCoda.Rd +++ b/BayesianTools/man/convertCoda.Rd @@ -7,22 +7,19 @@ convertCoda(sampler, names = NULL, info = NULL, likelihood = NULL) } \arguments{ -\item{sampler}{An object of class mcmc or mcmc.list} +\item{sampler}{an object of class mcmc or mcmc.list} -\item{names}{vector giving the parameter names (optional)} +\item{names}{a vector with parameter names (optional)} -\item{info}{matrix (or list with matrices for mcmc.list objects) with three coloumns containing log posterior, log likelihood and log prior of the sampler for each time step (optional; but see Details)} +\item{info}{a matrix (or list with matrices for mcmc.list objects) with three columns containing log posterior, log likelihood and log prior of the sampler for each time step (optional; but see Details)} -\item{likelihood}{likelihood function used in the sampling (see Details)} +\item{likelihood}{likelihood function used for sampling (see Details)} } \description{ -Function is used to make the plot and diagnostic functions -available for coda::mcmc objects +Function to support plotting and diagnostic functions for coda::mcmc objects. } \details{ -The parameter 'likelihood' is optional for most functions but can be needed e.g for -using the \code{\link{DIC}} function. +The parameter 'likelihood' is optional for most functions but can be needed e.g for \code{\link{DIC}} function. -Also the parameter info is optional for most uses. However for some functions (e.g. \code{\link{MAP}}) -the matrix or single coloumns (e.g. log posterior) are necessary for the diagnostics. +Also, the parameter information is typically optional for most uses. However, for certain functions (e.g. \code{\link{MAP}}), the matrix or single columns (e.g. log posterior) are necessary for diagnostics. } diff --git a/BayesianTools/man/createBayesianSetup.Rd b/BayesianTools/man/createBayesianSetup.Rd index 3225342..df98527 100644 --- a/BayesianTools/man/createBayesianSetup.Rd +++ b/BayesianTools/man/createBayesianSetup.Rd @@ -37,9 +37,9 @@ createBayesianSetup( \item{names}{optional vector with parameter names} -\item{parallelOptions}{list containing three lists. First "packages" determines the R packages necessary to run the likelihood function. Second "variables" the objects in the global environment needed to run the likelihood function and third "dlls" the DLLs needed to run the likelihood function (see Details and Examples).} +\item{parallelOptions}{list containing three lists.\itemize{ \item First, "packages" determines the R packages necessary to run the likelihood function.\item Second, "variables" - the objects in the global environment needed to run the likelihood function and \item Third, "dlls" is needed to run the likelihood function (see Details and Examples). }} -\item{catchDuplicates}{Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns.} +\item{catchDuplicates}{logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns.} \item{plotLower}{vector with lower limits for plotting} @@ -53,14 +53,14 @@ Creates a standardized collection of prior, likelihood and posterior functions, \details{ If prior is of class prior (e.g. create with \code{\link{createPrior}}), priorSampler, lower, upper and best will be ignored.\cr If prior is a function (log prior density), priorSampler (custom sampler), or lower/upper (uniform sampler) is required.\cr If prior is NULL, and lower and upper are passed, a uniform prior (see \code{\link{createUniformPrior}}) will be created with boundaries lower and upper. -For parallelization, Bayesiantools requies that the likelihood can evaluate several parameter vectors (supplied as a matrix) in parallel. +For parallelization, Bayesiantools requires that the likelihood can evaluate multiple parameter vectors (supplied as a matrix) in parallel. \itemize{ -\item parallel = T means that an automatic parallelization of the likelihood via a standard R socket cluster is attempted, using the function \code{\link{generateParallelExecuter}}. By default, of the N cores detected on the computer, N-1 cores are requested. Alternatively, you can provide a integer number to parallel, specifying the cores reserved for the cluster. When the cluster is cluster is created, a copy of your workspace, including DLLs and objects are exported to the cluster workers. Because this can be very inefficient, you can explicitly specify the packages, objects and DLLs that are to be exported via parallelOptions. Using parallel = T requires that the function to be parallelized is well encapsulate, i.e. can run on a shared memory / shared hard disk machine in parallel without interfering with each other. +\item parallel = T attempts to parallelize likelihood via a standard R socket cluster using the \code{\link{generateParallelExecuter}} function. By default, of the N cores detected on the computer, N-1 cores are requested. Alternatively, you can provide a integer number to parallel, specifying the cores reserved for the cluster. When the cluster is created, a copy of your workspace, including DLLs and objects are exported to the cluster workers. As this approach can be highly inefficient, it is recommended to explicitly specify the packages, objects and DLLs to export using parallelOptions. Using parallel = T requires that the function to be parallelized is well encapsulated, i.e. can run in parallel on a shared memory / shared hard disk machine in parallel without interfering with each other. } -If automatic parallelization cannot be done (e.g. because dlls are not thread-safe or write to shared disk), and only in this case, you should specify parallel = "external". In this case, it is assumed that the likelihood is programmed such that it accepts a matrix with parameters as columns and the different model runs as rows. It is then up to the user if and how to parallelize this function. This option gives most flexibility to the user, in particular for complicated parallel architecture or shared memory problems. +If automatic parallelization is not possible (e.g., because dlls are not thread-safe or write to shared disk), and only in this case, you should specify parallel = "external". In this case, it is assumed that the likelihood is programmed to accept a matrix with parameters as columns and the different model runs as rows. The user can then choose whether and how to parallelize this function. This option provides optimal flexibility for the user, especially regarding complicated parallel architectures or shared memory issues. -For more details on parallelization, make sure to read both vignettes, in particular the section on the likelihood in the main vignette, and the section on parallelization in the vignette on interfacing models. +For more details on parallelization, make sure to read both vignettes, especially the section on likelihood in the main vignette and the section on parallelization in the vignette on interfacing models. } \examples{ ll <- function(x) sum(dnorm(x, log = TRUE)) diff --git a/BayesianTools/man/createLikelihood.Rd b/BayesianTools/man/createLikelihood.Rd index 9b9d33b..b9ad678 100644 --- a/BayesianTools/man/createLikelihood.Rd +++ b/BayesianTools/man/createLikelihood.Rd @@ -14,17 +14,17 @@ createLikelihood( ) } \arguments{ -\item{likelihood}{Log likelihood density} +\item{likelihood}{log likelihood density} -\item{names}{Parameter names (optional)} +\item{names}{parameter names (optional)} -\item{parallel}{parallelization , either i) no parallelization --> F, ii) native R parallelization --> T / "auto" will select n-1 of your available cores, or provide a number for how many cores to use, or iii) external parallelization --> "external". External means that the likelihood is already able to execute parallel runs in form of a matrix with} +\item{parallel}{parallelization , either i) no parallelization --> F, ii) native R parallelization --> T / "auto" will select n-1 of your available cores, or provide a number for how many cores to use, or iii) external parallelization --> "external". External means that the likelihood is already able to execute parallel runs in the form of a matrix.} -\item{catchDuplicates}{Logical, determines whether unique parameter combinations should only be evaluated once. Only used when the likelihood accepts a matrix with parameter as columns.} +\item{catchDuplicates}{logical, determines whether unique parameter combinations should only be evaluated once. This is only applicable when the likelihood accepts a matrix with parameters as columns.} \item{sampler}{sampler} -\item{parallelOptions}{list containing two lists. First "packages" determines the R packages necessary to run the likelihood function. Second "objects" the objects in the global envirnment needed to run the likelihood function (for details see \code{\link{createBayesianSetup}}).} +\item{parallelOptions}{a list containing two lists. First, "packages" specifies the R packages necessary to run the likelihood function. Second, "objects" contains the objects in the global environment needed to run the likelihood function (for details see \code{\link{createBayesianSetup}}).} } \description{ Creates a standardized likelihood class#' diff --git a/BayesianTools/man/createMcmcSamplerList.Rd b/BayesianTools/man/createMcmcSamplerList.Rd index a81dc76..5f3ebde 100644 --- a/BayesianTools/man/createMcmcSamplerList.Rd +++ b/BayesianTools/man/createMcmcSamplerList.Rd @@ -7,10 +7,10 @@ createMcmcSamplerList(mcmcList) } \arguments{ -\item{mcmcList}{a list with each object being an mcmcSampler} +\item{mcmcList}{list of objects, each of which is an mcmcSampler} } \value{ -Object of class "mcmcSamplerList" +object of class "mcmcSamplerList" } \description{ Convenience function to create an object of class mcmcSamplerList from a list of mcmc samplers diff --git a/BayesianTools/man/createPosterior.Rd b/BayesianTools/man/createPosterior.Rd index d197c47..7686af3 100644 --- a/BayesianTools/man/createPosterior.Rd +++ b/BayesianTools/man/createPosterior.Rd @@ -9,13 +9,13 @@ createPosterior(prior, likelihood) \arguments{ \item{prior}{prior class} -\item{likelihood}{Log likelihood density} +\item{likelihood}{log likelihood density} } \description{ Creates a standardized posterior class } \details{ -Function is internally used in \code{\link{createBayesianSetup}} to create a standarized posterior class. +Function is internally used in \code{\link{createBayesianSetup}} to create a standardized posterior class. } \author{ Florian Hartig diff --git a/BayesianTools/man/createPrior.Rd b/BayesianTools/man/createPrior.Rd index 666dd9b..4718810 100644 --- a/BayesianTools/man/createPrior.Rd +++ b/BayesianTools/man/createPrior.Rd @@ -13,7 +13,7 @@ createPrior( ) } \arguments{ -\item{density}{Prior density} +\item{density}{prior density} \item{sampler}{Sampling function for density (optional)} @@ -27,7 +27,7 @@ createPrior( Creates a standardized prior class } \details{ -This is the general prior generator. It is highly recommended to not only implement the density, but also the sampler function. If this is not done, the user will have to provide explicit starting values for many of the MCMC samplers. Note the existing, more specialized prior function. If your prior can be created by those, they are preferred. Note also that priors can be created from an existing MCMC output from BT, or another MCMC sample, via \code{\link{createPriorDensity}}. +This is the general prior generator. It is highly recommended to implement both the density and sampler function. If not, the user will have to provide explicit starting values for many of the MCMC samplers. Note the existing, more specialized prior functions. It is recommended to use those specialized prior functions, if possible. Also note that priors can be created from an existing MCMC output from BT, or another MCMC sample, via \code{\link{createPriorDensity}}. } \note{ min and max truncate, but not re-normalize the prior density (so, if a pdf that integrated to one is truncated, the integral will in general be smaller than one). For MCMC sampling, this doesn't make a difference, but if absolute values of the prior density are a concern, one should provide a truncated density function for the prior. diff --git a/BayesianTools/man/createProposalGenerator.Rd b/BayesianTools/man/createProposalGenerator.Rd index f0e28d4..dbab345 100644 --- a/BayesianTools/man/createProposalGenerator.Rd +++ b/BayesianTools/man/createProposalGenerator.Rd @@ -100,6 +100,7 @@ testGenerator$returnProposal(testVector) x <- testGenerator$returnProposalMatrix(testMatrix) boxplot(x) table(x[,4]) + } \seealso{ \code{\link{updateProposalGenerator}} diff --git a/BayesianTools/man/createSmcSamplerList.Rd b/BayesianTools/man/createSmcSamplerList.Rd index 11208b2..735ce36 100644 --- a/BayesianTools/man/createSmcSamplerList.Rd +++ b/BayesianTools/man/createSmcSamplerList.Rd @@ -10,7 +10,7 @@ createSmcSamplerList(...) \item{...}{a list of MCMC samplers} } \value{ -a list of class smcSamplerList with each object being an smcSampler +a list of class smcSamplerList with objects of smcSampler } \description{ Convenience function to create an object of class SMCSamplerList from a list of mcmc samplers diff --git a/BayesianTools/man/generateCRvalues.Rd b/BayesianTools/man/generateCRvalues.Rd index 41e6302..6f0a81c 100644 --- a/BayesianTools/man/generateCRvalues.Rd +++ b/BayesianTools/man/generateCRvalues.Rd @@ -7,9 +7,9 @@ generateCRvalues(pCR, settings, Npop) } \arguments{ -\item{pCR}{Vector of crossover probabilities. Needs to be of length nCR.} +\item{pCR}{vector of crossover probabilities. Needs to be of length nCR.} -\item{settings}{settings list} +\item{settings}{list of settings} \item{Npop}{number of chains} } diff --git a/BayesianTools/man/getBlockSettings.Rd b/BayesianTools/man/getBlockSettings.Rd index 6c29a25..a0bd314 100644 --- a/BayesianTools/man/getBlockSettings.Rd +++ b/BayesianTools/man/getBlockSettings.Rd @@ -13,6 +13,6 @@ getBlockSettings(blockUpdate) list with block settings } \description{ -Transforms the original settings in settings used in the model runs +Transforms the original settings to settings used in the model runs } \keyword{internal} diff --git a/BayesianTools/man/getSample.Rd b/BayesianTools/man/getSample.Rd index 900847f..d72d774 100644 --- a/BayesianTools/man/getSample.Rd +++ b/BayesianTools/man/getSample.Rd @@ -148,15 +148,15 @@ getSample( \item{parametersOnly}{for a BT output, if F, likelihood, posterior and prior values are also provided in the output} -\item{coda}{works only for mcmc classes - provides output as a coda object. Note: if mcmcSamplerList contains mcmc samplers such as DE that have several chains, the internal chains will be collapsed. This may not be the desired behavior for all applications.} +\item{coda}{works only for mcmc classes - returns output as a coda object. Note: if mcmcSamplerList contains mcmc samplers such as DE that have several chains, the internal chains will be collapsed. This may not be desired for all applications.} -\item{start}{for mcmc samplers start value in the chain. For SMC samplers, start particle} +\item{start}{for mcmc samplers, start value in the chain. For SMC samplers, start particle} \item{end}{for mcmc samplers end value in the chain. For SMC samplers, end particle} -\item{thin}{thinning parameter. Either an integer determining the thinning intervall (default is 1) or "auto" for automatic thinning.} +\item{thin}{thinning parameter. Either an integer determining the thinning interval (default is 1) or "auto" for automatic thinning.} -\item{numSamples}{sample size (only used if thin = 1). If you want to use numSamples set thin to 1.} +\item{numSamples}{sample size (only used if thin = 1). If you want to use numSamples, set thin to 1.} \item{whichParameters}{possibility to select parameters by index} @@ -168,9 +168,9 @@ getSample( Extracts the sample from a bayesianOutput } \details{ -If thin is greater than the total number of samples in the sampler object the first and the last element (of each chain if a sampler with multiples chains is used) are sampled. If numSamples is greater than the total number of samples all samples are selected. In both cases a warning is displayed. +If thin is greater than the total number of samples in the sampler object, the first and the last element (of each chain if a sampler with multiples chains is used) are sampled. If numSamples is greater than the total number of samples all samples are selected. A warning will be displayed in both cases. -If thin and numSamples is passed, the function will use the thin argument if it is valid and greater than 1, else numSamples will be used. +If both thin and numSamples are provided, the function will use thin only if it is valid and greater than 1; otherwise, numSamples will be used. } \examples{ ll = function(x) sum(dnorm(x, log = TRUE)) diff --git a/BayesianTools/man/getVolume.Rd b/BayesianTools/man/getVolume.Rd index 9d63ed9..55e4659 100644 --- a/BayesianTools/man/getVolume.Rd +++ b/BayesianTools/man/getVolume.Rd @@ -7,9 +7,9 @@ getVolume(sampler, prior = F, method = "MVN", ...) } \arguments{ -\item{sampler}{an object of superclass bayesianOutput or any other class that has the getSample function implemented (e.g. Matrix)} +\item{sampler}{an object of superclass bayesianOutput or any other class that has implemented the getSample function (e.g. Matrix)} -\item{prior}{schould also prior volume be calculated} +\item{prior}{logical, should prior volume be calculated?} \item{method}{method for volume estimation. Currently, the only option is "MVN"} @@ -19,7 +19,7 @@ getVolume(sampler, prior = F, method = "MVN", ...) Calculate posterior volume } \details{ -The idea of this function is to provide an estimate of the "posterior volume", i.e. how "broad" the posterior is. One potential application is to the overall reduction of parametric uncertainty between different data types, or between prior and posterior. +The idea of this function is to provide an estimate of the "posterior volume", i.e. how "broad" the posterior is. One potential application is the overall reduction of parametric uncertainty between different data types, or between prior and posterior. Implemented methods for volume estimation: diff --git a/BayesianTools/man/makeObjectClassCodaMCMC.Rd b/BayesianTools/man/makeObjectClassCodaMCMC.Rd index ccabf2c..d442a60 100644 --- a/BayesianTools/man/makeObjectClassCodaMCMC.Rd +++ b/BayesianTools/man/makeObjectClassCodaMCMC.Rd @@ -9,14 +9,14 @@ makeObjectClassCodaMCMC(chain, start = 1, end = numeric(0), thin = 1) \arguments{ \item{chain}{mcmc Chain} -\item{start}{for mcmc samplers start value in the chain. For SMC samplers, start particle} +\item{start}{For MCMC samplers, the initial value in the chain. For SMC samplers, initial particle} -\item{end}{for mcmc samplers end value in the chain. For SMC samplers, end particle} +\item{end}{For MCMC samplers, the end value in the chain. For SMC samplers, end particle.} \item{thin}{thinning parameter} } \value{ -object of class coda::mcmc +object an object of class coda::mcmc } \description{ Helper function to change an object to a coda mcmc class, diff --git a/BayesianTools/man/marginalLikelihood.Rd b/BayesianTools/man/marginalLikelihood.Rd index dabc68e..684d15e 100644 --- a/BayesianTools/man/marginalLikelihood.Rd +++ b/BayesianTools/man/marginalLikelihood.Rd @@ -16,7 +16,7 @@ marginalLikelihood(sampler, numSamples = 1000, method = "Chib", ...) \item{...}{further arguments passed to \code{\link{getSample}}} } \value{ -A list with log of the marginal likelihood, as well as other diagnostics depending on the chose method +A list with log of the marginal likelihood, as well as other diagnostics depending on the chosen method } \description{ Calcluated the marginal likelihood from a set of MCMC samples @@ -32,11 +32,11 @@ Given that MLs are calculated for each model, you can get posterior weights (for In BT, we return the log ML, so you will have to exp all values for this formula. -It is well-known that the ML is VERY dependent on the prior, and in particular the choice of the width of uninformative priors may have major impacts on the relative weights of the models. It has therefore been suggested to not use the ML for model averaging / selection on uninformative priors. If you have no informative priors, and option is to split the data into two parts, use one part to generate informative priors for the model, and the second part for the model selection. See help for an example. +It is well-known that the ML is strongly dependent on the prior, and in particular the choice of the width of uninformative priors may have major impacts on the relative weights of the models. It has therefore been suggested to not use the ML for model averaging / selection on uninformative priors. If you have no informative priors, and option is to split the data into two parts, use one part to generate informative priors for the model, and the second part for the model selection. See help for an example. The marginalLikelihood function currently implements four ways to calculate the marginal likelihood. Be aware that marginal likelihood calculations are notoriously prone to numerical stability issues. Especially in high-dimensional parameter spaces, there is no guarantee that any of the implemented algorithms will converge reasonably fast. The recommended (and default) method is the method "Chib" (Chib and Jeliazkov, 2001), which is based on MCMC samples, with a limited number of additional calculations. Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions. -The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unrealiable and usually should not be used. +The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unreliable and usually should not be used. The third method is simply sampling from the prior. While in principle unbiased, it will only converge for a large number of samples, and is therefore numerically inefficient. @@ -121,7 +121,7 @@ exp(M1$ln.ML - M2$ln.ML) # it has therefore been suggested that ML should not be calculated on uninformative priors. But # what to do if there are no informative priors? -# one option is to calculate the fractional BF, which means that one splites the data in half, +# one option is to calculate the fractional BF, which means that one splits the data in half, # uses the first half to fit the model, and then use the posterior as a new (now informative) # prior for the ML - let's do this for the previous case @@ -190,7 +190,7 @@ exp(M1$ln.ML - M2$ln.ML) # Low dimensional case with narrow priors - all methods have low error # we use a truncated normal for the likelihood to make sure that the density -# integrates to 1 - makes it easier to calcuate the theoretical ML +# integrates to 1 - makes it easier to calculate the theoretical ML likelihood <- function(x) sum(msm::dtnorm(x, log = TRUE, lower = -1, upper = 1)) prior = createUniformPrior(lower = rep(-1,2), upper = rep(1,2)) bayesianSetup <- createBayesianSetup(likelihood = likelihood, prior = prior) @@ -225,7 +225,6 @@ marginalLikelihood(out, method = "HM", numSamples = 500)$ln.ML - theory marginalLikelihood(out, method = "Bridge", numSamples = 500)$ln.ML - theory - } \references{ Chib, Siddhartha, and Ivan Jeliazkov. "Marginal likelihood from the Metropolis-Hastings output." Journal of the American Statistical Association 96.453 (2001): 270-281. diff --git a/BayesianTools/man/mcmcMultipleChains.Rd b/BayesianTools/man/mcmcMultipleChains.Rd index 8cbdcec..e97e31b 100644 --- a/BayesianTools/man/mcmcMultipleChains.Rd +++ b/BayesianTools/man/mcmcMultipleChains.Rd @@ -7,9 +7,9 @@ mcmcMultipleChains(bayesianSetup, settings, sampler) } \arguments{ -\item{bayesianSetup}{Object of class "BayesianSetup"} +\item{bayesianSetup}{object of class "BayesianSetup"} -\item{settings}{list with settings for sampler} +\item{settings}{list containing settings for sampler} \item{sampler}{character, either "Metropolis" or "DE"} } diff --git a/BayesianTools/man/package-deprecated.Rd b/BayesianTools/man/package-deprecated.Rd index 5333bc4..84af775 100644 --- a/BayesianTools/man/package-deprecated.Rd +++ b/BayesianTools/man/package-deprecated.Rd @@ -9,7 +9,7 @@ createMixWithDefaults(pars, defaults, locations) \arguments{ \item{pars}{vector with new parameter values} -\item{defaults}{vector with defaukt parameter values} +\item{defaults}{vector with default parameter values} \item{locations}{indices of the new parameter values} } diff --git a/BayesianTools/man/plotSensitivity.Rd b/BayesianTools/man/plotSensitivity.Rd index 4546ed4..6f55c63 100644 --- a/BayesianTools/man/plotSensitivity.Rd +++ b/BayesianTools/man/plotSensitivity.Rd @@ -20,11 +20,10 @@ Performs a one-factor-at-a-time sensitivity analysis for the posterior of a give This function can also be used for sensitivity analysis of an arbitrary output - just create a BayesianSetup with this output. } \examples{ - ll <- testDensityBanana bayesianSetup <- createBayesianSetup(likelihood = ll, lower = rep(-10, 2), upper = rep(10, 2)) - plotSensitivity(bayesianSetup) + } \author{ Florian Hartig diff --git a/BayesianTools/man/runMCMC.Rd b/BayesianTools/man/runMCMC.Rd index 16f4f02..eacee0f 100644 --- a/BayesianTools/man/runMCMC.Rd +++ b/BayesianTools/man/runMCMC.Rd @@ -44,7 +44,7 @@ The MCMC samplers will have a number of additional settings, which are described Note that even if you specify parallel = T, this will only turn on internal parallelization of the samplers. The independent samplers controlled by nrChains are not evaluated in parallel, so if time is an issue it will be better to run the MCMCs individually and then combine them via \code{\link{createMcmcSamplerList}} into one joint object. -Note that DE and DREAM variants as well as SMC and T-walk require a population to start, which should be provided as a matrix. Default (NULL) sets the population size for DE to 3 x dimensions of parameters, for DREAM to 2 x dimensions of parameters and for DEzs and DREAMzs to three, sampled from the prior. Note also that the zs variants of DE and DREAM require two populations, the current population and the z matrix (a kind of memory) - if you want to set both, provide a list with startvalue$X and startvalue$Z. +Note that, DE and DREAM variants as well as SMC and T-walk require a population to start, which should be provided as a matrix. Default (NULL) sets the population size for DE to 3 x dimensions of parameters, for DREAM to 2 x dimensions of parameters and for DEzs and DREAMzs to three, sampled from the prior. Note also that the zs variants of DE and DREAM require two populations, the current population and the z matrix (a kind of memory) - if you want to set both, provide a list with startvalue$X and startvalue$Z. setting startValue for sampling with nrChains > 1 : if you want to provide different start values for the different chains, provide them as a list } diff --git a/BayesianTools/man/stopParallel.Rd b/BayesianTools/man/stopParallel.Rd index 54ca80b..b898737 100644 --- a/BayesianTools/man/stopParallel.Rd +++ b/BayesianTools/man/stopParallel.Rd @@ -11,7 +11,7 @@ stopParallel(bayesianSetup) } \description{ Function closes -the parallel executer (if available) +the parallel executor (if available) } \author{ Stefan Paul diff --git a/BayesianTools/man/testLinearModel.Rd b/BayesianTools/man/testLinearModel.Rd index c2c9563..00a0b73 100644 --- a/BayesianTools/man/testLinearModel.Rd +++ b/BayesianTools/man/testLinearModel.Rd @@ -15,6 +15,7 @@ testLinearModel(x, env = NULL) Fake model, returns a ax + b linear response to 2-param vector } \examples{ + x = c(1,2) y = testLinearModel(x) plot(y) diff --git a/BayesianTools/man/tracePlot.Rd b/BayesianTools/man/tracePlot.Rd index 50afc14..92e15f6 100644 --- a/BayesianTools/man/tracePlot.Rd +++ b/BayesianTools/man/tracePlot.Rd @@ -17,6 +17,8 @@ tracePlot(sampler, thin = "auto", ...) Trace plot for MCMC class } \examples{ + + # set up and run the MCMC ll <- function(x) sum(dnorm(x, log = TRUE)) setup <- createBayesianSetup(likelihood = ll, lower = c(-10, -10), upper = c(10,10)) diff --git a/BayesianTools/man/updateGroups.Rd b/BayesianTools/man/updateGroups.Rd index 2067ea3..188180a 100644 --- a/BayesianTools/man/updateGroups.Rd +++ b/BayesianTools/man/updateGroups.Rd @@ -7,9 +7,9 @@ updateGroups(chain, blockSettings) } \arguments{ -\item{chain}{MCMC chain including only the parameters (not logP,ll, logP)} +\item{chain}{MCMC chain including only the parameters (not logP, ll, logP)} -\item{blockSettings}{list with settings} +\item{blockSettings}{a list with settings} } \value{ groups diff --git a/BayesianTools/vignettes/BayesianTools.Rmd b/BayesianTools/vignettes/BayesianTools.Rmd index 84ee809..8ad9f5d 100644 --- a/BayesianTools/vignettes/BayesianTools.Rmd +++ b/BayesianTools/vignettes/BayesianTools.Rmd @@ -34,9 +34,7 @@ If you haven't installed the package yet, either run install.packages("BayesianTools") ``` -or follow the instructions at - to install a -development version or an older version. +Or follow the instructions on to install a development or an older version. Loading and citation @@ -80,15 +78,13 @@ ll <- generateTestDensityMultiNormal(sigma = "no correlation") bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper = rep(10, 3)) ``` -A more detailed description of the `BayesianSetup` will follow. +A more detailed description of the `BayesianSetup` will follow below. **Hint:** For an example of how to run these steps for a dynamic ecological model, see `?VSEM`. ## Run MCMC and SMC functions -After setup, you may need to run a calibration. The `runMCMC` function -serves as the main wrapper for all other implemented `MCMC`/`SMC` -functions. It always requires the following arguments: +After setup, you may need to run a calibration. The `runMCMC` function serves as the main wrapper for all other implemented `MCMC`/`SMC` functions. It always requires the following arguments: - `bayesianSetup` (alternatively, the log target function) - `sampler` name @@ -165,6 +161,35 @@ getSample(out, start = 100, end = NULL, thin = 5, whichParameters = 1:2) ### Which sampler to choose? +```{r, echo = T} +iter = 1000 +settings = list(iterations = iter, nrChains = 3, message = FALSE) +out <- runMCMC(bayesianSetup = bayesianSetup, sampler = "Metropolis", settings = settings) + +``` + +The result is an object of mcmcSamplerList, which should allow to do everything one can do with an mcmcSampler object (with slightly different output sometimes). + +```{r} +print(out) +summary(out) +``` + +For example, in the plot you now see 3 chains. + +```{r} +plot(out) +``` + +There are a few additional functions that may only be available for lists, for example convergence checks + +```{r} +#getSample(out, coda = F) +gelmanDiagnostics(out, plot = T) +``` + +#### Which sampler to choose? + The BT package provides a large class of different MCMC samplers, and it depends on the particular application which one is most suitable. In the absence of further information, we currently recommend the `DEz`s sampler. This is also the default in the `runMCMC` function. @@ -191,16 +216,8 @@ bayesianSetup = createBayesianSetup(likelihood = ll, lower = rep(-10, 3), upper Likelihoods are often costly to compute. If this is the case for you, you should think about parallelization possibilities. The 'createBayesianSetup' function has an input variable 'parallel', with the following options - F / FALSE means no parallelization should be used -- T / TRUE means to use R's automatic parallelization options (note: - this will not work if your likelihood is writing to file, or using - global variables or functions - see R's general help on - parallelization) -- "external", assumes that the likelihood is already parallelized. In - this case, the function needs to take a matrix with parameters as - columns, and rows as the different model runs you want to evaluate. - This is the most likely option to use if you have a complicated - setup (file I/O, HPC cluster) that cannot be handled with the - standard R parallelization. +- T / TRUE means to use R's automatic parallelization options (note: this will not work if your likelihood is writing to file, or using global variables or functions - see R's general help on parallelization) +- "external", assumes that the likelihood is already parallelized. In this case, the function needs to take a matrix with parameters as columns, and rows as the different model runs you want to evaluate. This is the most likely option to use if you have a complicated setup (file I/O, HPC cluster) that cannot be handled with the standard R parallelization. Algorithms in the `BayesianTools` package can take advantage of parallelization if this option is specified in the `BayesianSetup`. Note that parallelization is currently used by the following algorithms: `SMC`, `DEzs` and `DREAMzs` sampler. It can also be used through the `BayesianSetup` with the functions of the sensitivity package. @@ -337,11 +354,8 @@ This information can be passed by first creating an extra object, via `createPri You have 5 options to create a prior - Do not set a prior - in this case, an infinite prior will be created -- Set min/max values - a bounded flat prior and the corresponding - sampling function will be created -- Use one of the predefinded priors, see `?createPrior` for a list. - One of the options here is to use a previous MCMC output as the new - prior. Predefined priors will usually come with a sampling function +- Set min/max values - a bounded flat prior and the corresponding sampling function will be created +- Use one of the predefinded priors, see `?createPrior` for a list. One of the options here is to use a previous MCMC output as the new prior. Predefined priors will usually come with a sampling function - Use a user-defined prior, see `?createPrior` - Create a prior from a previous MCMC sample @@ -349,15 +363,10 @@ You have 5 options to create a prior When creating a user-defined prior, the following information can/should be passed to `createPrior` -- A log density function, as a function of a parameter vector x, same - syntax as the likelihood. -- In addition, you should consider providing a function that samples - from the prior, as many samplers (SMC, DE, DREAM) can use this - function for initial conditions. If you use one of the predefined - priors, the sampling function is already implemented. -- Lower / upper bounds (can be set on top of any prior to create - truncation) -- Additional information - best values, parameter names, ... +- A log density function, as a function of a parameter vector x, same syntax as the likelihood +- Additionally, you should consider providing a function that samples from the prior, because many samplers (SMC, DE, DREAM) can make use of this function for initial conditions. If you use one of the pre-defined priors, the sampling function is already implemented +- lower / upper boundaries (can be set on top of any prior, to create truncation) +- Additional info - best values, names of the parameters, ... #### Creating a prior from a previous MCMC sample @@ -442,8 +451,7 @@ The subsequent code provides an overview of the default settings of the MH `samp applySettingsDefault(sampler = "Metropolis") ``` -Here are some examples of how to apply different settings. Activate -individual options or combinations as demonstrated. +Here are some examples of how to apply different settings. Activate individual options or combinations as demonstrated. ### Standard MH MCMC @@ -596,6 +604,18 @@ An alternative to MCMCs are particle filters, also known as Sequential Monte-Car ### Rejection sampling +#chain = getSample(out, coda = T) +gelmanDiagnostics(out, plot = F) +``` + +## Non-MCMC sampling algorithms + +MCMCs sample the posterior space by creating a chain in parameter space. While this allows "learning" from past steps, it does not permit the parallel execution of a large number of posterior values at the same time. + +An alternative to MCMCs are particle filters, aka Sequential Monte-Carlo (SMC) algorithms. See Hartig, F.; Calabrese, J. M.; Reineking, B.; Wiegand, T. & Huth, A. Statistical inference for stochastic simulation models - theory and application Ecol. Lett., 2011, 14, 816-827 + +### Rejection sampling + The simplest option is to sample a large number of parameters and accept them according to their posterior value. This option can be emulated with the implemented SMC by setting iterations to 1. ```{r, results = 'hide', eval = F} @@ -622,13 +642,13 @@ There are a number of Bayesian methods for model selection and model comparison. - On the Bayes factor, see Kass, R. E. & Raftery, A. E. Bayes Factors J. Am. Stat. Assoc., Amer Statist Assn, 1995, 90, 773-795 -- For an overview of DIC and WAIC, see Gelman, A.; Hwang, J. & Vehtari, A. (2014) Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24, 997-1016-. On DIC, see also the original reference by Spiegelhalter, D. J.; Best, N. G.; Carlin, B. P. & van der Linde, A. (2002) Bayesian measures of model complexity and fit. J. Roy. Stat. Soc. B, 64, 583-639. +- For an overview of DIC and WAIC, see Gelman, A.; Hwang, J. & Vehtari, A. (2014) Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24, 997-1016-. On DIC, see also the original reference by Spiegelhalter, D. J.; Best, N. G.; Carlin, B. P. & van der Linde, A. (2002) Bayesian measures of model complexity and fit. J. Roy. Stat. Soc. B, 64, 583-639. The Bayes factor relies on the estimation of marginal likelihoods, which is numerically challenging. The BT package currently implements three methods -- The recommended method is the "Chib" method (Chib and Jeliazkov, 2001), which is based on MCMC samples but performs additional calculation. Although this is the current recommendation, note that there are some numerical issues with this algorithm that may limit its reliability for larger dimensions. +- The recommended method is the "Chib" method (Chib and Jeliazkov, 2001), which is based on MCMC samples but performs additional calculation. Although this is the current recommendation, note that there are some numerical issues with this algorithm that may limit its reliability for larger dimensions. -- The harmonic mean approximation is implemented for comparison purposes only. Note that this method is numerically unreliable and should not normally be used. +- The harmonic mean approximation, is implemented only for comparison. Note that the method is numerically unrealiable and usually should not be used. - The third method is simply sampling from the prior. While in principle unbiased, it converges only for a large number of samples and is therefore numerically inefficient. @@ -699,9 +719,7 @@ Assuming equal prior weights for all models, we can calculate the posterior weig exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML)) ``` -If the models have different model priors, then multiply by the prior probability of each of the models. - -### Comparing Models by DIC +If models have different model priors, multiply with the prior probabilities of each model. The Deviance Information Criterion is a commonly used method to summarize the fit of an MCMC chain. It can be calculated using