R/analysis_mixG.R

In rchan26/mixGaussTempering: What the Package Does (One Line, Title Case)

#' Plot kernel density of samples in hierarchical tree
#'
#' This function plots the kernel densities of the samples in the hierarchical tree (in red)
#' and the tempered target densities that they were targetting (in blue)
#'
#' @param hier_result hierarchical Monte Carlo fusion result
#' @param from,to the range over which the function will be plotted: integer (if want to have it constant throughout)
#'                or vector of length L, where L is the number of levels in the hierarchy
#' @param plot_rows number of rows in plot (should make is so plot_rows x plot_columns = L)
#' @param plot_columns number of rows in plot (should make is so plot_rows x plot_columns = L)
#' @param whole boolean value: defaults to F, determines whether or not to plot kde for all samples in the level or for each node (T)
#' @param bw the smoothing bandwidth used: integer (if want to have it constant throughout) or vector of length L,
#'           where L is the number of levels in the hierarchy
#'
#' @return none
#'
#' @examples
#' # setting variables
#' w_example <- c(0.35, 0.65)
#' m_example <- c(-3, 8)
#' s_example <- c(1, 1.5)
#' b_example <- 1/4
#'
#' # sampling from tempered density
#' nsamples <- 500000
#' base <- hmc_base_sampler_mixG(w_example, m_example, s_example, b_example, nsamples, 4)
#'
#' test <- parallel_h_fusion_TA_mixG(N_schedule = rep(10000, 2), time_schedule = 1,
#'                                   start_beta = b_example, C_schedule = rep(2, 2), L = 3, base_samples = base,
#'                                   weights = w_example, means = m_example, sds = s_example, study = T)
#'
#' # plot results with constant from and to
#' plot_levels_hier_mixG(test, weights = w_example, means = m_example, sds = s_example,
#'                       from = -15, to = 20, plot_rows = 3, plot_columns = 1)
#'
#' # plot results with specified bandwidths
#' plot_levels_hier_mixG(test, weights = w_example, means = m_example, sds = s_example,
#'                       from = -15, to = 20, plot_rows = 3, plot_columns = 1, bw = c(0.2, 0.3, 0.4))
#'
#' @export
plot_levels_hier_mixG <- function(hier_result, weights, means, sds, from, to, plot_rows, plot_columns, whole = F, bw = rep('nrd0', length(hier_result$samples))) {
  if (length(from)==1) {
    from <- rep(from, length(hier_result$samples))
  }
  if (length(to)==1) {
    to <- rep(to, length(hier_result$samples))
  }
  if (length(bw)==1) {
    bw <- rep(bw, length(hier_result$samples))
  }
  par(mfrow=c(plot_rows, plot_columns))
  base_beta <- hier_result$output_beta[length(hier_result$output_beta)]
  curve(dnorm_mix(x, weights, means, sds), from = from[1], to = to[1], col = 'blue',
        main = paste('beta =', 1/base_beta, '/', 1/base_beta), lwd = 2, ylab = 'pdf')
  lines(density(hier_result$samples[[1]], bw = bw[1]), col = 'red', lwd = 1)
  for (k in 2:length(hier_result$samples)) {
    curve(dnorm_mix_tempered(x, weights, means, sds, beta = hier_result$output_beta[k]), from = from[k], to = to[k],
          col = 'blue', main = paste('beta =', hier_result$output_beta[k]/base_beta, '/', 1/base_beta), lwd = 2, ylab = 'pdf')
    if (whole == T) {
      lines(density(unlist(hier_result$samples[[k]]), bw = bw[k]), col = 'red', lwd = 1)
    } else {
      for (j in 1:length(hier_result$samples[[k]])) {
        lines(density(hier_result$samples[[k]][[j]], bw = bw[k]), col = 'red', lwd = 1)
      }
    }
  }
  par(mfrow=c(1,1))
}

#' Plot kernel density of samples in sequential tree
#'
#' This function plots the kernel densities of the samples in the sequential tree (in red)
#' and the tempered target densities that they were targetting (in blue)
#'
#' @param hier_result sequential Monte Carlo fusion result
#' @param from,to the range over which the function will be plotted: integer (if want to have it constant throughout)
#'                or vector of length L, where L is the number of levels in the hierarchy
#' @param plot_rows number of rows in plot (should make is so plot_rows x plot_columns = L)
#' @param plot_columns number of rows in plot (should make is so plot_rows x plot_columns = L)
#' @param whole boolean value: defaults to F, determines whether or not to plot kde for all samples in the level or for each node (T)
#' @param bw the smoothing bandwidth used: integer (if want to have it constant throughout) or vector of length L,
#'           where L is the number of levels in the hierarchy
#'
#' @return none
#'
#' @examples
#' # setting variables
#' w_example <- c(0.35, 0.65)
#' m_example <- c(-3, 8)
#' s_example <- c(1, 1.5)
#' b_example <- 1/4
#'
#' # sampling from tempered density
#' nsamples <- 500000
#' base <- hmc_base_sampler_mixG(w_example, m_example, s_example, b_example, nsamples, 4)
#'
#' test <- parallel_sequential_fusion_TA_mixG(N_schedule = rep(10000, 3), global_T = 1,
#'                                            start_beta = b_example, base_samples = base,
#'                                            weights = w_example, means = m_example, sds = s_example, study = T)
#'
#' # plot results
#' plot_levels_seq_mixG(test, weights = w_example, means = m_example, sds = s_example,
#'                      from = -15, to = 20, plot_rows = 2, plot_columns = 2, bw = c(0.2, 0.3, 0.3, 0.35))
#'
#' @export
plot_levels_seq_mixG <- function(seq_result, weights, means, sds, from, to, plot_rows, plot_columns, bw = rep('nrd0', length(seq_result$samples))) {
  if (length(from)==1) {
    from <- rep(from, length(seq_result$samples))
  }
  if (length(to)==1) {
    to <- rep(to, length(seq_result$samples))
  }
  if (length(bw)==1) {
    bw <- rep(bw, length(seq_result$samples))
  }
  par(mfrow=c(plot_rows, plot_columns))
  base_beta <- seq_result$output_beta[length(seq_result$output_beta)]
  curve(dnorm_mix(x, weights, means, sds), from = from[1], to = to[1], col = 'blue',
        main = paste('beta =', 1/base_beta, '/', 1/base_beta), lwd = 2, ylab = 'pdf')
  lines(density(seq_result$samples[[1]], bw = bw[1]), col = 'red', lwd = 1)
  for (k in 2:length(seq_result$samples)) {
    curve(dnorm_mix_tempered(x, weights, means, sds, beta = seq_result$output_beta[k]), from = from[k], to = to[k],
          col = 'blue', main = paste('beta =', seq_result$output_beta[k]/base_beta, '/', 1/base_beta), lwd = 2, ylab = 'pdf')
    if (length(seq_result$samples[[k]])==1) {
      lines(density(seq_result$samples[[k]], bw = bw[k]), col = 'red', lwd = 1)
    } else {
      lines(density(unlist(seq_result$samples[[k]]), bw = bw[k]), col = 'red', lwd = 1)
    }
  }
  par(mfrow=c(1,1))
}