Nothing
#' @title Calculate marginal likelihood by thermodynamic integration (LTI)
#' @description Calculate the marginal likelihood from a logfile generated by \code{\link{mcmc_bite}} with thermodynamic integration (Lartillot and Philippe, 2006) or stepping stone (Xie et al., 2011).
#' @param mcmc.log the output file of a \code{\link{mcmc_bite}} run
#' @param burnin number or proportion of iteration to delete
#' @param method one of "TI" for thermodynamic integration and "SS" for stepping stone integration (the default)
#' @export
#' @return a length one numeric double giving the marginal likelihood of the model.
#' @author Theo Gaboriau and Simon Joly
#' @examples
#' ## Load test data
#' data(Anolis_traits)
#' data(Anolis_tree)
#' data(Anolis_map)
#'
#' ## Run a MCMC chain with thermodynamic Integration
#' set.seed(300)
#' my.jive <- make_jive(Anolis_tree, Anolis_traits[,-3],
#' model.priors = list(mean="BM", logvar="OU"))
#' bite_ex <- tempdir()
#' logfile <- sprintf("%s/my.jive_mcmc_TI.log", bite_ex)
#' mcmc_bite(my.jive, log.file=logfile, ncat=10, sampling.freq=10,
#' print.freq=100, ngen=1000, burnin=0)
#'
#' ## import the results in R
#' res <- read.csv(logfile, header = TRUE, sep = "\t")
#'
#' mlikTI <- marginal_lik(res, burnin = 0.1, method = "TI")
#' mlikTI
#'
#' mlikSS <- marginal_lik(res, burnin = 0.1, method = "SS")
#' mlikSS
#' @encoding UTF-8
marginal_lik <- function(mcmc.log, burnin = 0, method = "SS") {
if (!("temperature" %in% colnames(mcmc.log))) stop('No \'temperature\' column in the log file')
temp <- unique(mcmc.log$temperature)
if(burnin>0){
if(burnin < 1) burnin <- sapply(temp, function(t) quantile(mcmc.log$iter[mcmc.log$temperature == t], burnin))
else burnin <- sapply(temp, function(t) min(mcmc.log$iter[mcmc.log$temperature == t])) + burnin
burn <- unlist(lapply(1:length(temp), function(t) mcmc.log[mcmc.log$temperature == temp[t],"iter"] <= burnin[t]))
} else {
burn <- rep(FALSE, nrow(mcmc.log))
}
# The marginal likelihood
lik.ti = 0
# Calculate the intervals in temperature
intervals <- abs(diff(unique(mcmc.log$temperature)))
if(method == "TI"){
# Get mean marginal likelihood for each temperature
lik <- sapply(temp, function(t) mean(rowSums(mcmc.log[mcmc.log$temperature == t & !burn, grepl("prior", colnames(mcmc.log))], na.rm = TRUE)))
if (length(lik)<=1) warning('Only one temerature for the thermodynamic integration')
# Estimate the integral
for (i in 1:length(intervals)) {
lik.ti = lik.ti + ((lik[i]+lik[i+1])/2 * intervals[i])
}
} else if(method == "SS"){
## Estimate RSS for each k and sum it
for(i in 1:length(intervals)) {
betaliks <- rowSums(mcmc.log[mcmc.log$temperature == temp[i] & !burn, grepl("prior", colnames(mcmc.log))])
Lmax <- max(betaliks)
n <- length(betaliks)
lik.ti = lik.ti + intervals[i]*Lmax + log(1/n * sum(exp(intervals[i]*(betaliks - Lmax))))
}
} else {
stop(sprintf("%s is not a supported method"), method)
}
return(lik.ti)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.