Bootstrap i jackknife

Author

Natalia Bąchor

Wstęp

Projekt dotyczy metod resamplingu (ponownego próbkowania), które stanowią narzędzie współczesnej statystyki obliczeniowej. Skupia się na dwóch klasycznych podejściach nieparametrycznych: metodzie Bootstrap oraz metodzie Jackknife.

Głównym celem stosowania tych algorytmów w środowisku R jest estymacja właściwości rozkładu próbkowego danych statystyk w sytuacjach, gdy tradycyjne podejście asymptotyczne zawodzi z powodu małej liczebności próby lub braku znajomości teoretycznego rozkładu populacji.

Rozwiń kod źródłowy 
library(ggplot2)
library(knitr)
set.seed(616)

1. Metoda Jackknife

Metoda Jackknife polega na sekwencyjnym usuwaniu dokładnie jednej obserwacji z oryginalnej próby o wielkości \(n\) i ponownym obliczaniu statystyki dla podprób o wielkości \(n-1\).

Przykład estymatora odchylenia standardowego w małej próbie

Rozwiń kod źródłowy 
n <- 15
dane_oryginalne <- rexp(n, rate = 0.5)

theta_hat <- sd(dane_oryginalne)

jack_estimates <- numeric(n)
for(i in 1:n) {
  jack_estimates[i] <- sd(dane_oryginalne[-i])
}

theta_jack_mean <- mean(jack_estimates)

bias_jack <- (n - 1) * (theta_jack_mean - theta_hat)

se_jack <- sqrt(((n - 1) / n) * sum((jack_estimates - theta_jack_mean)^2))

wyniki_jack <- data.frame(
  Metryka = c("Statystyka z próby", "Estymowane obciążenie Jackknife", "Estymowany błąd standardowy Jackknife"),
  Wartosc = c(theta_hat, bias_jack, se_jack)
)
kable(wyniki_jack, caption = "Wyniki analizy obciążenia i błędu metodą Jackknife")
Wyniki analizy obciążenia i błędu metodą Jackknife
Metryka Wartosc
Statystyka z próby 1.5638318
Estymowane obciążenie Jackknife -0.0375341
Estymowany błąd standardowy Jackknife 0.3424812

2.1 Wizualizacja replikacji Jackknife

Jackknife generuje dokładnie \(n\) wartości (po jednej dla każdej odrzuconej obserwacji), wykres prezentuje rozkład tych wartości i pokazuje wpływ poszczególnych punktów z próby na ostateczny wynik.

Rozwiń kod źródłowy 
df_jack <- data.frame(ID = 1:n, Wartosc = jack_estimates)

ggplot(df_jack, aes(x = ID, y = Wartosc)) +
  geom_bar(stat = "identity", fill = "#e67e22", width = 0.6) +
  geom_hline(yintercept = theta_hat, color = "darkred", linetype = "dashed", linewidth = 1) +
  annotate("text", x = 2, y = theta_hat * 1.02, label = "Wartość oryginalna (theta_hat)", color = "darkred", hjust = 0) +
  labs(title = "Wpływ usuwania pojedynczych obserwacji (Replikacje Jackknife)",
       x = "Indeks usuniętej obserwacji (i)", y = "Wartość estymatora dla próby bez i-tej obserwacji") +
  theme_minimal()

2. Metoda Bootstrap

Metoda Bootstrap, stanowi bezpośrednie rozwinięcie idei Jackknife. Bootstrap opiera się na losowaniu ze zwracaniem, próby o dokładnie takiej samej liczebności \(n\) jak próba oryginalna. Pozwala to na wyznaczenie pełnego empirycznego rozkładu statystyki i budowę m.in. percentylowych przedziałów ufności.

2.1 Symulacja i rozkład empiryczny Bootstrap

Generujemy rozkład empiryczny statystyki poprzez wykonanie \(B = 10000\) powtórzeń procedury losowania ze zwracaniem.

Rozwiń kod źródłowy 
B <- 10000     
boot_estimates <- numeric(B)

for (i in 1:B) {
  proba_boot <- sample(dane_oryginalne, size = n, replace = TRUE)
  boot_estimates[i] <- sd(proba_boot)
}

bias_boot <- mean(boot_estimates) - theta_hat
se_boot <- sd(boot_estimates)

ci_lower <- quantile(boot_estimates, 0.025)
ci_upper <- quantile(boot_estimates, 0.975)

cat("Estymowane obciążenie Bootstrap wynosi:", bias_boot, "\n")
Estymowane obciążenie Bootstrap wynosi: -0.07859855
Rozwiń kod źródłowy 
cat("Estymowany błąd standardowy Bootstrap wynosi:", se_boot, "\n")
Estymowany błąd standardowy Bootstrap wynosi: 0.2827635

2.2 Wizualizacja rozkładu Bootstrap

Poniższy wykres przedstawia gęstość empirycznego rozkładu wskaźnika, wraz z zaznaczoną wartością oryginalną \(\hat{\theta}\) oraz granicami 95% przedziału ufności.

Rozwiń kod źródłowy 
df_boot <- data.frame(Statystyka = boot_estimates)

ggplot(df_boot, aes(x = Statystyka)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "#3498db", color = "white", alpha = 0.8) +
  geom_density(color = "#2980b9", linewidth = 1) +
  geom_vline(xintercept = theta_hat, color = "darkred", linetype = "dashed", linewidth = 1.2) +
  geom_vline(xintercept = c(ci_lower, ci_upper), color = "#2ecc71", linetype = "dotted", linewidth = 1.2) +
  annotate("text", x = theta_hat + 0.1, y = 1.0, 
           label = paste("Obserwowana\nstatystyka T0 =", round(theta_hat, 3)), 
           color = "darkred", hjust = 0) +
  labs(title = "Rozkład empiryczny statystyki z symulacji Bootstrap",
       x = "Wartość odchylenia standardowego z replikacji", y = "Gęstość") +
  theme_minimal()

3. Porównanie efektywności metod: Hipotezy i stabilność

Podczas gdy Jackknife jest metodą w pełni powtarzalną (zawsze daje ten sam wynik dla danej próby), Bootstrap oferuje większą elastyczność i dokładność w przypadku statystyk nieliniowych lub niedoróżniczkowalnych (np. mediana).

Załóżmy, że chcemy porównać stabilność estymacji błędu standardowego dla mediany przy użyciu obu podejść na danych z asymetrycznego rozkładu.

Rozwiń kod źródłowy 
theta_hat_med <- median(dane_oryginalne)

# Jackknife dla mediany
jack_med <- numeric(n)
for(i in 1:n) { jack_med[i] <- median(dane_oryginalne[-i]) }
se_jack_med <- sqrt(((n - 1) / n) * sum((jack_med - mean(jack_med))^2))

# Bootstrap dla mediany
boot_med <- numeric(B)
for(i in 1:B) { boot_med[i] <- median(sample(dane_oryginalne, size = n, replace = TRUE)) }
se_boot_med <- sd(boot_med)

wyniki_porownania <- data.frame(
  Metoda = c("Jackknife (deterministyczna)", "Bootstrap (symulacyjna B=10000)"),
  Estymowany_Blad_Standardowy_Mediany = c(se_jack_med, se_boot_med)
)
kable(wyniki_porownania, caption = "Porównanie oszacowania błędu standardowego mediany")
Porównanie oszacowania błędu standardowego mediany
Metoda Estymowany_Blad_Standardowy_Mediany
Jackknife (deterministyczna) 1.0810284
Bootstrap (symulacyjna B=10000) 0.7180927

Wizualizacja porównawcza rozkładów replikacji

Wizualizacja róznicy między metodą Jackknife (zależna tylko od liczby obserwacji \(n\)), a metodą Bootstrap (rozkładu ciągłego).

Rozwiń kod źródłowy 
df_jack_plot <- data.frame(Wartosc = jack_med, Metoda = "Jackknife")
df_boot_plot <- data.frame(Wartosc = boot_med, Metoda = "Bootstrap")
df_combined <- rbind(df_jack_plot, df_boot_plot[1:500, ]) 

ggplot(df_combined, aes(x = Wartosc, fill = Metoda)) +
  geom_histogram(position = "dodge", bins = 20, color = "white") +
  scale_fill_manual(values = c("Jackknife" = "#e74c3c", "Bootstrap" = "#9b59b6")) +
  labs(title = "Porównanie charakteru replikacji metod Jackknife i Bootstrap dla mediany",
       x = "Wartość mediany w podpróbach", y = "Liczność (skalowana)") +
  theme_minimal()

Podsumowanie

Prezentowane techniki pokazują odejście od restrykcyjnych założeń parametrycznych na rzecz czystej mocy obliczeniowej. Metoda Jackknife świetnie sprawdza się do szybkiej i powtarzalnej redukcji obciążenia gładkich statystyk w bardzo małych próbach. Z kolei nieparametryczny Bootstrap, dzięki symulacji Monte Carlo, pozwala na pełne odzwierciedlenie niepewności estymacji i precyzyjne wyznaczanie przedziałów ufności bez względu na stopień skomplikowania formy matematycznej analizowanej statystyki.

Źródła