Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
### 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")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment