$$\newcommand{\E}{\mathbb{E}}\newcommand{\Var}{\operatorname{Var}}\newcommand{\Cov}{\operatorname{Cov}}\newcommand{\iid}{\stackrel{\text{iid}}{\sim}}\newcommand{\ind}{\perp\!\!\!\perp}\newcommand{\Normal}{\mathcal{N}}\newcommand{\Uniform}{\mathcal{U}}\newcommand{\Exp}{\operatorname{Exp}}\newcommand{\Pois}{\operatorname{Pois}}\newcommand{\Ber}{\operatorname{Ber}}\newcommand{\Bin}{\operatorname{Bin}}\newcommand{\Gam}{\operatorname{Gamma}}\newcommand{\Betad}{\operatorname{Beta}}\newcommand{\argmax}{\mathop{\operatorname{arg\,max}}}\newcommand{\argmin}{\mathop{\operatorname{arg\,min}}}\newcommand{\bx}{\mathbf{x}}\newcommand{\by}{\mathbf{y}}\newcommand{\bz}{\mathbf{z}}\newcommand{\bu}{\mathbf{u}}\newcommand{\bX}{\mathbf{X}}\newcommand{\bL}{\mathbf{L}}\newcommand{\bI}{\mathbf{I}}\newcommand{\bmu}{\boldsymbol{\mu}}\newcommand{\bSigma}{\boldsymbol{\Sigma}}\newcommand{\btheta}{\boldsymbol{\theta}}\newcommand{\R}{\mathbb{R}}\newcommand{\one}{\mathbf{1}}$$

Lab 5: Bootstrap Confidence Intervals

STAT 536 — Statistical Computing

Author

Soumojit Das

Published

Fall 2026

Learning Objectives

By the end of today’s lab, you will:

  • Implement the percentile, basic, bootstrap-t, and BCa intervals from the algorithm boxes below
  • Estimate the coverage of an interval method by simulation, on both sides of the interval
  • Compare the bootstrap intervals with the \(t\) interval and with an exact interval on one problem
  • Recognize a statistic the nonparametric bootstrap fails for, and repair it with the parametric bootstrap

Setup

Figure labels are typeset by LaTeX. tikz_png() draws a ggplot through tikzDevice, compiles it with pdflatex, and returns a PNG. It needs a LaTeX installation (TinyTeX is enough) and installs tikzDevice, pdftools, and png on first use.

library(ggplot2)
library(patchwork)
set.seed(536)
theme_set(theme_bw(base_size = 14))

pal <- c(red = "#D62828", blue = "#2E86AB", green = "#06A77D", orange = "#F77F00", dark = "#023047", grey = "#6a6a6a")

find_pdflatex <- function() {
  cand <- c(Sys.which("pdflatex"),
            path.expand("~/Library/TinyTeX/bin/universal-darwin/pdflatex"),
            file.path(Sys.getenv("APPDATA"), "TinyTeX/bin/windows/pdflatex.exe"),
            path.expand("~/.TinyTeX/bin/x86_64-linux/pdflatex"))
  cand <- cand[nzchar(cand) & file.exists(cand)]
  if (!length(cand)) stop("pdflatex not found; install TinyTeX with tinytex::install_tinytex()")
  cand[[1]]
}

tikz_png <- function(p, png_file, width = 9, height = 5, dpi = 150) {
  for (pkg in c("tikzDevice", "pdftools", "png"))
    if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
  options(tikzLatex = find_pdflatex(),
          tikzDocumentDeclaration = r"(\documentclass[10pt]{article})",
          tikzLatexPackages = c(r"(\usepackage{tikz})",
                                r"(\usepackage[active,tightpage,psfixbb]{preview})",
                                r"(\PreviewEnvironment{pgfpicture})",
                                r"(\setlength\PreviewBorder{0pt})",
                                r"(\usepackage[T1]{fontenc})", r"(\usepackage{lmodern})",
                                r"(\usepackage{amsmath})"))
  work <- file.path(tempdir(), "tikz_png"); dir.create(work, showWarnings = FALSE)
  tex <- file.path(work, "fig.tex"); pdf <- file.path(work, "fig.pdf")
  tikzDevice::tikz(tex, width = width, height = height, standAlone = TRUE)
  print(p); dev.off()
  ok <- system2(getOption("tikzLatex"),
                c("-interaction=batchmode", "-halt-on-error",
                  paste0("-output-directory=", shQuote(work)), shQuote(tex)),
                stdout = FALSE, stderr = FALSE)
  if (ok != 0 || !file.exists(pdf)) stop("pdflatex failed; see ", file.path(work, "fig.log"))
  dir.create(dirname(png_file), recursive = TRUE, showWarnings = FALSE)
  png::writePNG(pdftools::pdf_render_page(pdf, page = 1, dpi = dpi), png_file)
  png_file
}

# tick labels in math mode, so minus signs and decimals are set by LaTeX
tex_num <- function(v) ifelse(is.na(v), NA,
  paste0("$", format(v, scientific = FALSE, drop0trailing = TRUE, trim = TRUE), "$"))
tex_theme <- theme(text = element_text(size = 17), plot.title = element_text(hjust = 0.5),
                   legend.position = "top", legend.key.width = unit(1.2, "cm"))

B <- 999          # resamples; (B + 1) * alpha / 2 is an integer, so the quantile cuts fall on replicates
alpha <- 0.05     # nominal level 1 - alpha

Notation and Algorithms

Symbol Meaning
\(F\); \(\hat F_n\) the unknown distribution of the data; the empirical distribution, mass \(1/n\) on each \(X_i\)
\(\theta = T(F)\); \(\hat\theta = T(\hat F_n)\) the parameter; the same rule \(T\) on the data
\(X_1^*, \ldots, X_n^*\); \(\hat\theta^*\); \(\hat\theta^{*(b)}\) a resample, \(n\) draws with replacement; the rule \(T\) on it; on resample \(b\), \(b = 1, \ldots, B\)
\(\hat\theta^*_{(q)}\) the \(q\)-quantile of the replicates \(\hat\theta^{*(1)}, \ldots, \hat\theta^{*(B)}\)
\(P_F\), \(\E_F\), \(\Var_F\) probability, mean, variance over samples from \(F\)
\(P^*\), \(\E^*\), \(\Var^*\) probability, mean, variance over the resampling (the index draws), data fixed
\(\widehat{\text{SE}}\); \(\widehat{\text{SE}}^{*(b)}\) the standard-error rule on the data; on resample \(b\)

A star means computed from a resample; \((b)\) says which one. The nominal level is \(1 - \alpha\) with \(\alpha = 0.05\) throughout.

NoteAlgorithm 1 (Bootstrap replicates)

Input: data \(x_1, \ldots, x_n\), the rule \(T\), and \(B\).

For \(b = 1, \ldots, B\):

  1. draw \(U_1, \ldots, U_n \iid \Uniform\{1, \ldots, n\}\)
  2. set \(X_i^* = x_{U_i}\)
  3. compute \(\hat\theta^{*(b)} = T(X_1^*, \ldots, X_n^*)\)

Output: \(\hat\theta^{*(1)}, \ldots, \hat\theta^{*(B)}\).

NoteAlgorithm 2 (Percentile interval)

\[\text{CI}_{\text{perc}} = \bigl[\hat\theta^*_{(\alpha/2)},\ \hat\theta^*_{(1-\alpha/2)}\bigr].\]

NoteAlgorithm 3 (Basic interval)

\[\text{CI}_{\text{basic}} = \bigl[2\hat\theta - \hat\theta^*_{(1-\alpha/2)},\ 2\hat\theta - \hat\theta^*_{(\alpha/2)}\bigr].\]

NoteAlgorithm 4 (Bootstrap-t interval)

Input: the replicates, \(\widehat{\text{SE}}\) on the data, and \(\widehat{\text{SE}}^{*(b)}\) on each resample \(b\) (the same rule applied to resample \(b\)).

  1. \(t^{*(b)} = \bigl(\hat\theta^{*(b)} - \hat\theta\bigr) \big/ \widehat{\text{SE}}^{*(b)}\) for \(b = 1, \ldots, B\).
  2. \(t^*_{(q)}\): the \(q\)-quantile of \(t^{*(1)}, \ldots, t^{*(B)}\).
  3. \(\text{CI}_{t} = \bigl[\hat\theta - t^*_{(1-\alpha/2)}\,\widehat{\text{SE}},\ \hat\theta - t^*_{(\alpha/2)}\,\widehat{\text{SE}}\bigr]\).

For the mean: \(\widehat{\text{SE}} = S/\sqrt{n}\) and \(\widehat{\text{SE}}^{*(b)} = S^{*(b)}/\sqrt{n}\), with \(S^{*(b)}\) the standard deviation of resample \(b\).

Make sure:

  • each \(t^{*(b)}\) is divided by its own \(\widehat{\text{SE}}^{*(b)}\), from resample \(b\)
  • each \(t^{*(b)}\) is centered at \(\hat\theta\)
  • the upper quantile \(t^*_{(1-\alpha/2)}\) sets the lower endpoint, and the lower quantile \(t^*_{(\alpha/2)}\) sets the upper endpoint
NoteAlgorithm 5 (BCa interval)
  1. \(\hat z_0 = \Phi^{-1}\bigl(\#\{\hat\theta^{*(b)} < \hat\theta\}/B\bigr)\), with \(\Phi\) the standard normal CDF.
  2. For the mean, \(\hat a = \sum_i (x_i - \bar x)^3 \Big/ 6\Bigl[\sum_i (x_i - \bar x)^2\Bigr]^{3/2}\). For a general \(T\), \(\hat a\) comes from the jackknife (next week).
  3. With \(z_q = \Phi^{-1}(q)\): \(\alpha_1 = \Phi\!\left(\hat z_0 + \dfrac{\hat z_0 + z_{\alpha/2}}{1 - \hat a(\hat z_0 + z_{\alpha/2})}\right)\), \(\alpha_2 = \Phi\!\left(\hat z_0 + \dfrac{\hat z_0 + z_{1-\alpha/2}}{1 - \hat a(\hat z_0 + z_{1-\alpha/2})}\right)\).
  4. \(\text{CI}_{\text{BCa}} = \bigl[\hat\theta^*_{(\alpha_1)},\ \hat\theta^*_{(\alpha_2)}\bigr]\).
NoteAlgorithm 6 (Parametric bootstrap)

Input: a model \(\{F_\theta\}\) and the estimate \(\hat\theta\).

For \(b = 1, \ldots, B\):

  1. draw \(X_1^*, \ldots, X_n^* \iid F_{\hat\theta}\)
  2. compute \(\hat\theta^{*(b)} = T(X_1^*, \ldots, X_n^*)\)

Algorithms 2 to 5 then apply unchanged.

NoteDefinition (Coverage)

An interval method is a rule that turns a sample into an interval \([L, U]\). Its coverage under \(F\) at sample size \(n\) is \[P_F(L \leq \theta \leq U),\] the probability, over samples from \(F\), that the interval contains \(\theta\). The nominal level \(1 - \alpha\) is what the method promises. The actual coverage is this probability. An equal-tailed method also promises \(P_F(\theta < L) = P_F(\theta > U) = \alpha/2\). One sample gives one interval. It contains \(\theta\) or it does not. Coverage is a property of the method. Algorithm 7 gives one way to estimate it.

NoteAlgorithm 7 (Coverage by simulation)

Input: an \(F\) with known \(\theta\), a sample size \(n\), an interval method, and \(R\).

For \(r = 1, \ldots, R\):

  1. draw a dataset of size \(n\) from \(F\)
  2. build \([L_r, U_r]\)
  3. record \(C_r = \mathbf{1}\{L_r \leq \theta \leq U_r\}\), \(A_r = \mathbf{1}\{\theta < L_r\}\), \(D_r = \mathbf{1}\{\theta > U_r\}\)

Output: coverage \(\hat C = \frac{1}{R}\sum_r C_r\) with standard error \(\sqrt{\hat C(1 - \hat C)/R}\), and the miss rates \(\frac{1}{R}\sum_r A_r\) and \(\frac{1}{R}\sum_r D_r\), each nominally \(\alpha/2\).

Part 1: Bootstrap Intervals for a Skewed Mean

One sample of \(n = 17\) from \(\Exp(1)\). The parameter is the mean, \(\theta = \mu = 1\). We treat it as unknown. The rule \(T\) is the sample mean, so \(\hat\theta = \bar x\) and \(\widehat{\text{SE}} = S/\sqrt{n}\).

Two intervals need no bootstrap.

  • The \(t\) interval \(\bar x \pm t_{n-1,\,0.975}\,S/\sqrt{n}\). It assumes a normal sample: its pivot \((\bar x - \mu)/(S/\sqrt{n})\) has the \(t_{n-1}\) distribution only then. Here it is an approximation.
  • The exact interval. It assumes the exponential model, which supplies the pivot \(2n\bar x/\mu \sim \chi^2_{2n}\): \[\Bigl[\,2n\bar x \big/ \chi^2_{2n,\,0.975},\ \ 2n\bar x \big/ \chi^2_{2n,\,0.025}\,\Bigr].\]

We build the four bootstrap intervals from Algorithms 2 to 5 and compare all six by coverage (Definition, Algorithm 7).

Exercise 1.1: Resample

Implement Algorithm 1 for the mean with one matrix. Column \(b\) holds resample \(b\).

  1. Fill in the index draw in boot_matrix().
  2. The replicates are the column means.
  3. \(\widehat{\text{SE}}^{*(b)}\) is the column standard deviation over \(\sqrt{n}\).
set.seed(536)
x <- rexp(17)                  # the sample; its true mean is 1
n <- length(x)
theta_hat <- mean(x)
se_hat <- sd(x) / sqrt(n)

boot_matrix <- function(x, B) {
  n <- length(x)
  idx <- ___________                                # TODO: n * B indices, uniform on 1:n, with replacement
  matrix(x[idx], nrow = n)                          # column b is resample b
}

M <- boot_matrix(x, B)
theta_star <- ___________                           # TODO: the B replicates of the mean
se_star <- ___________                              # TODO: the SE rule on each resample (column sd over sqrt(n))
print(c(theta_hat = theta_hat, se_hat = se_hat, se_boot = sd(theta_star)))

Report.

  1. \(\hat\theta\) and \(\widehat{\text{SE}}\).
  2. The bootstrap standard error, the standard deviation of the replicates.

Exercise 1.2: Four Bootstrap Intervals

Fill in the four functions from Algorithms 2 to 5. Each takes the replicates and returns the two endpoints. The \(t\) interval and the exact interval are given. Then run the figure chunk.

set.seed(536)
ci_perc <- function(theta_star) {
  ___________                                       # TODO: Algorithm 2
}
ci_basic <- function(theta_star, theta_hat) {
  ___________                                       # TODO: Algorithm 3
}
ci_boott <- function(theta_star, se_star, theta_hat, se_hat) {
  t_star <- ___________                             # TODO: Algorithm 4, step 1
  q <- ___________                                  # TODO: the two quantiles of t_star, upper first
  unname(theta_hat - q * se_hat)
}
ci_bca <- function(theta_star, theta_hat, x) {
  z0 <- ___________                                 # TODO: Algorithm 5, step 1
  d <- x - mean(x)
  a <- ___________                                  # TODO: step 2
  z <- qnorm(c(alpha / 2, 1 - alpha / 2))
  adj <- ___________                                # TODO: step 3, both alpha_1 and alpha_2 at once
  unname(quantile(theta_star, adj))
}
ci_t <- function(x) {
  n <- length(x)
  mean(x) + c(-1, 1) * qt(1 - alpha / 2, n - 1) * sd(x) / sqrt(n)
}
ci_exact <- function(x) {
  n <- length(x)
  2 * n * mean(x) / qchisq(c(1 - alpha / 2, alpha / 2), 2 * n)
}

all_intervals <- function(x, B) {
  n <- length(x)
  M <- boot_matrix(x, B)
  theta_star <- colMeans(M)
  se_star <- apply(M, 2, sd) / sqrt(n)
  rbind(percentile    = ci_perc(theta_star),
        basic         = ci_basic(theta_star, mean(x)),
        `bootstrap-t` = ci_boott(theta_star, se_star, mean(x), sd(x) / sqrt(n)),
        BCa           = ci_bca(theta_star, mean(x), x),
        `t interval`  = ci_t(x),
        exact         = ci_exact(x))
}
ci <- all_intervals(x, B)
colnames(ci) <- c("lower", "upper")
print(round(ci, 3))
df <- data.frame(method = factor(rownames(ci), levels = rev(rownames(ci))), lower = ci[, 1], upper = ci[, 2])
p <- ggplot(df) +
  geom_vline(xintercept = 1, color = pal[["dark"]], linewidth = 1) +
  geom_vline(xintercept = theta_hat, color = pal[["grey"]], linetype = "dashed") +
  geom_segment(aes(x = lower, xend = upper, y = method, yend = method), color = pal[["blue"]], linewidth = 1.6) +
  scale_x_continuous(labels = tex_num) +
  labs(x = r"(value ($\mu = 1$ solid, $\hat\theta$ dashed))", y = NULL,
       title = r"(Six $95\%$ intervals for the mean of one $\mathrm{Exp}(1)$ sample, $n = 17$)") +
  tex_theme
knitr::include_graphics(tikz_png(p, "figures/lab05/six-intervals.png", width = 9, height = 4.5))

Report.

  1. The six intervals.
  2. The figure.
  3. Which intervals contain \(\mu = 1\).
  4. Which interval reaches furthest to the right.

Exercise 1.3: Coverage

Run Algorithm 7 with \(F = \Exp(1)\), \(\theta = 1\), \(n = 17\), and \(R = 2000\) datasets. The code chunk below is complete. It calls your all_intervals(), so all six methods run at once.

set.seed(536)
coverage_study <- function(R, n, B, make_ci, theta) {
  hit <- below <- above <- NULL
  for (r in seq_len(R)) {
    x_r <- rexp(n)
    ci_r <- make_ci(x_r, B)
    if (is.null(hit)) hit <- below <- above <- setNames(numeric(nrow(ci_r)), rownames(ci_r))
    hit   <- hit   + (ci_r[, 1] <= theta & theta <= ci_r[, 2])
    below <- below + (theta < ci_r[, 1])
    above <- above + (theta > ci_r[, 2])
  }
  data.frame(method = names(hit), coverage = hit / R, mc_se = sqrt(hit / R * (1 - hit / R) / R),
             theta_below_L = below / R, theta_above_U = above / R, row.names = NULL)
}
cov1 <- coverage_study(R = 2000, n = 17, B = B, make_ci = all_intervals, theta = 1)
knitr::kable(cov1, digits = 3)

Report.

  1. The coverage table.
  2. Which methods cover within Monte Carlo error of \(0.95\).
  3. How each bootstrap interval compares with the \(t\) interval. That is the interval you would report without a bootstrap.
  4. On which side the percentile and basic intervals miss.

Part 2: Where the Bootstrap Fails, and a Repair

\(X_1, \ldots, X_n \iid \Uniform(0, \theta)\) with \(\hat\theta = X_{(n)} = \max_i X_i\), the lecture’s counterexample, at \(n = 50\). We generate the data with \(\theta = 1\) and treat \(\theta\) as unknown. From the lecture: \(n(\theta - X_{(n)})/\theta \xrightarrow{d} \Exp(1)\). The nonparametric bootstrap puts mass \(1 - (1 - 1/n)^n \to 0.632\) on \(X^*_{(n)} = X_{(n)}\).

Exercise 2.1: The Spike

Bootstrap the maximum: Algorithm 1 with \(T = \max\).

  1. Fill in the bootstrap maxima.
  2. Fill in the proportion of replicates equal to \(\hat\theta\).
set.seed(536)
n <- 50
u <- runif(n)                       # theta = 1, treated as unknown
theta_hat <- max(u)
max_star <- ___________             # TODO: the B bootstrap maxima (boot_matrix, then the column maxima)
spike <- ___________                # TODO: the proportion of replicates equal to theta_hat
print(c(theta_hat = theta_hat, spike = spike, theory = 1 - (1 - 1 / n)^n))

Report.

  1. \(\hat\theta\).
  2. The proportion of replicates equal to \(\hat\theta\), next to \(1 - (1 - 1/n)^n\).

Exercise 2.2: The Parametric Bootstrap

The parametric bootstrap assumes a model for \(F\), fits it, and resamples from the fitted model (Algorithm 6). Part 1 never assumed a model. This one does, and a wrong model passes its error to every resample. Here the model is \(\Uniform(0, \theta)\) and the fit is \(\hat\theta = X_{(n)}\).

  1. Resample from \(\Uniform(0, X_{(n)})\) instead of from the data.
  2. Report the proportion of replicates equal to \(\hat\theta\).
  3. Run the figure chunk. It draws both bootstrap distributions of \(n(\hat\theta - \hat\theta^*)/\hat\theta\) against the \(\Exp(1)\) limit.
set.seed(536)
max_param <- ___________            # TODO: B resamples of size n from Uniform(0, theta_hat), the maximum of each
spike_param <- mean(max_param == theta_hat)
print(c(spike_param = spike_param))
stat_np <- n * (theta_hat - max_star) / theta_hat
stat_p  <- n * (theta_hat - max_param) / theta_hat
xmax <- max(c(stat_np, stat_p, 6))
pmf <- as.data.frame(table(stat_np) / B, stringsAsFactors = FALSE)
names(pmf) <- c("val", "prop"); pmf$val <- as.numeric(pmf$val)
p1 <- ggplot(pmf) +
  geom_segment(aes(x = val, xend = val, y = 0, yend = prop), color = pal[["blue"]], linewidth = 1.3) +
  geom_point(aes(val, prop), color = pal[["blue"]], size = 1.6) +
  annotate("text", x = xmax, y = max(pmf$prop), hjust = 1, vjust = 1, color = pal[["red"]], size = 5,
           label = sprintf(r"(spike at $0$: mass $%.3f$)", spike)) +
  coord_cartesian(xlim = c(0, xmax)) +
  labs(x = r"($n(\hat\theta - \hat\theta^*)/\hat\theta$)", y = "proportion", title = "Nonparametric bootstrap: resample the data")
p2 <- ggplot(data.frame(s = stat_p), aes(s)) +
  geom_histogram(aes(y = after_stat(density)), breaks = seq(0, xmax, length.out = 31), fill = pal[["grey"]], color = "white") +
  stat_function(fun = function(z) dexp(z, 1), color = pal[["dark"]], linewidth = 1.1) +
  coord_cartesian(xlim = c(0, xmax)) +
  labs(x = r"($n(\hat\theta - \hat\theta^*)/\hat\theta$)", y = "density",
       title = r"(Parametric bootstrap; curve: the $\mathrm{Exp}(1)$ limit)")
p <- (p1 + p2) & scale_x_continuous(labels = tex_num) & scale_y_continuous(labels = tex_num) & tex_theme &
  theme(plot.title = element_text(size = 14))
knitr::include_graphics(tikz_png(p, "figures/lab05/max-two-ways.png", width = 10, height = 4.5))

Report.

  1. The proportion of parametric replicates equal to \(\hat\theta\).
  2. The figure.
  3. One sentence: which bootstrap distribution does the \(\Exp(1)\) curve describe?

Exercise 2.3: Coverage under Both Schemes

Here \(X_{(n)}/\theta \sim \Betad(n, 1)\), with \(P(X_{(n)}/\theta \leq u) = u^n\), so \[P\Bigl((\alpha/2)^{1/n} \leq X_{(n)}/\theta \leq (1 - \alpha/2)^{1/n}\Bigr) = 1 - \alpha \qquad\Longrightarrow\qquad \Bigl[\,X_{(n)} \big/ (1 - \alpha/2)^{1/n},\ \ X_{(n)} \big/ (\alpha/2)^{1/n}\,\Bigr]\] is an exact interval. Run Algorithm 7 with \(R = 2000\) for five intervals: percentile and basic under nonparametric resampling, percentile and basic under parametric resampling, and the exact interval. The code chunk below is complete.

set.seed(536)
max_intervals <- function(u, B) {
  n <- length(u); theta_hat <- max(u)
  a <- apply(boot_matrix(u, B), 2, max)                                   # nonparametric replicates
  p <- apply(matrix(runif(n * B, 0, theta_hat), nrow = n), 2, max)       # parametric replicates
  rbind(`percentile, nonparametric` = ci_perc(a),
        `basic, nonparametric`      = ci_basic(a, theta_hat),
        `percentile, parametric`    = ci_perc(p),
        `basic, parametric`         = ci_basic(p, theta_hat),
        exact                       = theta_hat / c((1 - alpha / 2)^(1 / n), (alpha / 2)^(1 / n)))
}
coverage_max <- function(R, n, B, theta = 1) {
  hit <- below <- above <- NULL
  for (r in seq_len(R)) {
    ci_r <- max_intervals(runif(n, 0, theta), B)
    if (is.null(hit)) hit <- below <- above <- setNames(numeric(nrow(ci_r)), rownames(ci_r))
    hit   <- hit   + (ci_r[, 1] <= theta & theta <= ci_r[, 2])
    below <- below + (theta < ci_r[, 1])
    above <- above + (theta > ci_r[, 2])
  }
  data.frame(method = names(hit), coverage = hit / R, mc_se = sqrt(hit / R * (1 - hit / R) / R),
             theta_below_L = below / R, theta_above_U = above / R, row.names = NULL)
}
cov2 <- coverage_max(R = 2000, n = 50, B = B)
knitr::kable(cov2, digits = 3)

Report.

  1. The table.
  2. One sentence per row on why its coverage is what it is.
  3. Which of the two failures the parametric bootstrap repairs, and which it does not.

Exercise 2.4: What Fails Means

On the one sample of Exercise 2.2 the basic intervals under both schemes were within \(0.01\) of the exact interval. In Part 1 the bootstrap intervals differed visibly from the exact one. One sample is not the criterion. Coverage is, and how it moves with \(n\). The code chunk below reruns both coverage studies at four times the sample size, with \(R = 500\).

set.seed(536)
cov1_big <- coverage_study(R = 500, n = 68, B = B, make_ci = all_intervals, theta = 1)
cov2_big <- coverage_max(R = 500, n = 200, B = B)
knitr::kable(cov1_big, digits = 3)
knitr::kable(cov2_big, digits = 3)
print(c(spike_theory_200 = 1 - (1 - 1 / 200)^200))

Report.

  1. The two tables.
  2. One sentence on what separates a bootstrap that fails from one that is imprecise at small \(n\).

Summary and Reflection

Takeaways

  1. One resampling primitive, four intervals. All four bootstrap intervals are functions of the same \(B\) replicates. The bootstrap-t also needs \(\widehat{\text{SE}}^{*(b)}\) on each resample, the same rule applied to resample \(b\).
  2. Coverage is the judge. Nominal \(95\%\) is a promise about the procedure. What counts is actual coverage at the \(n\) you have, on each side. Simulation is one way to estimate it, with standard error \(\sqrt{\hat C(1 - \hat C)/R}\).
  3. The benchmark is the \(t\) interval. A bootstrap interval earns its place by beating the interval you would report without it, at the sample size you have.
  4. Reflection has a direction. The basic interval is the percentile interval reflected through \(\hat\theta\). Whether that helps depends on which way the sampling distribution is skewed.
  5. What you resample determines what you assume. The nonparametric bootstrap draws from \(\hat F_n\). The parametric bootstrap draws from the fitted model. Each failure of the bootstrap traces back to what the resampling can and cannot produce.
  6. Fails means more data does not help. A consistent bootstrap’s coverage error shrinks with \(n\). An inconsistent one’s does not.

Self-Assessment Questions

  1. Why does the bootstrap-t need a standard error on each resample, and what happens if you use the data’s \(\widehat{\text{SE}}\) for every \(t^{*(b)}\)?
  2. Two intervals both cover \(0.90\) of the time at \(n = 17\). One misses \(0.05\) on each side, the other \(0.01\) below and \(0.09\) above. Which would you rather report a lower confidence limit from, and why?
  3. With \(R = 2000\), can you tell a coverage of \(0.90\) from \(0.91\)? From \(0.87\)? Give the standard error.
  4. For the uniform maximum, what bounds every bootstrap replicate of \(\hat\theta\), under either resampling scheme, and what does that imply for the percentile interval?
  5. State one situation in which you would trust the parametric bootstrap over the nonparametric one, and one in which you would not.
  6. Two bootstrap intervals both cover \(0.87\) at \(n = 50\). What one further simulation tells you which of the two, if either, is failing?