Constants in the code below are named as follows:
Draw a random sample from an exponential distribution with mean 900 to represent the lifetime of 60 watt light bulbs.
lambda <- 1/900
pop <- rexp(10000, rate = lambda)
The distribution of the lifetimes in this sample is as follows:
plot(density(pop), main = "Distribution of lifetimes of 60 Watt lightbulbs",
xlab = "")
The t-interval method assumes that the underlying distribution is approximately normal. If we were to generate 800 samples of sample size 30 from this exponential distribution and the calculate 95% t-intervals, what is the proportion of intervals that include the true mean 900? Is this what you would expect? Explain.
sims = 800
n = 15
tsamples <- replicate(sims, rexp(n, rate = lambda))
mu = 1/lambda
results <- as.numeric(sims)
t.int <- matrix(FALSE, sims, 2)
for (i in 1:sims) {
t.int[i, ] <- t.test(tsamples[, i], conf.level = 0.95)$conf.int
results[i] <- t.int[i, 1] < mu & t.int[i, 2] > mu
}
sum(results)/sims
## [1] 0.9125
Since the distribution is exponential, but assumed to be normal, I would expect the percentage of intervals containing the true mean to be fewer than 95% of the total intervals. As you can see, only 0.9125 of the intervals contain the true mean.
A 95% bootstrap percentile confidence interval is calculated by resampling one original sample many times and then using the 2.5 and 75.5 percentile of the resulting means.
This method does not have any assumptions about the underlying distribution. If we were to calculate 800 95% confidence interval based on the bootstrap percentile method, what is the proportion of intervals that include the true mean 900?
sims <- 800
n <- 15 # small sample
boot.ci.pctle <- matrix(FALSE, sims, 2)
results <- as.numeric(sims)
for (i in 1:sims) {
pop <- rexp(10000, rate = lambda) # this is the population
x <- sample(pop, n) # this is the original sample that will be resampled
mu <- mean(pop) # this is the true mean
B <- 1000 # this is the number of resamples
bootsamples <- matrix(sample(x, size = n * B, replace = T), nrow = n, ncol = B)
bootEst <- apply(bootsamples, 2, mean)
boot.ci.pctle[i, ] <- quantile(bootEst, probs = c(0.025, 0.975))
results[i] <- boot.ci.pctle[i, 1] < mu & boot.ci.pctle[i, 2] > mu
}
sum(results)/sims
## [1] 0.9012
The “studentized” bootstrap confidence intervals are calculated by resampling one original sample many times, and by adjusting the means of the bootstrap estimates before choosing the corresponding percentiles. This method does not have any assumptions about the underlying distribution.
If we were to calculate 800 95% confidence interval based on the studentized bootstrap method, what is the proportion of intervals that include the true mean 900?
sims <- 800
n <- 15 # small sample
boot.ci.student <- matrix(FALSE, sims, 2)
results <- as.numeric(sims)
for (i in 1:sims) {
pop <- rexp(10000, rate = lambda) # this is the population
x <- sample(pop, n) # this is the original sample that will be resampled
mu <- mean(pop) # this is the true mean
B <- 1000 # this is the number of resamples
bootsamples <- matrix(sample(x, size = n * B, replace = T), nrow = n, ncol = B)
bootEst <- (apply(bootsamples, 2, mean) - mean(x))/(apply(bootsamples, 2,
sd)/sqrt(n))
boot.ci.student[i, ] <- mean(x) - quantile(bootEst, probs = c(0.975, 0.025)) *
sd(x)/sqrt(n)
results[i] <- boot.ci.student[i, 1] < mu & boot.ci.student[i, 2] > mu
}
sum(results)/sims
## [1] 0.9475
As we can see from the data, only 0.9475 of the intervals contain the true mean.
Compare the coverage of the three methods: t-interval, bootstrap percentile method, and studentized bootstrap method. This was calculated from a skewed distribution (exponential) with a small sample size (n=15). What proportion of coverage would you expect for each of these three methods if the underlying distribution was normal? For example, generate the data from rnorm(15, mean=900, sd=50)?
Comparing the three methods, I found that the studentized bootstrap method was closest to 95%, the t-interval was second best and the bootstrap percentile method was the worst - lowest percentage. However, the t-interval and bootstrap percentile method were very close. If there was a normal distribution underlying the data, I would expect the t-interval method to have about 95% of the intervals covering the true mean, but I would expect both bootstrap methods to have roughly the same coverage as they do now since neither of those methods make any assumptions about the underlysing distribution in the initial case.
For example: for t-intervals method:
sims <- 800
n <- 15
tsamples <- replicate(sims, rnorm(10000, mean = 900, sd = 50))
mu <- 900
results <- as.numeric(sims)
t.int <- matrix(FALSE, sims, 2)
for (i in 1:sims) {
t.int[i, ] <- t.test(tsamples[, i], conf.level = 0.95)$conf.int
results[i] <- t.int[i, 1] < mu & t.int[i, 2] > mu
}
sum(results)/sims
## [1] 0.96
So 0.96 intervals are covered.
For the bootstrap percentile method:
sims <- 800
n <- 15
boot.ci.pctle <- matrix(FALSE, sims, 2)
results <- as.numeric(sims)
for (i in 1:sims) {
pop <- rnorm(10000, mean = 900, sd = 50) # this is the population
x <- sample(pop, n) # this is the original sample that will be resampled
mu <- mean(pop) # this is the true mean
B <- 1000 # this is the number of resamples
bootsamples <- matrix(sample(x, size = n * B, replace = T), nrow = n, ncol = B)
bootEst <- apply(bootsamples, 2, mean)
boot.ci.pctle[i, ] <- quantile(bootEst, probs = c(0.025, 0.975))
results[i] <- boot.ci.pctle[i, 1] < mu & boot.ci.pctle[i, 2] > mu
}
sum(results)/sims
## [1] 0.8912
Thus, 0.8912 of the intervals cover the true mean.
For the studentized bootstrap method:
sims <- 800
n <- 15 # small sample
boot.ci.student <- matrix(FALSE, sims, 2)
results <- as.numeric(sims)
for (i in 1:sims) {
pop <- rnorm(10000, mean = 900, sd = 50) # this is the population
x <- sample(pop, n) # this is the original sample that will be resampled
mu <- mean(pop) # this is the true mean
B <- 1000 # this is the number of resamples
bootsamples <- matrix(sample(x, size = n * B, replace = T), nrow = n, ncol = B)
bootEst <- (apply(bootsamples, 2, mean) - mean(x))/(apply(bootsamples, 2,
sd)/sqrt(n))
boot.ci.student[i, ] <- mean(x) - quantile(bootEst, probs = c(0.975, 0.025)) *
sd(x)/sqrt(n)
results[i] <- boot.ci.student[i, 1] < mu & boot.ci.student[i, 2] > mu
}
sum(results)/sims
## [1] 0.9563
So 0.9562 of the intervals contain the true mean.