update

BayesianTools/R/MAP.R

@@ -1,20 +1,20 @@
-#' 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 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 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/blockUpdate.R

@@ -1,90 +1,90 @@
-#' Determine the groups of correlated parameters
-#' @author Stefan Paul
-#' @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){
- 
- settings <- getBlockSettings(blockSettings)
- blockUpdateType <- settings$blockUpdateType
- 
- switch(blockUpdateType,
- "correlation" = {
- ## (Pair wise) Correlation in the parameters
- cormat <- abs(cor(chain[,1:(ncol(chain)-3),sample(1:dim(chain)[3],1)]))
- diag(cormat) <- 0
- # Correct for NA and Inf values as this could cause error in as.dist()
- cormat[c(which(is.na(cormat)),which(cormat == Inf),which(cormat == -Inf)) ] <- 0
- tree <- hclust(as.dist(1-cormat)) # get tree based on distance(dissimilarity = 1-cor).
- cT <- cutree(tree, k = settings$k, h = settings$h) # get groups. With h we can manipulate the strength of the interaction.
- },
- "user" = { 
- cT <- settings$groups
- },
- "random" = {
- pool <- c(1:settings$k, sample(1:settings$k, (ncol(chain)-3-settings$k)))
- cT <- sample(pool)
- }
- )
- 
- pSel <- settings$pSel
- if(is.null(pSel) && is.null(settings$pGroup)) pSel = rep(1,ncol(chain)-3)
- return(list(cT = cT, pGroup = settings$pGroup, pSel = pSel))
-}
-
-
-#' Determine the parameters in the block update
-#' @param blockSettings settings for block update
-#' @return vector containing the parameter to be updated
-#' @keywords internal
-getBlock <- function(blockSettings){
- groups <- blockSettings$cT
- pGroup <- blockSettings$pGroup
- pSel <- blockSettings$pSel
- 
- 
- nGroups = max(groups)
- if(nGroups == 1) return(1:length(groups))
- if (is.null(pGroup)) pGroup = rep(1,nGroups)
- if(length(pSel) > nGroups) pSel <- pSel[1:nGroups]
- pSel = c(pSel, rep(0,nGroups - length(pSel)))
- groupsToSample = sample.int(nGroups, 1, prob = pSel)
- 
- selectedGroups = sample.int(nGroups,groupsToSample, prob = pGroup[1:nGroups])
- GroupMember <- which(is.element(groups,selectedGroups))
- return(GroupMember)
- 
-}
-
-
-#' getblockSettings
-#' @description Transforms the original settings to settings used in the model runs
-#' @param blockUpdate input settings
-#' @return list with block settings
-#' @keywords internal
-getBlockSettings <- function(blockUpdate){
- 
- h <- k <- pSel <- pGroup <- groups <- NULL 
- blockUpdateType <- blockUpdate[[1]]
- 
- switch(blockUpdateType,
- "correlation" = {
- h <- blockUpdate$h
- k <- blockUpdate$k
- pSel <- blockUpdate$pSel
- pGroup <- blockUpdate$pGroup 
- },
- "random"={
- k <- blockUpdate$k 
- },
- "user"= {
- groups <- blockUpdate$groups
- pSel <- blockUpdate$pSel
- pGroup <- blockUpdate$pGroup 
- })
- 
- return(list(blockUpdateType = blockUpdateType, h = h, k = k, pSel = pSel,
- pGroup = pGroup, groups = groups))
- }
- 
+#' Determine the groups of correlated parameters
+#' @author Stefan Paul
+#' @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){
+
+ settings <- getBlockSettings(blockSettings)
+ blockUpdateType <- settings$blockUpdateType
+
+ switch(blockUpdateType,
+ "correlation" = {
+ ## (Pair wise) Correlation in the parameters
+ cormat <- abs(cor(chain[,1:(ncol(chain)-3),sample(1:dim(chain)[3],1)]))
+ diag(cormat) <- 0
+ # Correct for NA and Inf values as this could cause error in as.dist()
+ cormat[c(which(is.na(cormat)),which(cormat == Inf),which(cormat == -Inf)) ] <- 0
+ tree <- hclust(as.dist(1-cormat)) # get tree based on distance(dissimilarity = 1-cor).
+ cT <- cutree(tree, k = settings$k, h = settings$h) # get groups. With h we can manipulate the strength of the interaction.
+ },
+ "user" = { 
+ cT <- settings$groups
+ },
+ "random" = {
+ pool <- c(1:settings$k, sample(1:settings$k, (ncol(chain)-3-settings$k)))
+ cT <- sample(pool)
+ }
+ )
+
+ pSel <- settings$pSel
+ if(is.null(pSel) && is.null(settings$pGroup)) pSel = rep(1,ncol(chain)-3)
+ return(list(cT = cT, pGroup = settings$pGroup, pSel = pSel))
+}
+
+
+#' Determine the parameters in the block update
+#' @param blockSettings settings for block update
+#' @return vector containing the parameter to be updated
+#' @keywords internal
+getBlock <- function(blockSettings){
+ groups <- blockSettings$cT
+ pGroup <- blockSettings$pGroup
+ pSel <- blockSettings$pSel
+
+
+ nGroups = max(groups)
+ if(nGroups == 1) return(1:length(groups))
+ if (is.null(pGroup)) pGroup = rep(1,nGroups)
+ if(length(pSel) > nGroups) pSel <- pSel[1:nGroups]
+ pSel = c(pSel, rep(0,nGroups - length(pSel)))
+ groupsToSample = sample.int(nGroups, 1, prob = pSel)
+
+ selectedGroups = sample.int(nGroups,groupsToSample, prob = pGroup[1:nGroups])
+ GroupMember <- which(is.element(groups,selectedGroups))
+ return(GroupMember)
+
+}
+
+
+#' getblockSettings
+#' @description Transforms the original settings to settings used in the model runs
+#' @param blockUpdate input settings
+#' @return list with block settings
+#' @keywords internal
+getBlockSettings <- function(blockUpdate){
+
+ h <- k <- pSel <- pGroup <- groups <- NULL 
+ blockUpdateType <- blockUpdate[[1]]
+
+ switch(blockUpdateType,
+ "correlation" = {
+ h <- blockUpdate$h
+ k <- blockUpdate$k
+ pSel <- blockUpdate$pSel
+ pGroup <- blockUpdate$pGroup 
+ },
+ "random"={
+ k <- blockUpdate$k 
+ },
+ "user"= {
+ groups <- blockUpdate$groups
+ pSel <- blockUpdate$pSel
+ pGroup <- blockUpdate$pGroup 
+ })
+
+ return(list(blockUpdateType = blockUpdateType, h = h, k = k, pSel = pSel,
+ pGroup = pGroup, groups = groups))
+ }
+

BayesianTools/R/codaFunctions.R

@@ -1,59 +1,59 @@
-#' Function to combine chains
-#' 
-#' @param x a list of MCMC chains
-#' @param merge should chains be merged? (T or F)
-#' @return combined chains
-#' 
-#' @note to combine several chains to a single McmcSamplerList, see \code{\link{createMcmcSamplerList}}
-#' 
-#' @keywords internal
-combineChains <- function(x, merge = T){
- 
- if(merge == T){
- temp1 = as.matrix(x[[1]])
- 
- names = colnames(temp1)
- 
- sel = seq(1, by = length(x), len = nrow(temp1) )
- 
- out = matrix(NA, nrow = length(x) * nrow(temp1), ncol = ncol(temp1))
- out[sel, ] = temp1
- if (length(x) > 1){
- for (i in 2:length(x)){
- out[sel+i-1, ] = as.matrix(x[[i]])
- } 
- } 
- 
- colnames(out) = names
- 
- } else{
- 
- out = as.matrix(x[[1]])
- if (length(x) > 1){
- for (i in 2:length(x)){
- out = rbind(out, as.matrix(x[[i]]))
- } 
- } 
- }
- 
- return(out)
-}
-
-
-
-#' Helper function to change an object to a coda mcmc class,
-#'
-#' @param chain mcmc Chain
-#' @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 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){
- attr(chain, "mcpar") <- c(start, end, thin)
- attr(chain, "class") <- "mcmc"
- chain
-}
-
-
+#' Function to combine chains
+#' 
+#' @param x a list of MCMC chains
+#' @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}}
+#' 
+#' @keywords internal
+combineChains <- function(x, merge = T){
+
+ if(merge == T){
+ temp1 = as.matrix(x[[1]])
+
+ names = colnames(temp1)
+
+ sel = seq(1, by = length(x), len = nrow(temp1) )
+
+ out = matrix(NA, nrow = length(x) * nrow(temp1), ncol = ncol(temp1))
+ out[sel, ] = temp1
+ if (length(x) > 1){
+ for (i in 2:length(x)){
+ out[sel+i-1, ] = as.matrix(x[[i]])
+ } 
+ } 
+
+ colnames(out) = names
+
+ } else{
+
+ out = as.matrix(x[[1]])
+ if (length(x) > 1){
+ for (i in 2:length(x)){
+ out = rbind(out, as.matrix(x[[i]]))
+ } 
+ } 
+ }
+
+ return(out)
+}
+
+
+
+#' Helper function to change an object to a coda mcmc class,
+#'
+#' @param chain mcmc Chain
+#' @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 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){
+ attr(chain, "mcpar") <- c(start, end, thin)
+ attr(chain, "class") <- "mcmc"
+ chain
+}
+
+