bootstrapFunctions.R

### BOOTSTRAP FUNCTIONS ###

# Install required packages if necessary.
#install.packages("MCMCpack", "rootSolve")

# Load required packages.
library(ggplot2)

# General function for computing boostrap samples for a numerical statistic. 
# Input: input data, number of bootstrap samples B, bootstrap sample size n,
#        function for computing statistic func.
# Output: numeric vector of bootstrap estimates.
getBootstrapEstimates = function(data, B, n, func) {
  samples = numeric(B)
  for (i in 1:B) {
    bootstrapSample = sample(data, n, replace = T)
    samples[i] = func(bootstrapSample)
  }
  
  return(samples)
}

# General function for extracting confidence intervals from bootstrap estimates.
# Input: a bootrapped sample, interval center point, desired confidence.
# Output: confidence interval centered around center.
getConfidenceInterval = function(bootstrapSample, center, confidence) {
  B = length(bootstrapSample)
  
  # Sort bootstrap sample values with center point.
  sortedEstimatesWithCenter = sort(c(bootstrapSample, center))
  
  # Get a quantile value for the center point relative to the bootstrap sample.
  centerQuantile = which(sortedEstimatesWithCenter == center)[1] / (B + 1)
  
  # Compute lower and upper bounds of confidence interval about the center point.
  upperBound = min(centerQuantile + (confidence / 2), 1)
  lowerBound = max(centerQuantile - (confidence / 2), 0)
  
  # Create confidence interval about the center point.
  confidenceInterval = quantile(
    sortedEstimatesWithCenter,
    c(lowerBound, upperBound))
  
  return(confidenceInterval)
}

### BOOTSTRAP PLOTTING FUNCTION ###

getBootstrapPlots = function(
  bootstrapSample, 
  center, 
  confidence,
  title,
  upper = 1200) {
  # Obtain mean and confidence interval.
  confidenceInterval = getConfidenceInterval(bootstrapSample, center,confidence)
  
  # Create data frames for plotting purposes.
  bootstrapSampleDF = data.frame(bootstrapSample)
  names(bootstrapSampleDF)[1] = "x"
  confidenceDF = data.frame(
    c(
      rep(confidenceInterval[1], 2), 
      rep(confidenceInterval[2], 2), 
      confidenceInterval[1]),
    c(0, upper, upper, 0, 0))
  names(confidenceDF) = c("x", "y")
  
  # Create plot of confidence interval centered around empirical quantile estimate
  # bootstrap sample mean.
  
  bootstrapPlot <- ggplot() + 
    labs(
      title = title, 
      x ="90% Quantile Estimate", 
      y = NULL) +
    theme(plot.title = element_text(face = "bold", size = 25)) +
    theme_bw() + 
    geom_histogram(
      aes(x = x),
      data = bootstrapSampleDF,
      fill = "orange", 
      alpha = 1, 
      colour = "grey",
      lwd = 0.2,
      binwidth = 3) + 
    geom_polygon(
      aes(x = x, y = y),
      data = confidenceDF,
      fill = "red",
      lwd = 0,
      alpha = 0.5) + 
    geom_vline(xintercept = center, colour = "blue")
}