Purpose of this practical aid

This RPubs-ready document develops three central computational ideas in advanced data analysis:

  1. Bootstrap: approximate a sampling distribution by repeated sampling with replacement from the empirical distribution.
  2. Jackknife: study bias, standard error, and influence by systematically deleting observations.
  3. Expectation-Maximization: maximize an observed-data likelihood by alternating conditional expectation and maximization.

Every major idea is demonstrated through simulation, real data, plots, tables, diagnostics, and student tasks.

1 Learning outcomes

After completing the practicals, a student should be able to:

  1. implement nonparametric, parametric, residual, and pairs bootstrap methods;
  2. estimate standard errors, bias, and confidence intervals;
  3. implement the delete-one jackknife from first principles;
  4. use jackknife replicates and pseudovalues for influence analysis;
  5. derive and code EM updates for missing or latent data;
  6. monitor EM convergence and recognize local maxima or degeneracy;
  7. compare methods through simulation rather than relying only on formulas.

2 Part I: Bootstrap

2.1 1. Bootstrap principle

Let \(X_1,\ldots,X_n\) be an observed random sample and let

\[ \widehat\theta = T(X_1,\ldots,X_n) \]

estimate a population parameter \(\theta\). The empirical distribution is

\[ \widehat F_n(x)=\frac{1}{n}\sum_{i=1}^n I(X_i\le x). \]

A nonparametric bootstrap sample is obtained by drawing

\[ X_1^*,\ldots,X_n^* \overset{\text{i.i.d.}}{\sim} \widehat F_n, \]

which is equivalent to sampling the observed data with replacement. The bootstrap statistic is

\[ \widehat\theta^*=T(X_1^*,\ldots,X_n^*). \]

Repeating this \(B\) times produces

\[ \widehat\theta_1^*,\ldots,\widehat\theta_B^*. \]

The bootstrap standard error is

\[ \widehat{\mathrm{SE}}_{\mathrm{boot}}(\widehat\theta) = \left[ \frac{1}{B-1} \sum_{b=1}^{B} \left(\widehat\theta_b^*-\overline{\theta}^*\right)^2 \right]^{1/2}. \]

bootstrap_1d <- function(x, statistic, B = 3000, seed = 2026, ...) {
  set.seed(seed)
  n <- length(x)
  theta_hat <- statistic(x, ...)
  theta_star <- replicate(
    B,
    statistic(sample(x, size = n, replace = TRUE), ...)
  )

  list(
    estimate = theta_hat,
    replicates = theta_star,
    bias = mean(theta_star) - theta_hat,
    se = sd(theta_star),
    percentile_ci = unname(quantile(theta_star, c(0.025, 0.975))),
    basic_ci = 2 * theta_hat -
      rev(unname(quantile(theta_star, c(0.025, 0.975))))
  )
}

2.2 2. Simulation 1: Bootstrap standard error of the mean

set.seed(12)
x_normal <- rnorm(40, mean = 10, sd = 3)

boot_mean <- bootstrap_1d(x_normal, mean, B = 5000)

mean_summary <- tibble(
  Quantity = c("Observed mean", "Bootstrap bias",
               "Bootstrap SE", "Formula-based SE",
               "Percentile CI lower", "Percentile CI upper"),
  Value = c(
    boot_mean$estimate,
    boot_mean$bias,
    boot_mean$se,
    sd(x_normal) / sqrt(length(x_normal)),
    boot_mean$percentile_ci
  )
)

pretty_table(mean_summary, caption =
  "Bootstrap analysis of the mean for simulated normal data")
Bootstrap analysis of the mean for simulated normal data
Quantity Value
Observed mean 9.6497
Bootstrap bias 0.0046
Bootstrap SE 0.4165
Formula-based SE 0.4227
Percentile CI lower 8.8453
Percentile CI upper 10.4977 
tibble(value = boot_mean$replicates) |>
  ggplot(aes(value)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 38, fill = "#5B8FF9",
                 color = "white", alpha = 0.9) +
  geom_density(color = "#D62728", linewidth = 1.2) +
  geom_vline(xintercept = boot_mean$estimate,
             color = "#111827", linewidth = 1,
             linetype = 2) +
  labs(
    title = "Bootstrap distribution of the sample mean",
    subtitle = "Simulated N(10, 3²) sample, n = 40, B = 5,000",
    x = expression(bar(X)^"*"),
    y = "Density"
  )

2.3 3. Simulation 2: Mean versus median under skewness

set.seed(15)
x_exp <- rexp(35, rate = 0.5)

boot_exp_mean <- bootstrap_1d(x_exp, mean, B = 5000)
boot_exp_median <- bootstrap_1d(x_exp, median, B = 5000)

boot_compare <- bind_rows(
  tibble(Statistic = "Mean",
         Estimate = boot_exp_mean$estimate,
         Bias = boot_exp_mean$bias,
         SE = boot_exp_mean$se,
         Lower = boot_exp_mean$percentile_ci[1],
         Upper = boot_exp_mean$percentile_ci[2]),
  tibble(Statistic = "Median",
         Estimate = boot_exp_median$estimate,
         Bias = boot_exp_median$bias,
         SE = boot_exp_median$se,
         Lower = boot_exp_median$percentile_ci[1],
         Upper = boot_exp_median$percentile_ci[2])
)

pretty_table(boot_compare,
             caption = "Bootstrap comparison under exponential skewness")
Bootstrap comparison under exponential skewness
Statistic Estimate Bias SE Lower Upper
Mean 1.8305 0.0087 0.3345 1.2429 2.5404
Median 1.2149 0.0008 0.1962 0.8265 1.5745 
bind_rows(
  tibble(value = boot_exp_mean$replicates, Statistic = "Mean"),
  tibble(value = boot_exp_median$replicates, Statistic = "Median")
) |>
  ggplot(aes(value, fill = Statistic)) +
  geom_density(alpha = 0.55) +
  facet_wrap(~ Statistic, scales = "free") +
  scale_fill_manual(values = c("Mean" = "#FF9D4D",
                               "Median" = "#36B37E")) +
  labs(
    title = "Bootstrap distributions under severe skewness",
    subtitle = "The mean and median have visibly different uncertainty",
    x = "Bootstrap statistic",
    y = "Density"
  ) +
  guides(fill = "none")

2.4 4. Bootstrap confidence intervals

The most common intervals are:

Normal interval \[ \widehat\theta \pm z_{1-\alpha/2} \widehat{\mathrm{SE}}_{\mathrm{boot}}. \]

Percentile interval \[ \left[ \widehat\theta_{(\alpha/2)}^*, \widehat\theta_{(1-\alpha/2)}^* \right]. \]

Basic interval \[ \left[ 2\widehat\theta-\widehat\theta_{(1-\alpha/2)}^*, 2\widehat\theta-\widehat\theta_{(\alpha/2)}^* \right]. \]

normal_ci <- boot_exp_mean$estimate +
  qnorm(c(0.025, 0.975)) * boot_exp_mean$se

ci_table <- tibble(
  Method = c("Normal", "Percentile", "Basic"),
  Lower = c(normal_ci[1],
            boot_exp_mean$percentile_ci[1],
            boot_exp_mean$basic_ci[1]),
  Upper = c(normal_ci[2],
            boot_exp_mean$percentile_ci[2],
            boot_exp_mean$basic_ci[2])
)

pretty_table(ci_table,
             caption = "Three bootstrap confidence intervals for the mean")
Three bootstrap confidence intervals for the mean
Method Lower Upper
Normal 1.1749 2.4861
Percentile 1.2429 2.5404
Basic 1.1206 2.4181 
ci_table |>
  mutate(Method = factor(Method, levels = rev(Method))) |>
  ggplot(aes(y = Method, x = Lower, xend = Upper,
             yend = Method, color = Method)) +
  geom_segment(linewidth = 3, lineend = "round") +
  geom_point(aes(x = Lower), size = 3) +
  geom_point(aes(x = Upper), size = 3) +
  geom_vline(xintercept = boot_exp_mean$estimate,
             linetype = 2, color = "grey25") +
  scale_color_manual(values = c(
    "Normal" = "#5B8FF9",
    "Percentile" = "#36B37E",
    "Basic" = "#FF6B6B"
  )) +
  labs(
    title = "Bootstrap confidence-interval comparison",
    x = "Parameter value",
    y = NULL
  ) +
  guides(color = "none")

2.5 5. Simulation 3: Coverage study

This simulation checks how often nominal 95% intervals cover the true exponential mean \(2\).

set.seed(25)

one_coverage_run <- function(n = 25, B = 800) {
  x <- rexp(n, rate = 0.5)
  b <- bootstrap_1d(x, mean, B = B, seed = sample.int(1e7, 1))

  normal <- b$estimate + qnorm(c(0.025, 0.975)) * b$se
  c(
    Normal = normal[1] <= 2 && 2 <= normal[2],
    Percentile = b$percentile_ci[1] <= 2 && 2 <= b$percentile_ci[2],
    Basic = b$basic_ci[1] <= 2 && 2 <= b$basic_ci[2]
  )
}

coverage_raw <- replicate(250, one_coverage_run())
coverage <- rowMeans(coverage_raw)

coverage_table <- tibble(
  Method = names(coverage),
  Empirical_Coverage = as.numeric(coverage),
  Nominal_Coverage = 0.95
)

pretty_table(coverage_table,
             caption = "Empirical 95% coverage in an exponential-mean experiment")
Empirical 95% coverage in an exponential-mean experiment
Method Empirical_Coverage Nominal_Coverage
Normal 0.920 0.95
Percentile 0.924 0.95
Basic 0.900 0.95 
coverage_table |>
  ggplot(aes(Method, Empirical_Coverage, fill = Method)) +
  geom_col(width = 0.68) +
  geom_hline(yintercept = 0.95, linetype = 2,
             linewidth = 1, color = "#B91C1C") +
  scale_y_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0, 1)) +
  scale_fill_manual(values = c(
    "Normal" = "#5B8FF9",
    "Percentile" = "#36B37E",
    "Basic" = "#FF9D4D"
  )) +
  labs(
    title = "Do nominal 95% intervals really cover 95%?",
    subtitle = "250 simulated samples; dashed line is nominal coverage",
    x = NULL, y = "Empirical coverage"
  ) +
  guides(fill = "none")

2.6 6. Real data 1: Bootstrap correlation in mtcars

cor_stat <- function(data, indices) {
  d <- data[indices, ]
  cor(d$mpg, d$wt)
}

set.seed(102)
boot_mtcars <- boot(mtcars, statistic = cor_stat, R = 5000)

mtcars_boot_summary <- tibble(
  Quantity = c("Observed correlation",
               "Bootstrap bias",
               "Bootstrap SE"),
  Value = c(boot_mtcars$t0,
            mean(boot_mtcars$t) - boot_mtcars$t0,
            sd(boot_mtcars$t))
)

pretty_table(mtcars_boot_summary,
             caption = "Bootstrap analysis of correlation: mpg versus wt")
Bootstrap analysis of correlation: mpg versus wt
Quantity Value
Observed correlation -0.8677
Bootstrap bias -0.0028
Bootstrap SE 0.0340 
tibble(correlation = as.numeric(boot_mtcars$t)) |>
  ggplot(aes(correlation)) +
  geom_histogram(bins = 35, fill = "#7C3AED",
                 color = "white", alpha = 0.9) +
  geom_vline(xintercept = boot_mtcars$t0,
             linetype = 2, linewidth = 1.1) +
  labs(
    title = "Bootstrap distribution of correlation",
    subtitle = "Real data: mtcars, mpg versus wt",
    x = "Bootstrap correlation",
    y = "Frequency"
  )

2.7 7. Real data 2: Pairs bootstrap for regression

slope_stat <- function(data, indices) {
  coef(lm(mpg ~ wt + hp, data = data[indices, ]))["wt"]
}

set.seed(103)
boot_slope <- boot(mtcars, statistic = slope_stat, R = 5000)

fit_mtcars <- lm(mpg ~ wt + hp, data = mtcars)

reg_table <- tibble(
  Method = c("Ordinary least squares",
             "Pairs bootstrap"),
  Estimate = c(coef(fit_mtcars)["wt"], boot_slope$t0),
  Standard_Error = c(
    summary(fit_mtcars)$coef["wt", "Std. Error"],
    sd(boot_slope$t)
  )
)

pretty_table(reg_table,
             caption = "Model-based and pairs-bootstrap uncertainty for wt")
Model-based and pairs-bootstrap uncertainty for wt
Method Estimate Standard_Error
Ordinary least squares -3.8778 0.6327
Pairs bootstrap -3.8778 0.7111 
tibble(slope = as.numeric(boot_slope$t)) |>
  ggplot(aes(slope)) +
  geom_density(fill = "#22C55E", alpha = 0.55,
               color = "#166534", linewidth = 1.1) +
  geom_vline(xintercept = boot_slope$t0,
             linetype = 2, linewidth = 1) +
  labs(
    title = "Pairs-bootstrap distribution of the weight coefficient",
    subtitle = "Model: mpg ~ wt + hp",
    x = "Bootstrap coefficient of wt",
    y = "Density"
  )

2.8 8. Residual bootstrap for regression

residual_bootstrap <- function(model, B = 3000, seed = 2026) {
  set.seed(seed)
  yhat <- fitted(model)
  e <- resid(model) - mean(resid(model))
  dat <- model.frame(model)
  response_name <- names(dat)[1]

  replicate(B, {
    y_star <- yhat + sample(e, replace = TRUE)
    dat[[response_name]] <- y_star
    coef(update(model, data = dat))
  })
}

resid_boot <- residual_bootstrap(fit_mtcars, B = 4000)

coef_compare <- tibble(
  Coefficient = rownames(t(resid_boot)),
  Estimate = coef(fit_mtcars),
  Residual_Bootstrap_SE = apply(resid_boot, 1, sd)
)

pretty_table(coef_compare,
             caption = "Residual-bootstrap standard errors")
Residual-bootstrap standard errors
Estimate Residual_Bootstrap_SE
37.2273 1.5282
-3.8778 0.5978
-0.0318 0.0085

3 Part II: Jackknife

3.1 9. Definition and implementation

For the estimator

\[ \widehat\theta=T(X_1,\ldots,X_n), \]

the \(i\)-th delete-one replicate is

\[ \widehat\theta_{(-i)} = T(X_1,\ldots,X_{i-1},X_{i+1},\ldots,X_n). \]

Let

\[ \overline\theta_{(-\cdot)} = \frac1n\sum_{i=1}^n\widehat\theta_{(-i)}. \]

The jackknife bias estimate is

\[ \widehat{\mathrm{Bias}}_{\mathrm{jack}} = (n-1) \left( \overline\theta_{(-\cdot)}-\widehat\theta \right). \]

The bias-corrected jackknife estimator is

\[ \widehat\theta_{\mathrm{jack}} = n\widehat\theta-(n-1)\overline\theta_{(-\cdot)}. \]

The jackknife variance estimate is

\[ \widehat{\mathrm{Var}}_{\mathrm{jack}} = \frac{n-1}{n} \sum_{i=1}^n \left( \widehat\theta_{(-i)} -\overline\theta_{(-\cdot)} \right)^2. \]

jackknife <- function(data, statistic, ...) {
  n <- if (is.data.frame(data) || is.matrix(data)) nrow(data) else length(data)

  theta_hat <- statistic(data, ...)

  theta_delete <- map_dbl(seq_len(n), function(i) {
    reduced <- if (is.data.frame(data) || is.matrix(data)) {
      data[-i, , drop = FALSE]
    } else {
      data[-i]
    }
    statistic(reduced, ...)
  })

  theta_bar <- mean(theta_delete)
  bias <- (n - 1) * (theta_bar - theta_hat)
  variance <- (n - 1) / n *
    sum((theta_delete - theta_bar)^2)
  pseudovalues <- n * theta_hat - (n - 1) * theta_delete
  influence <- (n - 1) * (theta_hat - theta_delete)

  list(
    estimate = theta_hat,
    delete_one = theta_delete,
    replicate_mean = theta_bar,
    bias = bias,
    bias_corrected = theta_hat - bias,
    variance = variance,
    se = sqrt(variance),
    pseudovalues = pseudovalues,
    influence = influence
  )
}

3.2 10. Simulation 4: Exact jackknife SE for a mean

set.seed(206)
x_jack <- rnorm(30, 5, 2)

jk_mean <- jackknife(x_jack, mean)

jack_mean_table <- tibble(
  Quantity = c("Mean", "Jackknife bias",
               "Jackknife SE", "Usual SE"),
  Value = c(jk_mean$estimate,
            jk_mean$bias,
            jk_mean$se,
            sd(x_jack) / sqrt(length(x_jack)))
)

pretty_table(jack_mean_table,
             caption = "Jackknife and usual standard error for a mean")
Jackknife and usual standard error for a mean
Quantity Value
Mean 5.2081
Jackknife bias 0.0000
Jackknife SE 0.2912
Usual SE 0.2912 
tibble(
  Deleted = seq_along(jk_mean$delete_one),
  Replicate = jk_mean$delete_one
) |>
  ggplot(aes(Deleted, Replicate)) +
  geom_line(color = "#2563EB", linewidth = 0.8) +
  geom_point(aes(color = abs(Replicate - mean(Replicate))),
             size = 2.8) +
  geom_hline(yintercept = jk_mean$estimate,
             linetype = 2) +
  scale_color_gradient(low = "#93C5FD", high = "#DC2626") +
  labs(
    title = "Delete-one means",
    subtitle = "Points farther from the dashed line have greater influence",
    x = "Deleted observation",
    y = "Leave-one-out mean",
    color = "Absolute change"
  )

3.3 11. Simulation 5: Bias correction for variance

biased_variance <- function(z) mean((z - mean(z))^2)

set.seed(207)
x_var <- rnorm(18, 0, 3)
jk_var <- jackknife(x_var, biased_variance)

variance_table <- tibble(
  Estimator = c("Biased variance, divisor n",
                "Jackknife bias-corrected",
                "Usual unbiased variance"),
  Value = c(jk_var$estimate,
            jk_var$bias_corrected,
            var(x_var))
)

pretty_table(variance_table,
             caption = "Jackknife correction recovers the unbiased variance")
Jackknife correction recovers the unbiased variance
Estimator Value
Biased variance, divisor n 6.8551
Jackknife bias-corrected 7.2583
Usual unbiased variance 7.2583

3.4 12. Simulation 6: Ratio estimator and influential observations

set.seed(208)
ratio_data <- tibble(
  x = rgamma(60, shape = 5, rate = 1),
  y = 2.4 * x + rnorm(60, sd = 2)
)

ratio_stat <- function(d) mean(d$y) / mean(d$x)
jk_ratio <- jackknife(ratio_data, ratio_stat)

ratio_summary <- tibble(
  Quantity = c("Original estimate", "Estimated bias",
               "Bias-corrected estimate", "Jackknife SE"),
  Value = c(jk_ratio$estimate, jk_ratio$bias,
            jk_ratio$bias_corrected, jk_ratio$se)
)

pretty_table(ratio_summary,
             caption = "Jackknife analysis of a ratio-of-means estimator")
Jackknife analysis of a ratio-of-means estimator
Quantity Value
Original estimate 2.3860
Estimated bias -0.0006
Bias-corrected estimate 2.3866
Jackknife SE 0.0538 
tibble(
  Observation = seq_along(jk_ratio$influence),
  Influence = jk_ratio$influence
) |>
  ggplot(aes(Observation, Influence,
             fill = abs(Influence))) +
  geom_col() +
  scale_fill_gradient(low = "#A7F3D0", high = "#B91C1C") +
  geom_hline(yintercept = 0, linewidth = 0.6) +
  labs(
    title = "Jackknife influence values",
    subtitle = "Ratio-of-means simulation",
    x = "Observation",
    y = "Empirical influence",
    fill = "|Influence|"
  )

3.5 13. Real data 3: Jackknife correlation for iris

iris2 <- iris |>
  select(Sepal.Length, Petal.Length)

iris_cor <- function(d) cor(d$Sepal.Length, d$Petal.Length)
jk_iris <- jackknife(iris2, iris_cor)

iris_summary <- tibble(
  Quantity = c("Correlation", "Bias",
               "Bias-corrected correlation", "SE"),
  Value = c(jk_iris$estimate, jk_iris$bias,
            jk_iris$bias_corrected, jk_iris$se)
)

pretty_table(iris_summary,
             caption = "Jackknife correlation analysis for iris")
Jackknife correlation analysis for iris
Quantity Value
Correlation 0.8718
Bias -0.0001
Bias-corrected correlation 0.8719
SE 0.0173 
iris_influence <- iris |>
  mutate(
    Row = row_number(),
    Influence = jk_iris$influence
  )

top_iris <- iris_influence |>
  slice_max(abs(Influence), n = 8) |>
  select(Row, Species, Sepal.Length, Petal.Length, Influence)

pretty_table(top_iris,
             caption = "Most influential iris observations")
Most influential iris observations
Row Species Sepal.Length Petal.Length Influence
107 virginica 4.9 4.5 -1.1356
15 setosa 5.8 1.2 -0.8556
16 setosa 5.7 1.5 -0.5130
119 virginica 7.7 6.9 0.4308
85 versicolor 5.4 4.5 -0.4305
122 virginica 5.6 4.9 -0.4137
19 setosa 5.7 1.7 -0.4096
114 virginica 5.7 5.0 -0.3537 
iris_influence |>
  ggplot(aes(Row, Influence, color = Species)) +
  geom_hline(yintercept = 0, color = "grey50") +
  geom_point(size = 2.2, alpha = 0.85) +
  scale_color_manual(values = c(
    "setosa" = "#5B8FF9",
    "versicolor" = "#F6BD16",
    "virginica" = "#E8684A"
  )) +
  labs(
    title = "Observation-level influence on iris correlation",
    x = "Row",
    y = "Jackknife influence"
  )

3.6 14. Real data 4: Jackknife regression influence

wt_slope <- function(d)
  unname(coef(lm(mpg ~ wt + hp, data = d))["wt"])

jk_reg <- jackknife(mtcars, wt_slope)

reg_jk_table <- tibble(
  Quantity = c("Slope estimate", "Jackknife bias",
               "Bias-corrected slope", "Jackknife SE",
               "Model-based SE"),
  Value = c(
    jk_reg$estimate,
    jk_reg$bias,
    jk_reg$bias_corrected,
    jk_reg$se,
    summary(fit_mtcars)$coef["wt", "Std. Error"]
  )
)

pretty_table(reg_jk_table,
             caption = "Jackknife uncertainty for the wt coefficient")
Jackknife uncertainty for the wt coefficient
Quantity Value
Slope estimate -3.8778
Jackknife bias 0.0381
Bias-corrected slope -3.9160
Jackknife SE 0.7564
Model-based SE 0.6327 
tibble(
  Car = rownames(mtcars),
  Influence = jk_reg$influence
) |>
  mutate(Car = reorder(Car, Influence)) |>
  ggplot(aes(Car, Influence, fill = Influence > 0)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = c("#EF4444", "#3B82F6")) +
  labs(
    title = "Which cars influence the weight coefficient?",
    x = NULL,
    y = "Jackknife influence",
    fill = "Positive influence"
  )

3.7 15. Jackknife failure for a nonsmooth statistic

x_med <- c(1, 2, 3, 3, 3, 4, 5)

jk_med <- jackknife(x_med, median)
boot_med <- bootstrap_1d(x_med, median, B = 5000)

failure_table <- tibble(
  Method = c("Delete-one jackknife", "Nonparametric bootstrap"),
  Standard_Error = c(jk_med$se, boot_med$se)
)

pretty_table(failure_table,
             caption = "Median example: the delete-one jackknife can fail")
Median example: the delete-one jackknife can fail
Method Standard_Error
Delete-one jackknife 0.0000
Nonparametric bootstrap 0.5416 
tibble(median = boot_med$replicates) |>
  ggplot(aes(median)) +
  geom_bar(fill = "#EC4899", color = "white") +
  labs(
    title = "Bootstrap distribution of the median",
    subtitle = "Jackknife SE is zero, but bootstrap variation is positive",
    x = "Bootstrap median",
    y = "Frequency"
  )

4 Part III: Expectation-Maximization Algorithm

4.1 16. Missing-data likelihood framework

Let \(Y\) be observed data and \(Z\) be missing or latent data. The complete-data likelihood is

\[ L_c(\theta;Y,Z) \]

and the observed-data likelihood is

\[ L(\theta;Y) = \sum_z L_c(\theta;Y,z) \]

for discrete \(Z\), or the corresponding integral for continuous \(Z\).

At iteration \(t\), the EM algorithm performs:

E-step \[ Q(\theta\mid\theta^{(t)}) = E_{\theta^{(t)}} \left[ \log L_c(\theta;Y,Z)\mid Y \right]. \]

M-step \[ \theta^{(t+1)} = \arg\max_\theta Q(\theta\mid\theta^{(t)}). \]

Under standard conditions,

\[ \ell(\theta^{(t+1)};Y) \ge \ell(\theta^{(t)};Y). \]

4.2 17. Simulation 7: Missing Bernoulli observations

Suppose \(X_i\sim\mathrm{Bernoulli}(p)\), but some outcomes are missing. If the number of observed successes is \(s\), the number of observed failures is \(f\), and \(m\) values are missing, then:

\[ E(S_{\mathrm{mis}}\mid\text{observed},p^{(t)})=m p^{(t)}. \]

The M-step is

\[ p^{(t+1)} = \frac{s+mp^{(t)}}{s+f+m}. \]

em_bernoulli <- function(successes, failures, missing,
                         p0 = 0.5, tol = 1e-10,
                         max_iter = 1000) {
  p <- p0
  history <- tibble(iteration = 0, p = p)

  for (iter in seq_len(max_iter)) {
    expected_missing_successes <- missing * p
    p_new <- (successes + expected_missing_successes) /
      (successes + failures + missing)

    history <- bind_rows(
      history,
      tibble(iteration = iter, p = p_new)
    )

    if (abs(p_new - p) < tol) break
    p <- p_new
  }

  list(estimate = p_new, history = history,
       iterations = iter)
}

bern_em <- em_bernoulli(
  successes = 38,
  failures = 42,
  missing = 20,
  p0 = 0.15
)

pretty_table(
  tibble(
    EM_Estimate = bern_em$estimate,
    Iterations = bern_em$iterations,
    Complete_case_Estimate = 38 / (38 + 42)
  ),
  caption = "EM for Bernoulli data with missing outcomes"
)
EM for Bernoulli data with missing outcomes
EM_Estimate Iterations Complete_case_Estimate
0.475 15 0.475 
bern_em$history |>
  ggplot(aes(iteration, p)) +
  geom_line(color = "#7C3AED", linewidth = 1.2) +
  geom_point(color = "#7C3AED", size = 2) +
  labs(
    title = "EM convergence for a missing Bernoulli problem",
    x = "Iteration",
    y = expression(p^{(t)})
  )

4.3 18. Simulation 8: Two-component Gaussian mixture

Assume

\[ f(x) = \pi\phi(x;\mu_1,\sigma_1^2) + (1-\pi)\phi(x;\mu_2,\sigma_2^2). \]

Introduce latent indicators

\[ Z_i= \begin{cases} 1,&X_i\text{ belongs to component 1},\\ 0,&X_i\text{ belongs to component 2}. \end{cases} \]

The E-step responsibilities are

\[ \tau_i^{(t)} = P(Z_i=1\mid X_i,\theta^{(t)}) = \frac{ \pi^{(t)} \phi(X_i;\mu_1^{(t)},\sigma_1^{2(t)}) }{ \pi^{(t)} \phi(X_i;\mu_1^{(t)},\sigma_1^{2(t)}) + (1-\pi^{(t)}) \phi(X_i;\mu_2^{(t)},\sigma_2^{2(t)}) }. \]

The M-step updates are

\[ \pi^{(t+1)} = \frac1n\sum_i\tau_i^{(t)}, \]

\[ \mu_1^{(t+1)} = \frac{\sum_i\tau_i^{(t)}X_i} {\sum_i\tau_i^{(t)}}, \qquad \mu_2^{(t+1)} = \frac{\sum_i(1-\tau_i^{(t)})X_i} {\sum_i(1-\tau_i^{(t)})}, \]

with analogous weighted variance updates.

em_gaussian2 <- function(x,
                         pi0 = 0.5,
                         mu0 = quantile(x, c(0.3, 0.7)),
                         sigma0 = rep(sd(x), 2),
                         tol = 1e-8,
                         max_iter = 1000) {
  pi1 <- pi0
  mu1 <- mu0[1]
  mu2 <- mu0[2]
  s1 <- sigma0[1]
  s2 <- sigma0[2]

  loglik_history <- numeric()

  for (iter in seq_len(max_iter)) {
    # E-step
    d1 <- pi1 * dnorm(x, mu1, s1)
    d2 <- (1 - pi1) * dnorm(x, mu2, s2)
    tau <- d1 / pmax(d1 + d2, .Machine$double.eps)

    # M-step
    pi_new <- mean(tau)
    mu1_new <- sum(tau * x) / sum(tau)
    mu2_new <- sum((1 - tau) * x) / sum(1 - tau)
    s1_new <- sqrt(sum(tau * (x - mu1_new)^2) / sum(tau))
    s2_new <- sqrt(sum((1 - tau) * (x - mu2_new)^2) /
                     sum(1 - tau))

    ll <- sum(log(
      pi_new * dnorm(x, mu1_new, s1_new) +
        (1 - pi_new) * dnorm(x, mu2_new, s2_new)
    ))
    loglik_history <- c(loglik_history, ll)

    change <- max(abs(c(
      pi_new - pi1,
      mu1_new - mu1,
      mu2_new - mu2,
      s1_new - s1,
      s2_new - s2
    )))

    pi1 <- pi_new
    mu1 <- mu1_new
    mu2 <- mu2_new
    s1 <- max(s1_new, 1e-6)
    s2 <- max(s2_new, 1e-6)

    if (change < tol) break
  }

  list(
    pi = pi1,
    mu = c(mu1, mu2),
    sigma = c(s1, s2),
    responsibility = tau,
    loglik = loglik_history,
    iterations = iter
  )
}
set.seed(302)
n_mix <- 500
z_true <- rbinom(n_mix, 1, 0.38)
x_mix <- ifelse(
  z_true == 1,
  rnorm(n_mix, -2.2, 0.8),
  rnorm(n_mix, 2.3, 1.15)
)

mix_fit <- em_gaussian2(
  x_mix,
  pi0 = 0.5,
  mu0 = c(-1, 1),
  sigma0 = c(1.5, 1.5)
)

mix_table <- tibble(
  Parameter = c("Mixing probability",
                "Mean 1", "Mean 2",
                "SD 1", "SD 2",
                "Iterations"),
  Estimate = c(mix_fit$pi, mix_fit$mu,
               mix_fit$sigma, mix_fit$iterations)
)

pretty_table(mix_table,
             caption = "EM estimates for a simulated Gaussian mixture")
EM estimates for a simulated Gaussian mixture
Parameter Estimate
Mixing probability 0.3848
Mean 1 -2.2233
Mean 2 2.2590
SD 1 0.7384
SD 2 1.1927
Iterations 26.0000 
grid_x <- seq(min(x_mix) - 1, max(x_mix) + 1,
              length.out = 600)

density_df <- tibble(
  x = grid_x,
  Component_1 =
    mix_fit$pi *
    dnorm(grid_x, mix_fit$mu[1], mix_fit$sigma[1]),
  Component_2 =
    (1 - mix_fit$pi) *
    dnorm(grid_x, mix_fit$mu[2], mix_fit$sigma[2])
) |>
  mutate(Mixture = Component_1 + Component_2)

ggplot(tibble(x = x_mix), aes(x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 42, fill = "#CBD5E1",
                 color = "white") +
  geom_line(data = density_df,
            aes(x, Component_1, color = "Component 1"),
            linewidth = 1.2) +
  geom_line(data = density_df,
            aes(x, Component_2, color = "Component 2"),
            linewidth = 1.2) +
  geom_line(data = density_df,
            aes(x, Mixture, color = "Fitted mixture"),
            linewidth = 1.5) +
  scale_color_manual(values = c(
    "Component 1" = "#2563EB",
    "Component 2" = "#F97316",
    "Fitted mixture" = "#111827"
  )) +
  labs(
    title = "Gaussian-mixture fit obtained by EM",
    x = "Observed value",
    y = "Density",
    color = NULL
  )

tibble(
  Iteration = seq_along(mix_fit$loglik),
  LogLikelihood = mix_fit$loglik
) |>
  ggplot(aes(Iteration, LogLikelihood)) +
  geom_line(color = "#059669", linewidth = 1.25) +
  geom_point(color = "#059669", size = 1.7) +
  labs(
    title = "Observed-data log-likelihood across EM iterations",
    subtitle = "The EM monotonicity property is visible",
    y = "Log-likelihood"
  )

4.4 19. Multiple starting values and local optima

starts <- tibble(
  start = 1:12,
  pi0 = runif(12, 0.15, 0.85),
  mu1 = runif(12, min(x_mix), median(x_mix)),
  mu2 = runif(12, median(x_mix), max(x_mix)),
  sd1 = runif(12, 0.5, 2.5),
  sd2 = runif(12, 0.5, 2.5)
)

multi_fit <- pmap_dfr(starts, function(start, pi0, mu1, mu2, sd1, sd2) {
  fit <- em_gaussian2(
    x_mix,
    pi0 = pi0,
    mu0 = c(mu1, mu2),
    sigma0 = c(sd1, sd2)
  )
  tibble(
    start = start,
    final_loglik = tail(fit$loglik, 1),
    iterations = fit$iterations,
    pi = fit$pi,
    mu1 = fit$mu[1],
    mu2 = fit$mu[2]
  )
})

pretty_table(multi_fit,
             caption = "EM results from multiple initial values")
EM results from multiple initial values
start final_loglik iterations pi mu1 mu2
1 -1025.453 21 0.3848 -2.2233 2.259
2 -1025.453 44 0.3848 -2.2233 2.259
3 -1025.453 39 0.3848 -2.2233 2.259
4 -1025.453 45 0.3848 -2.2233 2.259
5 -1025.453 30 0.3848 -2.2233 2.259
6 -1025.453 28 0.3848 -2.2233 2.259
7 -1025.453 26 0.3848 -2.2233 2.259
8 -1025.453 33 0.3848 -2.2233 2.259
9 -1025.453 35 0.3848 -2.2233 2.259
10 -1025.453 46 0.3848 -2.2233 2.259
11 -1025.453 33 0.3848 -2.2233 2.259
12 -1025.453 41 0.3848 -2.2233 2.259 
multi_fit |>
  ggplot(aes(factor(start), final_loglik,
             fill = final_loglik)) +
  geom_col() +
  scale_fill_gradient(low = "#93C5FD", high = "#7C3AED") +
  labs(
    title = "Final log-likelihood from different EM starts",
    x = "Starting-value experiment",
    y = "Final log-likelihood",
    fill = "Log-likelihood"
  )

4.5 20. Real data 5: Gaussian mixture for faithful eruptions

The faithful dataset contains eruption durations and waiting times for the Old Faithful geyser. The eruption duration is strongly bimodal.

data(faithful)
faithful_fit <- em_gaussian2(
  faithful$eruptions,
  pi0 = 0.5,
  mu0 = c(2, 4.5),
  sigma0 = c(0.5, 0.5)
)

faithful_table <- tibble(
  Parameter = c("Mixing probability",
                "Short-eruption mean",
                "Long-eruption mean",
                "Short-eruption SD",
                "Long-eruption SD",
                "Iterations"),
  Estimate = c(
    faithful_fit$pi,
    faithful_fit$mu,
    faithful_fit$sigma,
    faithful_fit$iterations
  )
)

pretty_table(faithful_table,
             caption = "Two-component EM fit to Old Faithful eruptions")
Two-component EM fit to Old Faithful eruptions
Parameter Estimate
Mixing probability 0.3484
Short-eruption mean 2.0186
Long-eruption mean 4.2733
Short-eruption SD 0.2356
Long-eruption SD 0.4371
Iterations 31.0000 
faith_grid <- seq(min(faithful$eruptions) - 0.2,
                  max(faithful$eruptions) + 0.2,
                  length.out = 500)

faith_density <- tibble(
  x = faith_grid,
  Short =
    faithful_fit$pi *
    dnorm(faith_grid, faithful_fit$mu[1],
          faithful_fit$sigma[1]),
  Long =
    (1 - faithful_fit$pi) *
    dnorm(faith_grid, faithful_fit$mu[2],
          faithful_fit$sigma[2])
) |>
  mutate(Mixture = Short + Long)

ggplot(faithful, aes(eruptions)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 30, fill = "#DDE7F0",
                 color = "white") +
  geom_line(data = faith_density,
            aes(x, Short, color = "Short component"),
            linewidth = 1.2) +
  geom_line(data = faith_density,
            aes(x, Long, color = "Long component"),
            linewidth = 1.2) +
  geom_line(data = faith_density,
            aes(x, Mixture, color = "Fitted mixture"),
            linewidth = 1.5) +
  scale_color_manual(values = c(
    "Short component" = "#0EA5E9",
    "Long component" = "#F97316",
    "Fitted mixture" = "#111827"
  )) +
  labs(
    title = "EM decomposition of Old Faithful eruption duration",
    x = "Eruption duration (minutes)",
    y = "Density",
    color = NULL
  )

4.6 21. Real data 6: Soft clustering of faithful

faithful_classified <- faithful |>
  mutate(
    Probability_short = faithful_fit$responsibility,
    Assigned_component =
      if_else(Probability_short >= 0.5,
              "Short eruption", "Long eruption")
  )

ggplot(faithful_classified,
       aes(eruptions, waiting,
           color = Assigned_component,
           size = abs(Probability_short - 0.5))) +
  geom_point(alpha = 0.8) +
  scale_color_manual(values = c(
    "Short eruption" = "#0EA5E9",
    "Long eruption" = "#F97316"
  )) +
  scale_size_continuous(
    name = "Classification certainty",
    range = c(1.5, 4.5)
  ) +
  labs(
    title = "EM-based soft classification of Old Faithful eruptions",
    x = "Eruption duration",
    y = "Waiting time",
    color = "Assigned component"
  )

5 Part IV: Integrated Comparison

5.1 22. What question does each method answer?

method_comparison <- tibble(
  Method = c("Bootstrap", "Jackknife", "EM algorithm"),
  Main_question = c(
    "What is the sampling distribution or uncertainty of a statistic?",
    "How do deletion, bias, and individual observations affect an estimator?",
    "How can likelihood be maximized with missing or latent data?"
  ),
  Repeated_operation = c(
    "Resample observations with replacement",
    "Delete observations or groups",
    "Alternate conditional expectation and maximization"
  ),
  Principal_output = c(
    "SE, bias, CI, sampling distribution",
    "SE, bias correction, pseudovalues, influence",
    "MLEs, posterior responsibilities, completed sufficient statistics"
  ),
  Main_warning = c(
    "Resampling scheme must mimic data generation",
    "May fail for nonsmooth or dependent data",
    "May converge slowly or to a local maximum"
  )
)

pretty_table(method_comparison,
             caption = "Bootstrap, jackknife, and EM compared")
Bootstrap, jackknife, and EM compared
Method Main_question Repeated_operation Principal_output Main_warning
Bootstrap What is the sampling distribution or uncertainty of a statistic? Resample observations with replacement SE, bias, CI, sampling distribution Resampling scheme must mimic data generation
Jackknife How do deletion, bias, and individual observations affect an estimator? Delete observations or groups SE, bias correction, pseudovalues, influence May fail for nonsmooth or dependent data
EM algorithm How can likelihood be maximized with missing or latent data? Alternate conditional expectation and maximization MLEs, posterior responsibilities, completed sufficient statistics May converge slowly or to a local maximum

5.2 23. A combined simulation: uncertainty after EM

A useful advanced workflow is:

  1. fit a mixture model by EM;
  2. bootstrap the observations;
  3. refit the model in each bootstrap sample;
  4. use the resulting parameter distribution to estimate uncertainty.

Because mixture labels can switch, order the component means after every fit.

bootstrap_em <- function(x, B = 300, seed = 2026) {
  set.seed(seed)

  map_dfr(seq_len(B), function(b) {
    xb <- sample(x, replace = TRUE)

    fit <- try(
      em_gaussian2(
        xb,
        pi0 = 0.5,
        mu0 = quantile(xb, c(0.25, 0.75)),
        sigma0 = rep(sd(xb), 2)
      ),
      silent = TRUE
    )

    if (inherits(fit, "try-error")) {
      return(tibble(
        pi = NA_real_, mu_low = NA_real_,
        mu_high = NA_real_,
        sd_low = NA_real_, sd_high = NA_real_
      ))
    }

    ord <- order(fit$mu)

    tibble(
      pi = if (ord[1] == 1) fit$pi else 1 - fit$pi,
      mu_low = fit$mu[ord[1]],
      mu_high = fit$mu[ord[2]],
      sd_low = fit$sigma[ord[1]],
      sd_high = fit$sigma[ord[2]]
    )
  }) |>
    drop_na()
}

set.seed(401)
faithful_em_boot <- bootstrap_em(
  faithful$eruptions,
  B = 400
)

em_uncertainty <- faithful_em_boot |>
  pivot_longer(everything(),
               names_to = "Parameter",
               values_to = "Estimate") |>
  group_by(Parameter) |>
  summarise(
    Mean = mean(Estimate),
    SE = sd(Estimate),
    Lower = quantile(Estimate, 0.025),
    Upper = quantile(Estimate, 0.975),
    .groups = "drop"
  )

pretty_table(em_uncertainty,
             caption = "Bootstrap uncertainty for EM mixture parameters")
Bootstrap uncertainty for EM mixture parameters
Parameter Mean SE Lower Upper
mu_high 4.2742 0.0376 4.2017 4.3506
mu_low 2.0229 0.0312 1.9665 2.0889
pi 0.3516 0.0271 0.2966 0.4008
sd_high 0.4305 0.0333 0.3735 0.4964
sd_low 0.2384 0.0309 0.1907 0.3111 
faithful_em_boot |>
  select(mu_low, mu_high) |>
  pivot_longer(everything(),
               names_to = "Component",
               values_to = "Mean") |>
  mutate(Component = recode(
    Component,
    mu_low = "Short-eruption mean",
    mu_high = "Long-eruption mean"
  )) |>
  ggplot(aes(Mean, fill = Component)) +
  geom_density(alpha = 0.58) +
  facet_wrap(~ Component, scales = "free") +
  scale_fill_manual(values = c(
    "Short-eruption mean" = "#0EA5E9",
    "Long-eruption mean" = "#F97316"
  )) +
  labs(
    title = "Bootstrap uncertainty after EM fitting",
    subtitle = "Old Faithful mixture component means",
    x = "Bootstrap EM estimate",
    y = "Density"
  ) +
  guides(fill = "none")

6 Part V: Practical Assignments

6.1 24. Guided tasks

6.1.1 Task 1: Bootstrap a coefficient of variation

Generate 50 Gamma observations and estimate

\[ CV=\frac{S}{\bar X}. \]

Obtain the bootstrap bias, standard error, percentile interval, and basic interval.

6.1.2 Task 2: Bootstrap an odds ratio

Create a \(2\times2\) table from simulated binary exposure and outcome data. Bootstrap subjects, not table cells, and estimate the uncertainty of the odds ratio.

6.1.3 Task 3: Jackknife versus bootstrap for the median

Use at least three sample sizes and compare both standard errors. Explain why the methods disagree in small samples.

6.1.4 Task 4: Real-data influence analysis

For airquality, estimate the correlation between Ozone and Temp after removing missing values. Use jackknife influence values to identify the five most influential days.

6.1.5 Task 5: EM for a mixture under overlapping components

Simulate a two-normal mixture whose means are close together. Fit it from at least ten initial values. Discuss identifiability, slow convergence, and instability.

6.1.6 Task 6: EM plus bootstrap

Fit the Old Faithful mixture model, bootstrap it, and report standard errors and percentile intervals for both component means.

7 Part VI: Suggested Solutions

7.1 25. Solution skeletons

7.1.1 Solution to Task 1

set.seed(501)
x <- rgamma(50, shape = 4, rate = 2)

cv <- function(z) sd(z) / mean(z)
answer <- bootstrap_1d(x, cv, B = 5000)

answer$estimate
answer$bias
answer$se
answer$percentile_ci
answer$basic_ci

tibble(CV = answer$replicates) |>
  ggplot(aes(CV)) +
  geom_density(fill = "#5B8FF9", alpha = 0.6)

7.1.2 Solution to Task 2

set.seed(502)
n <- 400
exposure <- rbinom(n, 1, 0.45)
prob <- plogis(-1 + 1.1 * exposure)
outcome <- rbinom(n, 1, prob)

d <- tibble(exposure, outcome)

odds_ratio <- function(data) {
  tab <- table(data$exposure, data$outcome)
  (tab[2, 2] * tab[1, 1]) /
    (tab[2, 1] * tab[1, 2])
}

set.seed(503)
or_star <- replicate(5000, {
  idx <- sample(seq_len(nrow(d)), replace = TRUE)
  odds_ratio(d[idx, ])
})

quantile(or_star, c(0.025, 0.5, 0.975))
sd(or_star)

7.1.3 Solution to Task 3

compare_median <- function(n, reps = 250) {
  map_dfr(seq_len(reps), function(i) {
    x <- rexp(n)
    jk <- jackknife(x, median)
    bs <- bootstrap_1d(x, median, B = 1000,
                       seed = sample.int(1e7, 1))
    tibble(n = n, jackknife_se = jk$se,
           bootstrap_se = bs$se)
  })
}

result <- bind_rows(
  compare_median(15),
  compare_median(30),
  compare_median(80)
)

7.1.4 Solution to Task 4

aq <- airquality |>
  select(Ozone, Temp) |>
  drop_na()

aq_cor <- function(d) cor(d$Ozone, d$Temp)
jk_aq <- jackknife(aq, aq_cor)

top5 <- order(abs(jk_aq$influence),
              decreasing = TRUE)[1:5]

cbind(aq[top5, ],
      Influence = jk_aq$influence[top5])

7.1.5 Solution to Task 5

set.seed(505)
z <- rbinom(500, 1, 0.5)
x <- ifelse(z == 1,
            rnorm(500, 0, 1),
            rnorm(500, 1.1, 1))

# Repeat em_gaussian2() using many dispersed starts.
# Compare final log-likelihoods and estimates.

7.1.6 Solution to Task 6

Use the bootstrap_em() function from Section 23 and increase \(B\) to at least 1,000 for a final report.

8 Reproducibility information

sessionInfo()
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_India.utf8  LC_CTYPE=English_India.utf8   
## [3] LC_MONETARY=English_India.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=English_India.utf8    
## 
## time zone: Asia/Calcutta
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] boot_1.3-32   scales_1.4.0  knitr_1.51    tibble_3.3.1  purrr_1.2.2  
## [6] tidyr_1.3.2   dplyr_1.2.1   ggplot2_4.0.3
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3        cli_3.6.6          rlang_1.2.0        xfun_0.60         
##  [5] generics_0.1.4     S7_0.2.2           jsonlite_2.0.0     labeling_0.4.3    
##  [9] glue_1.8.1         htmltools_0.5.9    sass_0.4.10        rmarkdown_2.31    
## [13] grid_4.6.0         evaluate_1.0.5     jquerylib_0.1.4    fastmap_1.2.0     
## [17] yaml_2.3.12        lifecycle_1.0.5    compiler_4.6.0     RColorBrewer_1.1-3
## [21] pkgconfig_2.0.3    farver_2.1.2       digest_0.6.39      R6_2.6.1          
## [25] tidyselect_1.2.1   pillar_1.11.1      magrittr_2.0.5     bslib_0.12.0      
## [29] withr_3.0.2        tools_4.6.0        gtable_0.3.6       cachem_1.1.0

Final practical message

Bootstrap approximates repeated sampling. Jackknife approximates sensitivity to deletion. EM transforms latent-data likelihood maximization into a sequence of simpler conditional calculations. The methods solve different problems, but they can be combined—for example, bootstrap can quantify uncertainty after an EM fit.