R/analysis_exp_4.R

In rchan26/exp4FusionRCPP: Hierarchical and Sequential Monte Carlo Fusion for exp(-(x^4)/2)

#' 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
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/4, nsamples = 100000,  proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4)
#' test <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 2), time_schedule = rep(1, 2),
#'                                 start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
#'                                 base_samples = input_samples, study = T)
#'
#' # plot result
#' plot_levels_hier_exp_4(test, from = -3, to = 3, plot_rows = 3, plot_columns = 1)
#'
#' @export
plot_levels_hier_exp_4 <- function(hier_result, 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(target_density_exp_4(x), 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(tempered_target_density_exp_4(x, 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
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/4, nsamples = 100000,  proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4)
#' test <- parallel_sequential_fusion_TA_exp_4(N_schedule = rep(10000, 3), global_T = 1,
#'                                             start_beta = 1/4, base_samples = input_samples, study = T)
#'
#' # plot results
#' plot_levels_seq_exp_4(test, from = -3, to = 3, plot_rows = 2, plot_columns = 2)
#'
#' @export
plot_levels_seq_exp_4 <- function(seq_result, 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(target_density_exp_4(x), 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(tempered_target_density_exp_4(x, 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))
}


#' Difference density to the true target density exp(-(x^4)/2)
#'
#' This function obtains the 'difference density', i.e. the density |pi - tilde{pi}|, where tilde{pi} is the kernel density of the samples
#'
#' @param samples samples of exp(-(x^4)/2)
#' @param bw bandwidth for kernel density estimation (same bw as for density())
#'
#' @return difference density of kernel density of the samples and the true density, exp(-(x^4)/2)
#'
#' @examples
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/4, nsamples = 100000,  proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4)
#' test <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 2), time_schedule = rep(1, 2),
#'                                 start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
#'                                 base_samples = input_samples, study = T)
#'
#' # plot true density
#' curve(target_density_exp_4(x), -3, 3)
#' # plot density of samples
#' lines(density(test$samples[[1]]))
#' # plot density of the difference between samples obtained and true density
#' lines(diff_density_true_exp_4(test$samples[[1]]))
#'
#'  @export
diff_density_true_exp_4 <- function(samples, bw = "nrd0") {
  # calculate kernel density estimate for samples
  kde <- density(samples, bw = bw)

  # obtain x and y values from kde
  density_x <- kde$x
  density_y <- kde$y

  # for each value of density_x, we obtain the true value of the pdf at this point
  true_y <- target_density_exp_4(x = density_x)

  # calculate differences between true_y and density_y
  diff <- abs(true_y - density_y)
  return(list('x' = density_x, 'y' = diff))
}

#' Total variation between kernel density of samples to the true target density exp(-(x^4)/2)
#'
#' This function obtains the total variation, i.e. the integral of density |pi - tilde{pi}| / 2,
#' where tilde{pi} is the kernel density of the samples (note we divide by two so that this is bounded by 1)
#'
#' @param samples samples of exp(-(x^4)/2)
#' @param bw bandwidth for kernel density estimation (same bw as for density())
#'
#' @return the total variation between samples and the true density, exp(-(x^4)/2)
#'
#' @examples
#' input_samples <- base_rejection_sampler_exp_4(beta = 1/4, nsamples = 100000,  proposal_mean = 0, proposal_sd = 1.5, dominating_M = 1.4)
#' test <- parallel_h_fusion_exp_4(N_schedule = rep(10000, 2), time_schedule = rep(1, 2),
#'                                 start_beta = 1/4, C_schedule = rep(2, 2), L = 3,
#'                                 base_samples = input_samples, study = T)
#'
#' # plot true density
#' curve(target_density_exp_4(x), -3, 3)
#' # plot density of samples
#' lines(density(test$samples[[1]]))
#' # plot density of the difference between samples obtained and true density
#' lines(diff_density_true_exp_4(test$samples[[1]]))
#'
#' # print total variation
#' print(total_variation_true_exp_4(test$samples))
#'
#'  @export
total_variation_true_exp_4 <- function(samples, bw = "nrd0") {
  # calculate the difference between the kde and true density
  diff <- diff_density_true_exp_4(samples, bw)
  # calculating the total variation
  tv <- sfsmisc::integrate.xy(x = diff$x, fx = diff$y) / 2
  return(tv)
}


# diff_density_sampled_exp_4 <- function(samples, sampled_exp_4, bw = "nrd0") {
#   # if samples for exp_4 are not given, we use rejection sampling to obtain samples
#   if (missing(sampled_exp_4)) {
#     sampled_exp_4 <- sample_from_fc(N = length(samples), proposal_mean = 0, proposal_sd = 1, dominating_M = 1.35, beta = 1)
#   }
#
#   # samples that we wish to compare must have the same length - gives a warning
#   if (length(samples) != sampled_exp_4) {
#     warning('samples and sampled_exp_4 are not of the same length')
#   }
#
#   # calculate the kernel density estimates for the samples that we wish to compare
#   kde_samples <- density(samples, bw = bw)
#   kde_samples_exp_4 <- density(sampled_exp_4, bw = bw)
# }