Merge pull request #173 from florianhartig/0.1.5.1

0.1.5.1

Author: florianhartig

BayesianTools/DESCRIPTION

@@ -1,7 +1,7 @@
 Package: BayesianTools
 Title: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics
-Version: 0.1.5
-Date: 2018-07-16
+Version: 0.1.5.1
+Date: 2019-01-20
 Authors@R: c(person("Florian", "Hartig", email = "florian.hartig@biologie.uni-regensburg.de", role = c("aut", "cre")), person("Francesco", "Minunno", role = c("aut")), person("Stefan", " Paul", role = c("aut") ),  person("David", "Cameron", role = "ctb"), person("Tankred", "Ott", role = "ctb"), person("Maximilian", "Pichler", role = "ctb"))
 Author: Florian Hartig [aut, cre],
     Francesco Minunno [aut],

BayesianTools/NEWS

@@ -1,5 +1,25 @@
 NOTE: for more news about the package, see https://github.com/florianhartig/BayesianTools/releases
 
+# BayesianTools 0.1.6
+
+New features
+
+- calibrationTest + helper functions to test for  https://github.com/florianhartig/BayesianTools/pull/172
+
+Major changes
+
+- None
+
+Minor changes
+
+- Help improvements and clarifications
+
+Bugfixes 
+
+- Bugfix AR1 likelihood https://github.com/florianhartig/BayesianTools/issues/169
+- Bugfix restart MCMC https://github.com/florianhartig/BayesianTools/issues/158
+
+
 # BayesianTools 0.1.5
 
 Changes / improvements

BayesianTools/R/marginalLikelihood.R

@@ -14,18 +14,35 @@
 #' @param numSamples number of samples to use. How this works, and if it requires recalculating the likelihood, depends on the method
 #' @param method method to choose. Currently available are "Chib" (default), the harmonic mean "HM", sampling from the prior "prior", and bridge sampling "Bridge". See details
 #' @param ... further arguments passed to \code{\link{getSample}}
-#' @details The function currently implements four ways to calculate the marginal likelihood. The recommended (and default) method is the method "Chib" (Chib and Jeliazkov, 2001), which is based on MCMC samples, but performs additional calculations. Despite being the current recommendation, note there are some numeric issues with this algorithm that may limit reliability for larger dimensions.
+#' @details The marginal likelihood is the average likelihood across the prior space. It is used, for example, for Bayesian model selection and model averaging. 
+#' 
+#' It is defined as \deqn{ML = \int L(\Theta) p(\Theta) d\Theta}
+#' 
+#' Given that MLs are calculated for each model, you can get posterior weights (for model selection and/or model averaging) on the model by 
+#' 
+#' \deqn{P(M_i|D) = ML_i * p(M_i) / (\sum_i ML_i * p(M_i) )}
+#' 
+#' 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 Dormann et al., 2018, in particular the Appendix, 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. \cr
+#'  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 will only converge for a large number of samples, and is therefore numerically inefficient. \cr
+#' 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". 
+#' 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.
 #' 
-#' @note 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.
+#' @return A list with log of the marginal likelihood, as well as other diagnostics depending on the chose method
 #'    
 #' @example /inst/examples/marginalLikelihoodHelp.R
-#' @references Chib, Siddhartha, and Ivan Jeliazkov. "Marginal likelihood from the Metropolis-Hastings output." Journal of the American Statistical Association 96.453 (2001): 270-281.
+#' @references 
+#' 
+#' Chib, Siddhartha, and Ivan Jeliazkov. "Marginal likelihood from the Metropolis-Hastings output." Journal of the American Statistical Association 96.453 (2001): 270-281.
+#' 
+#' Dormann et al. 2018. Model averaging in ecology: a review of Bayesian, information-theoretic, and tactical approaches for predictive inference. Ecological Monographs
+#' 
 #' @seealso \code{\link{WAIC}}, \code{\link{DIC}}, \code{\link{MAP}}
 marginalLikelihood <- function(sampler, numSamples = 1000, method = "Chib", ...){
 

BayesianTools/R/mcmcDREAM.R

@@ -7,10 +7,10 @@
 ##' @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 gamma Kurtosis parameter Bayesian Inference Scheme. Default is 0.
-##' @param eps Ergodicity term. Default to 0.
-##' @param e Ergodicity term. Default to 5e-2.
-##' @param pCRupdate Update of crossover probabilities, default is TRUE
+##' @param gamma Kurtosis parameter Bayesian Inference Scheme
+##' @param eps Ergodicity term
+##' @param e Ergodicity term
+##' @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).

BayesianTools/R/mcmcDREAMzs.R

@@ -7,10 +7,10 @@
 ##' @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. Default is 0.
-##' @param eps Ergodicity term. Default to 0.
-##' @param e Ergodicity term. Default to 5e-2.
-##' @param pCRupdate Update of crossover probabilities, default is TRUE
+##' @param gamma Kurtosis parameter Bayesian Inference Scheme. 
+##' @param eps Ergodicity term
+##' @param e Ergodicity term
+##' @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)

BayesianTools/inst/examples/marginalLikelihoodHelp.R

@@ -1,3 +1,55 @@
+
+# Comparison of ML for two regression models
+
+sampleSize = 30
+x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
+y <-  1 * x + 1*x^2 + rnorm(n=sampleSize,mean=0,sd=10)
+#plot(x,y, main="Test Data")
+
+
+# linear and quadratic effect
+likelihood1 <- function(param){
+  pred = param[1] + param[2]*x + param[3] * x^2
+  singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[4]^2), log = TRUE)
+  return(sum(singlelikelihoods))  
+}
+
+# linear  effect
+likelihood2 <- function(param){
+  pred = param[1] + param[2]*x 
+  singlelikelihoods = dnorm(y, mean = pred, sd = 1/(param[3]^2), log = TRUE)
+  return(sum(singlelikelihoods))  
+}
+
+setUp1 <- createBayesianSetup(likelihood1, lower = c(-5,-5,-5,0.01), upper = c(5,5,5,30))
+
+setUp2 <- createBayesianSetup(likelihood2, lower = c(-5,-5,0.01), upper = c(5,5,30))
+
+out1 <- runMCMC(bayesianSetup = setUp1)
+M1 = marginalLikelihood(out1, start = 1000)
+
+out2 <- runMCMC(bayesianSetup = setUp2)
+M2 = marginalLikelihood(out2, start = 1000)
+
+
+### Bayes factor
+
+exp(M1$ln.ML - M2$ln.ML)
+
+# BF > 1 means the evidence is in favor of M1. See Kass, R. E. & Raftery, A. E. 
+# (1995) Bayes Factors. J. Am. Stat. Assoc., Amer Statist Assn, 90, 773-795.
+
+### Posterior weight
+
+exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML))
+
+# If models have different model priors, multiply with the prior probabilities of each model. 
+
+############################################################
+
+### Performance comparison ### 
+
+
 # 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 

BayesianTools/man/DREAM.Rd

@@ -1,15 +1,13 @@
 ---
-title: "Bayesian Tools"
+title: "Bayesian Tools - General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics"
 output: 
   rmarkdown::html_vignette:
     toc: true
 vignette: >
   %\VignetteIndexEntry{Manual for the BayesianTools R package}
   %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteEncoding{UTF-8}
-abstract: "The BayesianTools (BT) package supports model analysis (including sensitivity analysis and uncertainty analysis), Bayesian model calibration, as well as model selection and multi-model inference techniques for system models. 
-\n \n \n"
-
+  \usepackage[utf8]{inputenc}
+abstract: "The BayesianTools (BT) package supports model analysis (including sensitivity analysis and uncertainty analysis), Bayesian model calibration, as well as model selection and multi-model inference techniques for system models." 
 ---
 
 ```{r global_options, include=FALSE}
@@ -364,7 +362,7 @@ out <- runMCMC(bayesianSetup = bayesianSetup, settings = settings)
 
 # MCMC sampler reference
 
-## runMCMC
+## The runMCMC() function
 
 The runMCMC function is the central function for starting MCMC algorithms in the BayesianTools package. It requires a bayesianSetup, a choice of sampler (standard is DEzs), and optionally changes to the standard settings of the chosen sampler.
 
@@ -723,12 +721,18 @@ M2
 Bayes factor (need to reverse the log)
 
 ```{r}
-exp(M1$marginalLikelihod - M2$marginalLikelihod)
+exp(M1$ln.ML - M2$ln.ML)
 ```
 
 BF > 1 means the evidence is in favor of M1. See Kass, R. E. & Raftery, A. E. (1995) Bayes Factors. J. Am. Stat. Assoc., Amer Statist Assn, 90, 773-795.
 
-Note that we would have to multiply still with the model priors to arrive at Bayesian model weights. 
+Assuming equal prior weights on all models, we can calculate the posterior weight of M1 as 
+
+```{r}
+exp(M1$ln.ML) / ( exp(M1$ln.ML) + exp(M2$ln.ML))
+```
+
+If models have different model priors, multiply with the prior probabilities of each model. 
 
 
 ### Model comparison via DIC