Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update s3 class - fixed problems in R CMD check #261

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
2 changes: 1 addition & 1 deletion BayesianTools/DESCRIPTION
Original file line number Diff line number Diff line change
@@ -40,7 +40,7 @@ Suggests:
testthat
LinkingTo: Rcpp
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.2.1
RoxygenNote: 7.2.3
URL: https://github.com/florianhartig/BayesianTools
BugReports: https://github.com/florianhartig/BayesianTools/issues
VignetteBuilder: knitr
1 change: 0 additions & 1 deletion BayesianTools/NAMESPACE
Original file line number Diff line number Diff line change
@@ -27,7 +27,6 @@ S3method(print,smcSamplerList)
S3method(summary,mcmcSampler)
S3method(summary,mcmcSamplerList)
S3method(summary,smcSampler)
S3method(summary,smcSamplerList)
export(DE)
export(DEzs)
export(DIC)
17 changes: 13 additions & 4 deletions BayesianTools/R/classBayesianOutput.R
Original file line number Diff line number Diff line change
@@ -1,21 +1,27 @@
# 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
#' Extracts samples from a bayesianOutput
#' @description
#' This function extracts samples from a BT MCMC or SMC sampler after they have been run. It can also work on objects, for example coda::mcmc, matrix,
#'
#' It can return either a matrix or a coda mcmc object.
#'
#' @param sampler an object, typically of class mcmcSampler, mcmcSamplerList (created by \code{\link{runMCMC}}), 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 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 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.
#' @return a matrix or a coda object, depending on option
#' @author Florian Hartig
#' @export
getSample <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics = FALSE, ...) UseMethod("getSample")

@@ -219,6 +225,9 @@ getSample.MCMC <- function(sampler, parametersOnly = T, coda = F, start = 1, end
return(getSample(as.matrix(sampler$mvSamples), parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics = reportDiagnostics))
}


# getSample implementation for objects created by Bridge Sampling???

#' @rdname getSample
#' @author Tankred Ott
#' @export
20 changes: 17 additions & 3 deletions BayesianTools/R/classBayesianSetup.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@

#' Creates a standardized collection of prior, likelihood and posterior functions, including error checks etc.
#' @author Florian Hartig, Tankred Ott
#' @param likelihood log likelihood density function
#' @param prior either a prior class (see \code{\link{createPrior}}) or a log prior density function
#' @param likelihood either an object of class likelihood created by \code{\link{createLikelihood}} or a specialized likelihood function, or a log likelihood density function. If you provide a function, parallel, paralleOptions and catchDuplicates will be internally used to create the likelihood object
#' @param prior either an object of a prior created by \code{\link{createPrior}} or a log prior density function
#' @param priorSampler if a prior density (and not a prior class) is provided to prior, the optional prior sampling function can be provided here
#' @param lower vector with lower prior limits
#' @param upper vector with upper prior limits
@@ -14,6 +14,9 @@
#' @param plotLower vector with lower limits for plotting
#' @param plotUpper vector with upper limits for plotting
#' @param plotBest vector with best values for plotting
#'
#' @returns An object of class BayesianSetup. Implemented S3 generics include \code{\link{print.BayesianSetup}}.
#'
#' @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.
@@ -28,6 +31,7 @@
#' @seealso \code{\link{checkBayesianSetup}} \cr
#' \code{\link{createLikelihood}} \cr
#' \code{\link{createPrior}} \cr
#'
#' @example /inst/examples/classBayesianSetup.R
#'
#'
@@ -185,9 +189,19 @@ stopParallel <- function(bayesianSetup){
}


#' Print an object of BayesianSetup
#'
#' @description print function for objects of class BayesianSetup
#'
#'
#' @param x an object of class BayesianSetup
#' @param ... not implemented
#'
#' @author Maximilian Pichler
#'
#' @export

#' @seealso \code{\link{checkBayesianSetup}} \cr
#' \code{\link{createBayesianSetup}} \cr
print.BayesianSetup <- function(x, ...){
cat('BayesianSetup: \n\n')

11 changes: 6 additions & 5 deletions BayesianTools/R/classLikelihood.R
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
#' Creates a standardized likelihood class#'
#' Creates a standardized likelihood class
#'
#' @author Florian Hartig
#' @param likelihood Log likelihood density
#' @param names Parameter names (optional)
#' @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 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 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
29 changes: 25 additions & 4 deletions BayesianTools/R/classMcmcSampler.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# Functions for class mcmcSamper

# Functions for class mcmcSampler
#' @rdname getSample
#' @author Florian Hartig
#' @export
getSample.mcmcSampler <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics= F, ...){
@@ -113,10 +113,18 @@ getSample.mcmcSampler <- function(sampler, parametersOnly = T, coda = F, start =



#' Summmary of MCMC output
#' @description
#' Creates a summary table of a MCMC output
#' @param object object of class mcmcSampler or mcmcSamplerList
#' @param printCorrelation if set to TRUE, prints correlation among samples
#' @param ... not implemented
#' @method summary mcmcSampler
#' @author Stefan Paul
#' @export
#' @seealso \code{\link{getSample.mcmcSampler}}
summary.mcmcSampler <- function(object, printCorrelation = "auto", ...){

#codaChain = getSample(sampler, parametersOnly = parametersOnly, coda = T, ...)
#summary(codaChain)
#rejectionRate(sampler$codaChain)
@@ -209,9 +217,14 @@ summary.mcmcSampler <- function(object, printCorrelation = "auto", ...){
}
}

#' @author Florian Hartig
#' Prints MCMC output
#' @description
#' Prints MCMC output
#' @param x object of class mcmcSampler or mcmcSamplerList
#' @param ... additional options
#' @method print mcmcSampler
#' @export
#' @seealso \code{\link{getSample.mcmcSampler}}
print.mcmcSampler <- function(x, ...){
print("mcmcSampler - you can use the following methods to summarize, plot or reduce this class:")
print(methods(class ="mcmcSampler"))
@@ -220,9 +233,17 @@ print.mcmcSampler <- function(x, ...){
#effectiveSize(sampler$codaChain)
}

#' @author Florian Hartig


#' Plots of MCMC output
#' @description
#' Plots MCMC output
#' @param x object of class mcmcSampler or mcmcSamplerList
#' @param ... additional options passed to tracePlot
#' @method plot mcmcSampler
#' @author Florian Hartig
#' @export
#' @seealso \code{\link{getSample.mcmcSampler}}
plot.mcmcSampler <- function(x, ...){
tracePlot(x, ...)
}
8 changes: 7 additions & 1 deletion BayesianTools/R/classMcmcSamplerList.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
#' 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
#' @param mcmcList a list of objects of class 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.
for (i in 1:length(mcmcList)){
@@ -14,6 +15,7 @@ createMcmcSamplerList <- function(mcmcList){

#' @author Stefan Paul
#' @method summary mcmcSamplerList
#' @describeIn summary.mcmcSampler
#' @export
summary.mcmcSamplerList <- function(object, ...){
#codaChain = getSample(sampler, parametersOnly = parametersOnly, coda = T, ...)
@@ -100,8 +102,10 @@ summary.mcmcSamplerList <- function(object, ...){

}


#' @author Florian Hartig
#' @method print mcmcSamplerList
#' @describeIn print.mcmcSampler
#' @export
print.mcmcSamplerList <- function(x, ...){
print("mcmcSamplerList - you can use the following methods to summarize, plot or reduce this class:")
@@ -112,13 +116,15 @@ print.mcmcSamplerList <- function(x, ...){
}

#' @method plot mcmcSamplerList
#' @describeIn plot.mcmcSampler
#' @export
plot.mcmcSamplerList <- function(x, ...){
tracePlot(x, ...)
}

#' @author Florian Hartig
#' @export
#' @rdname getSample
getSample.mcmcSamplerList <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics, ...){

if(!is.null(numSamples)) numSamples = ceiling(numSamples/length(sampler))
2 changes: 1 addition & 1 deletion BayesianTools/R/classPosterior.R
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#' Creates a standardized posterior class
#' @author Florian Hartig
#' @param prior prior class
#' @param likelihood Log likelihood density
#' @param likelihood log likelihood density
#' @details Function is internally used in \code{\link{createBayesianSetup}} to create a standarized posterior class.
#' @export
createPosterior <- function(prior, likelihood){
8 changes: 6 additions & 2 deletions BayesianTools/R/classPrior.R
Original file line number Diff line number Diff line change
@@ -13,6 +13,7 @@
#' \code{\link{createUniformPrior}} \cr
#' \code{\link{createTruncatedNormalPrior}}\cr
#' \code{\link{createBayesianSetup}}\cr
#' \code{\link{print.prior}}\cr
#' @example /inst/examples/createPrior.R
createPrior <- function(density = NULL, sampler = NULL, lower = NULL, upper = NULL, best = NULL){

@@ -246,10 +247,13 @@ createPriorDensity <- function(sampler, method = "multivariate", eps = 1e-10, lo



#' Print object of prior class
#' @description print function for objects of class prior
#' @param x object of a prior
#' @param ... additional parameters
#' @author Maximilian Pichler

#' @export

#' @seealso \code{\link{createPrior}}
print.prior <- function(x, ...){
cat('Prior: \n\n')

4 changes: 3 additions & 1 deletion BayesianTools/R/classSMCSamplerList.R
Original file line number Diff line number Diff line change
@@ -15,21 +15,23 @@ createSmcSamplerList <- function(...){

#' @method summary smcSamplerList
#' @author Florian Hartig
#' @export
#' @describeIn summary.smcSampler
summary.smcSamplerList <- function(object, ...){
sample = getSample(object, parametersOnly = T, ...)
summary(sample)
}

#' @method print smcSamplerList
#' @author Florian Hartig
#' @describeIn print.smcSampler
#' @export
print.smcSamplerList <- function(x, ...){
print("smcSamplerList - you can use the following methods to summarize, plot or reduce this class:")
print(methods(class ="smcSamplerList"))
}

#' @method plot smcSamplerList
#' @describeIn plot.smcSampler
#' @export
plot.smcSamplerList <- function(x, ...){
marginalPlot(x, ...)
20 changes: 20 additions & 0 deletions BayesianTools/R/classSmcSampler.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,8 @@
# S3 Functions for class 'smcSampler'



#' @rdname getSample
#' @author Florian Hartig
#' @export
getSample.smcSampler <- function(sampler, parametersOnly = T, coda = F, start = 1, end = NULL, thin = 1, numSamples = NULL, whichParameters = NULL, reportDiagnostics = FALSE, ...){
@@ -43,6 +48,11 @@ getSample.smcSampler <- function(sampler, parametersOnly = T, coda = F, start =
} else return(out)
}

#' Summary for class 'smcSampler'
#' @description
#' Creates a summary table of a 'smcSampler' output
#' @param object object of class smcSampler or smcSamplerList
#' @param ... additional parameters
#' @author Florian Hartig
#' @method summary smcSampler
#' @export
@@ -52,12 +62,22 @@ summary.smcSampler<- function(object, ...){
summary(getSample(sampler, ...))
}

#' Plots of smcSampler output
#' @description
#' Plots smcSampler output
#' @param x object of class mcmcSampler or mcmcSamplerList
#' @param ... additional options passed to tracePlot
#' @method plot smcSampler
#' @export
plot.smcSampler<- function(x, ...){
marginalPlot(x, ...)
}

#' Print of smcSampler output
#' @description
#' Print smcSampler output
#' @param x object of class mcmcSampler or mcmcSamplerList
#' @param ... additional options
#' @author Florian Hartig
#' @method print smcSampler
#' @export
9 changes: 4 additions & 5 deletions BayesianTools/R/marginalLikelihood.R
Original file line number Diff line number Diff line change
@@ -1,14 +1,12 @@

# 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")
# @export
# marginalLikelihood <- function(x,lik,V,sampler$setup$likelihood$density,sampler$setup$prior$density,..., num.samples=1000,log=TRUE) UseMethod("marginalLikelihood")

#' Calcluated the marginal likelihood from a set of MCMC samples
#' @export
#' Calculated the marginal likelihood from a set of MCMC samples
#' @author Florian Hartig
#' @param sampler an MCMC or SMC sampler or list, or for method "Prior" also a BayesianSetup
#' @param numSamples number of samples to use. How this works, and if it requires recalculating the likelihood, depends on the method
@@ -44,6 +42,7 @@
#' 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}}
#' @export
marginalLikelihood <- function(sampler, numSamples = 1000, method = "Chib", ...){


4 changes: 3 additions & 1 deletion BayesianTools/R/mcmcRun.R
Original file line number Diff line number Diff line change
@@ -28,7 +28,9 @@
#'
#' setting startValue for sampling with nrChains > 1 : if you want to provide different start values for the different chains, provide them as a list
#'
#' @return The function returns an object of class mcmcSampler (if one chain is run) or mcmcSamplerList. Both have the superclass bayesianOutput. It is possible to extract the samples as a coda object or matrix with \code{\link{getSample}}.
#' @return The function returns an object of class mcmcSampler (if one chain is run) or mcmcSamplerList. Both have the superclass bayesianOutput. It is possible to extract the samples as a coda object or matrix with \code{\link{getSample}}. Other S3 classes that are implemented include \code{\link{summary.mcmcSampler}} or \code{\link{summary.mcmcSamplerList}}, \code{\link{print.mcmcSampler}} or \code{\link{print.mcmcSamplerList}}, \code{\link{plot.mcmcSampler}} or \code{\link{plot.mcmcSamplerList}}
#'
#'
#' It is also possible to summarize the posterior as a new prior via \code{\link{createPriorDensity}}.
#' @example /inst/examples/mcmcRun.R
#' @seealso \code{\link{createBayesianSetup}}
7 changes: 5 additions & 2 deletions BayesianTools/man/createBayesianSetup.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

12 changes: 6 additions & 6 deletions BayesianTools/man/createLikelihood.Rd
2 changes: 1 addition & 1 deletion BayesianTools/man/createMcmcSamplerList.Rd
2 changes: 1 addition & 1 deletion BayesianTools/man/createPosterior.Rd
1 change: 1 addition & 0 deletions BayesianTools/man/createPrior.Rd
58 changes: 53 additions & 5 deletions BayesianTools/man/getSample.Rd
4 changes: 2 additions & 2 deletions BayesianTools/man/marginalLikelihood.Rd
30 changes: 30 additions & 0 deletions BayesianTools/man/plot.mcmcSampler.Rd
24 changes: 24 additions & 0 deletions BayesianTools/man/plot.smcSampler.Rd
23 changes: 23 additions & 0 deletions BayesianTools/man/print.BayesianSetup.Rd
30 changes: 30 additions & 0 deletions BayesianTools/man/print.mcmcSampler.Rd
22 changes: 22 additions & 0 deletions BayesianTools/man/print.prior.Rd
27 changes: 27 additions & 0 deletions BayesianTools/man/print.smcSampler.Rd
3 changes: 2 additions & 1 deletion BayesianTools/man/runMCMC.Rd
32 changes: 32 additions & 0 deletions BayesianTools/man/summary.mcmcSampler.Rd
27 changes: 27 additions & 0 deletions BayesianTools/man/summary.smcSampler.Rd
98 changes: 98 additions & 0 deletions BayesianTools/vignettes/testVignette.Rmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
---
title: "Implementing BayesianTools on real data"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{Implementing BayesianTools on real data}
\usepackage[utf8]{inputenc}
%\VignetteEngine{knitr::rmarkdown}
abstract: "This tutorial discusses how to implement BayesianTools on real data"
author: Florian Hartig
editor_options:
chunk_output_type: console
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=5, fig.height=5, warning=FALSE, cache = F)
```

```{r, echo = F, message = F}
set.seed(123)
```

In this vignette, we want to illustrate how BayesianTools can be implemented on a real data.

In this case, we are using airquality data.

For the illustration purpose, we are fitting a linear regression model on the data - so that we can compare it with our result.


```{r}
fit = lm(Ozone ~ Wind, data = airquality)
summary(fit)
```

```{r}
dat = airquality[complete.cases(airquality), ]
```

Creating a likelihood function from the data.

```{r}
likelihood <- function(par){
# c("a", "b", "err-sd")
llObservation = sum(dnorm(par[1] * dat$Wind + par[2] - dat$Ozone, sd = par[3], log = T))
return(llObservation)
}
likelihood(c(4,10,6))
```

Creating a standardized collection of prior, likelihood and posterior function.

```{r, message=F}
library(BayesianTools)
setup <- createBayesianSetup(likelihood = likelihood,
lower = c(-100,-100,0),
upper = c(100,200,50),
names = c("a", "b", "err-sd"))
```

Creating required setup for wrapper function runMCMC to start MCMCs

```{r}
settings <- list(iterations = 10000, nrChains = 2)
res <- runMCMC(bayesianSetup = setup, settings = settings)
summary(res, start=1000)
```


```{r}
plot(res)
```

```{r}
marginalPlot(res, start = 1000, prior = T)
```

```{r}
gelmanDiagnostics(res)
```
```{r}
correlationPlot(res)
```

# Linear mixed model with the BT package





# Fitting a ecosystem model with the BT package