update

BayesianTools/R/MAP.R

@@ -1,20 +1,21 @@
-#' calculates the Maxiumum 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 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, ...){
-
- samples = getSample(bayesianOutput, parametersOnly = F, ...)
-
- if("mcmcSamplerList" %in% class(bayesianOutput)) nPars <- bayesianOutput[[1]]$setup$numPars
- else nPars = bayesianOutput$setup$numPars
-
- best = which.max(samples[,nPars + 1])
-
- return(list(parametersMAP = samples[best, 1:nPars], valuesMAP = samples[best, (nPars + 1):(nPars + 3)] ))
-
-}
+
+#' calculates the Maxiumum 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 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, ...){
+
+ samples = getSample(bayesianOutput, parametersOnly = F, ...)
+
+ if("mcmcSamplerList" %in% class(bayesianOutput)) nPars <- bayesianOutput[[1]]$setup$numPars
+ else nPars = bayesianOutput$setup$numPars
+
+ best = which.max(samples[,nPars + 1])
+
+ return(list(parametersMAP = samples[best, 1:nPars], valuesMAP = samples[best, (nPars + 1):(nPars + 3)] ))
+
+}
 

BayesianTools/R/classMcmcSamplerList.R

@@ -1,151 +1,152 @@
-#' Convenience function to create an object of class mcmcSamplerList from a list of mcmc samplers
-#' @author Florian Hartig
-#' @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.
- for (i in 1:length(mcmcList)){
- if (! ("mcmcSampler" %in% class(mcmcList[[i]])) ) stop("list objects are not of class mcmcSampler")
- }
- class(mcmcList) = c("mcmcSamplerList", "bayesianOutput")
- return(mcmcList)
-}
-
-#' @author Stefan Paul
-#' @method summary mcmcSamplerList
-#' @export
-summary.mcmcSamplerList <- function(object, ...){
- #codaChain = getSample(sampler, parametersOnly = parametersOnly, coda = T, ...)
- #summary(codaChain)
- #rejectionRate(sampler$codaChain)
- #effectiveSize(sampler$codaChain)
- #DIC(sampler)
- #max()
-
- sampler <- object
-
- DInf <- DIC(sampler)
- MAPvals <- round(MAP(sampler)$parametersMAP,3)
-
- gelDiag <- gelmanDiagnostics(sampler)
- psf <- round(gelDiag$psrf[,1], 3)
-
- mcmcsampler <- sampler[[1]]$settings$sampler
-
- runtime <- 0
- for(i in 1:length(sampler)) runtime <- runtime+sampler[[i]]$settings$runtime[3]
-
- correlations <- round(cor(getSample(sampler)),3)
-
-
- sampler <- getSample(sampler, parametersOnly = T, coda = T, ...)
- if("mcmc.list" %in% class(sampler)){
- nrChain <- length(sampler)
- nrIter <- nrow(sampler[[1]])
- conv <- round(gelDiag$mpsrf,3)
- npar <- ncol(sampler[[1]])
- lowerq <- upperq <- numeric(npar)
- medi <- numeric(npar)
- parnames <- colnames(sampler[[1]])
- for(i in 1:npar){
- tmp <- unlist(sampler[,i])
- tmp <- quantile(tmp, probs = c(0.025, 0.5, 0.975))
- lowerq[i] <- round(tmp[1],3)
- medi[i] <- round(tmp[2],3)
- upperq[i] <- round(tmp[3],3)
- }
-
- }else{
- nrChain <- 1
- nrIter <- nrow(sampler)
- npar <- ncol(sampler)
- conv <- "Only one chain; convergence cannot be determined!"
- medi <- numeric(npar)
- lowerq <- upperq <- numeric(npar)
- parnames <- colnames(sampler)
- for(i in 1:npar){
- tmp <- quantile(sampler[,i], probs = c(0.025, 0.5, 0.975))
- lowerq[i] <- round(tmp[1],3)
- medi[i] <- round(tmp[2],3)
- upperq[i] <- round(tmp[3],3)
- }
-
- }
-
- # output for parameter metrics
- parOutDF <- cbind(psf, MAPvals, lowerq, medi, upperq)
- colnames(parOutDF) <- c("psf", "MAP", "2.5%", "median", "97.5%")
- row.names(parOutDF) <- parnames
-
-
- cat(rep("#", 25), "\n")
- cat("## MCMC chain summary ##","\n")
- cat(rep("#", 25), "\n", "\n")
- cat("# MCMC sampler: ",mcmcsampler, "\n")
- cat("# Nr. Chains: ", nrChain, "\n")
- cat("# Iterations per chain: ", nrIter, "\n")
- cat("# Rejection rate: ", ifelse(object[[1]]$setup$numPars == 1, # this is a hack because coda::rejectionRate does not work for 1-d MCMC lists
- round(mean(sapply(sampler, coda::rejectionRate)),3),
- round(mean(coda::rejectionRate(sampler)),3) ), "\n")
- cat("# Effective sample size: ", round(mean(coda::effectiveSize(sampler)),0), "\n")
- cat("# Runtime: ", runtime, " sec.","\n", "\n")
- cat("# Parameters\n")
- print(parOutDF)
- cat("\n")
- cat("## DIC: ", round(DInf$DIC,3), "\n")
- cat("## Convergence" ,"\n", "Gelman Rubin multivariate psrf: ", conv, "\n","\n")
- cat("## Correlations", "\n")
- print(correlations)
-
-}
-
-#' @author Florian Hartig
-#' @method print mcmcSamplerList
-#' @export
-print.mcmcSamplerList <- function(x, ...){
- print("mcmcSamplerList - you can use the following methods to summarize, plot or reduce this class:")
- print(methods(class ="mcmcSamplerList"))
- #codaChain = getSample(sampler, coda = T, ...)
- #rejectionRate(sampler$codaChain)
- #effectiveSize(sampler$codaChain)
-}
-
-#' @method plot mcmcSamplerList
-#' @export
-plot.mcmcSamplerList <- function(x, ...){
- tracePlot(x, ...)
-}
-
-#' @author Florian Hartig
-#' @export
-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))
-
- if(coda == F){
- # out = NULL
- out <- rep(list(NA), length(sampler))
- for (i in 1:length(sampler)){
- # out = rbind(out, getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F))
- out[[i]] <- getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F)
- }
- out <- combineChains(out)
- }
-
- if(coda == T){
-
- out = list()
-
- for (i in 1:length(sampler)){
-
- out[[i]] = getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F)
- }
-
- if(inherits(out[[1]], "mcmc.list")) out = unlist(out, recursive = F)
- class(out) = "mcmc.list"
- out = out
- }
-
- return(out)
-}
+
+#' Convenience function to create an object of class mcmcSamplerList from a list of mcmc samplers
+#' @author Florian Hartig
+#' @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.
+ for (i in 1:length(mcmcList)){
+ if (! ("mcmcSampler" %in% class(mcmcList[[i]])) ) stop("list objects are not of class mcmcSampler")
+ }
+ class(mcmcList) = c("mcmcSamplerList", "bayesianOutput")
+ return(mcmcList)
+}
+
+#' @author Stefan Paul
+#' @method summary mcmcSamplerList
+#' @export
+summary.mcmcSamplerList <- function(object, ...){
+ #codaChain = getSample(sampler, parametersOnly = parametersOnly, coda = T, ...)
+ #summary(codaChain)
+ #rejectionRate(sampler$codaChain)
+ #effectiveSize(sampler$codaChain)
+ #DIC(sampler)
+ #max()
+
+ sampler <- object
+
+ DInf <- DIC(sampler)
+ MAPvals <- round(MAP(sampler)$parametersMAP,3)
+
+ gelDiag <- gelmanDiagnostics(sampler)
+ psf <- round(gelDiag$psrf[,1], 3)
+
+ mcmcsampler <- sampler[[1]]$settings$sampler
+
+ runtime <- 0
+ for(i in 1:length(sampler)) runtime <- runtime+sampler[[i]]$settings$runtime[3]
+
+ correlations <- round(cor(getSample(sampler)),3)
+
+
+ sampler <- getSample(sampler, parametersOnly = T, coda = T, ...)
+ if("mcmc.list" %in% class(sampler)){
+ nrChain <- length(sampler)
+ nrIter <- nrow(sampler[[1]])
+ conv <- round(gelDiag$mpsrf,3)
+ npar <- ncol(sampler[[1]])
+ lowerq <- upperq <- numeric(npar)
+ medi <- numeric(npar)
+ parnames <- colnames(sampler[[1]])
+ for(i in 1:npar){
+ tmp <- unlist(sampler[,i])
+ tmp <- quantile(tmp, probs = c(0.025, 0.5, 0.975))
+ lowerq[i] <- round(tmp[1],3)
+ medi[i] <- round(tmp[2],3)
+ upperq[i] <- round(tmp[3],3)
+ }
+
+ }else{
+ nrChain <- 1
+ nrIter <- nrow(sampler)
+ npar <- ncol(sampler)
+ conv <- "Only one chain; convergence cannot be determined!"
+ medi <- numeric(npar)
+ lowerq <- upperq <- numeric(npar)
+ parnames <- colnames(sampler)
+ for(i in 1:npar){
+ tmp <- quantile(sampler[,i], probs = c(0.025, 0.5, 0.975))
+ lowerq[i] <- round(tmp[1],3)
+ medi[i] <- round(tmp[2],3)
+ upperq[i] <- round(tmp[3],3)
+ }
+
+ }
+
+ # output for parameter metrics
+ parOutDF <- cbind(psf, MAPvals, lowerq, medi, upperq)
+ colnames(parOutDF) <- c("psf", "MAP", "2.5%", "median", "97.5%")
+ row.names(parOutDF) <- parnames
+
+
+ cat(rep("#", 25), "\n")
+ cat("## MCMC chain summary ##","\n")
+ cat(rep("#", 25), "\n", "\n")
+ cat("# MCMC sampler: ",mcmcsampler, "\n")
+ cat("# Nr. Chains: ", nrChain, "\n")
+ cat("# Iterations per chain: ", nrIter, "\n")
+ cat("# Rejection rate: ", ifelse(object[[1]]$setup$numPars == 1, # this is a hack because coda::rejectionRate does not work for 1-d MCMC lists
+ round(mean(sapply(sampler, coda::rejectionRate)),3),
+ round(mean(coda::rejectionRate(sampler)),3) ), "\n")
+ cat("# Effective sample size: ", round(mean(coda::effectiveSize(sampler)),0), "\n")
+ cat("# Runtime: ", runtime, " sec.","\n", "\n")
+ cat("# Parameters\n")
+ print(parOutDF)
+ cat("\n")
+ cat("## DIC: ", round(DInf$DIC,3), "\n")
+ cat("## Convergence" ,"\n", "Gelman Rubin multivariate psrf: ", conv, "\n","\n")
+ cat("## Correlations", "\n")
+ print(correlations)
+
+}
+
+#' @author Florian Hartig
+#' @method print mcmcSamplerList
+#' @export
+print.mcmcSamplerList <- function(x, ...){
+ print("mcmcSamplerList - you can use the following methods to summarize, plot or reduce this class:")
+ print(methods(class ="mcmcSamplerList"))
+ #codaChain = getSample(sampler, coda = T, ...)
+ #rejectionRate(sampler$codaChain)
+ #effectiveSize(sampler$codaChain)
+}
+
+#' @method plot mcmcSamplerList
+#' @export
+plot.mcmcSamplerList <- function(x, ...){
+ tracePlot(x, ...)
+}
+
+#' @author Florian Hartig
+#' @export
+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))
+
+ if(coda == F){
+ # out = NULL
+ out <- rep(list(NA), length(sampler))
+ for (i in 1:length(sampler)){
+ # out = rbind(out, getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F))
+ out[[i]] <- getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F)
+ }
+ out <- combineChains(out)
+ }
+
+ if(coda == T){
+
+ out = list()
+
+ for (i in 1:length(sampler)){
+
+ out[[i]] = getSample(sampler[[i]], parametersOnly = parametersOnly, coda = coda, start = start, end = end, thin = thin, numSamples = numSamples, whichParameters = whichParameters, reportDiagnostics= F)
+ }
+
+ if(inherits(out[[1]], "mcmc.list")) out = unlist(out, recursive = F)
+ class(out) = "mcmc.list"
+ out = out
+ }
+
+ return(out)
+}