diff --git a/DESCRIPTION b/DESCRIPTION index 774f6fb..39b27d5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,7 @@ Description: drjacoby is an R package for performing Bayesian inference via License: MIT + file LICENSE Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 LinkingTo: Rcpp Imports: diff --git a/NAMESPACE b/NAMESPACE index 43db110..50753fa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ export(plot_scatter) export(plot_trace) export(run_mcmc) export(sample_chains) +export(tune_rungs) import(dplyr) import(ggplot2) importFrom(GGally,ggpairs) @@ -25,6 +26,7 @@ importFrom(coda,geweke.diag) importFrom(coda,mcmc) importFrom(grDevices,grey) importFrom(magrittr,"%>%") +importFrom(stats,approx) importFrom(stats,cor) importFrom(stats,pnorm) importFrom(stats,quantile) diff --git a/R/main.R b/R/main.R index d548e98..8539bb9 100644 --- a/R/main.R +++ b/R/main.R @@ -609,6 +609,58 @@ deploy_chain <- function(args) { return(ret) } +#------------------------------------------------ +#' @title Tune mcmc rungs +#' +#' @description Takes a completed mcmc object that has been run with parallel +#' tempering, and uses the acceptance rates to recalculate the optimal number +#' and distribution of rungs. +#' +#' @param x an object of class \code{drjacoby_output}, with the MCMC completed +#' and run with parallel tempering over multiple rungs. +#' @param target_acceptance the target acceptance rate between rungs. Higher +#' values will lead to more rungs being needed. +#' +#' @importFrom stats approx +#' @export + +tune_rungs <- function(x, target_acceptance = 0.23) { + + # avoid "no visible binding" note + link <- value <- NULL + + # check inputs + assert_class(x, "drjacoby_output") + assert_neq(length(x$diagnostics$mc_accept), 1, message = "run_mcmc() must have been run with multiple temperature rungs") + assert_single_bounded(target_acceptance) + + # get acceptance rates and check all non-zero + acceptance <- x$diagnostics$mc_accept |> + group_by(link) |> + summarise(value = mean(value)) |> + pull(value) + + assert_gr(acceptance, 0, message = "acceptance rates must be positive for all pairs of rungs") + n_rungs <- length(acceptance) + 1 + + # calculate cumulative barrier Lambda + Lambda <- sum(1 - acceptance) + + # calculate new number of rungs + n_rungs_new <- ceiling(Lambda / (1 - target_acceptance)) + 1 + + # calculate new distribution of rungs by linear interpolation + beta_old <- x$parameters$beta_manual + y_old <- c(0, cumsum(1 - acceptance)) + #plot(beta_old, y_old) + y_new <- seq(0, Lambda, l = n_rungs_new) + result <- approx(y_old, beta_old, xout = y_new) + #points(result$y, result$x, col = 2) + beta_new <- result$y + + return(beta_new) +} + #------------------------------------------------ # update progress bar # pb_list = list of progress bar objects @@ -625,5 +677,7 @@ update_progress <- function(pb_list, name, i, max_i, close = TRUE) { } } +#------------------------------------------------ # Deal with user input cpp not being defined globalVariables(c("create_xptr")) + diff --git a/R_ignore/deploy.R b/R_ignore/deploy.R index b8d0a05..02ea2b0 100644 --- a/R_ignore/deploy.R +++ b/R_ignore/deploy.R @@ -42,7 +42,17 @@ mcmc <- run_mcmc(data = list(x = x), loglike = r_loglike, logprior = r_logprior, burnin = 1e3, - samples = 1e3) + samples = 1e3, + chains = 2, + beta_manual = beta_tuned) + +plot_mc_acceptance(mcmc) + +beta_tuned + +beta_tuned + +beta_tuned <- interleave_with_halfway(beta_tuned) plot_trace(mcmc, show = "mu", phase = "burnin") diff --git a/man/tune_rungs.Rd b/man/tune_rungs.Rd new file mode 100644 index 0000000..1f1e6f1 --- /dev/null +++ b/man/tune_rungs.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/main.R +\name{tune_rungs} +\alias{tune_rungs} +\title{Tune mcmc rungs} +\usage{ +tune_rungs(x, target_acceptance = 0.23) +} +\arguments{ +\item{x}{an object of class \code{drjacoby_output}, with the MCMC completed +and run with parallel tempering over multiple rungs.} + +\item{target_acceptance}{the target acceptance rate between rungs. Higher +values will lead to more rungs being needed.} +} +\description{ +Takes a completed mcmc object that has been run with parallel + tempering, and uses the acceptance rates to recalculate the optimal number + and distribution of rungs. +}