Zaawansowane Metody Repróbkowania w Analizie Statystycznej

Przewodnik metodologiczny z implementacją w języku R

Author

Julia Wydra

Published

May 20, 2026

1 Wstęp do metod repróbkowania

Metody repróbkowania (*resampling*) stanowią fundament współczesnej statystyki obliczeniowej. Główną ideą tych technik jest wielokrotne wykorzystanie informacji zawartych w jednej, dostępnej próbie badawczej w celu oceny stabilności, błędów standardowych oraz rozkładów stosowanych estymatorów. Wybór odpowiedniej metody zależy od struktury danych oraz właściwości matematycznych samego estymatora.

**Ważna zasada metodologiczna:** Metody takie jak Bootstrap i Jackknife nie generują nowych informacji o populacji, a jedynie pozwalają w sposób efektywny wydobyć informacje o zmienności estymatora ukryte wewnątrz posiadanej próby losowej.

2 Głęboka Analiza: Metoda Bootstrap

Metoda bootstrap, wprowadzona przez Bradleya Efrona w 1979 roku, zrewolucjonizowała podejście do wnioskowania nieparametrycznego. Opiera się na tzw. **zasadzie substytucji** (*plug-in principle*). Jeśli nie znamy prawdziwego rozkładu populacji \[F\] zastępujemy go rozkładem empirycznym \[\hat{F}_n\] przypisując prawdopodobieństwo \[1/n\] każdej z obserwacji w próbie.

2.1 Algorytm i Teoretyczne Podstawy

  1. Dysponujemy próbą pierwotną \[{X} = (x_1, x_2, \dots, x_n)\]

  2. Generujemy $\(B\) prób bootstrapowych \[\mathbf{X}^* = (x_1^*, x_2^*, \dots, x_n^*)\] poprzez **losowanie ze zwracaniem** o identycznej liczności \[n\]

  3. Dla każdej próby obliczamy replikację interesującej nas statystyki \(\theta_b^* = \hat{\theta}(\mathbf{X}^*)\)

    Błąd standardowy bootstrapowy wyznaczany jest jako odchylenie standardowe z uzyskanych \[B\] replikacji:

    $$\widehat{SE}_B(\hat{\theta}) = \sqrt{\frac{1}{B-1}\sum_{b=1}^B \left(\theta_b^* - \bar{\theta}^*\right)^2}$$

    Gdzie \(\bar{\theta}^* = \frac{1}{B}\sum_{b=1}^B \theta_b^*\)

    2.2 Szacowanie Wariancji i Błędu Standardowego Estymatora

    Kluczowym celem bootstrapu jest wyznaczenie charakterystyk rozkładu estymatora \(\hat{\theta}\), kiedy jego analityczna postać wariancji jest nieznana. Empiryczna wariancja bootstrapowa estymatora wyraża się wzorem:

    $$\widehat{Var}_B(\hat{\theta}) = \frac{1}{B-1}\sum_{b=1}^B \left(\theta_b^* - \bar{\theta}^*\right)^2$$

    Gdzie błąd standardowy \(\widehat{SE}_B(\hat{\theta})\) jest pierwiastkiem kwadratowym z tej wariancji, a średnia z replikacji bootstrapowych wynosi:

    $$\bar{\theta}^* = \frac{1}{B}\sum_{b=1}^B \theta_b^*$$

    2.3 Szacowanie Obciążenia (Bias) i Korekta Estymatora

    Bootstrap pozwala również ocenić, jak bardzo nasz estymator z próby \(\hat{\theta}\) jest przesunięty względem parametru populacyjnego. Obciążenie bootstrapowe definiujemy jako:

    $$\widehat{Bias}_B(\hat{\theta}) = \bar{\theta}^* - \hat{\theta}$$

    Dzięki temu możemy wyznaczyć estymator skorygowany o obciążenie (\(\hat{\theta}_{Cal}\)):

    $$\hat{\theta}_{Cal} = \hat{\theta} - \widehat{Bias}_B(\hat{\theta}) = 2\hat{\theta} - \bar{\theta}^*$$

    2.4 Konstrukcja Przedziałów Ufności

    W zależności od zachowania rozkładu empirycznego wyróżniamy kilka metod budowy przedziałów ufności na poziomie istotności \(\alpha\):

    1. Przedział Percentylowy: Wykorzystuje kwantyle rozkładu empirycznego \(\theta^*\).

      $$[\theta^*_{(\alpha/2)}, \theta^*_{(1-\alpha/2)}]$$

    2. Przedział Podstawowy (Pivotowy): Bierze pod uwagę potencjalne obciążenie estymatora, odwracając relację kwantyli.

      $$[2\hat{\theta} - \theta^*_{(1-\alpha/2)}, 2\hat{\theta} - \theta^*_{(\alpha/2)}]$$

Wyjaśnienie inwersji kwantyli: W przedziale podstawowym (pivotowym) zauważyć można odwrócenie pozycji kwantyli — kwantyl rzędu (1−α/2) modyfikuje lewą (dolną) granicę przedziału. Wynika to stąd, że duże dodatnie błędy w symulacji bootstrapowej (θ∗≫θ^) oznaczają, iż estymator z próby ma tendencję do przeszacowania parametru, co wymaga silniejszej korekty w dół przy wyznaczaniu dolnej granicy rzeczywistego przedziału ufności.

### Implementacja i Symulacja w R

Proces symulacji komputerowej podzielimy na bloki kodu: generowanie danych wejściowych, wykonanie algorytmu ręcznego oraz weryfikację za pomocą pakietu systemowego wraz z wykresem.
set.seed(2026)

dane_analiza <- rexp(60, rate = 0.25)

summary(dane_analiza)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3168  1.0230  3.3406  4.8167  5.7794 25.8587 
B <- 2000
n <- length(dane_analiza)
boot_manualny <- numeric(B)

for(b in 1:B) {
  proba_b <- sample(dane_analiza, size = n, replace = TRUE)
  boot_manualny[b] <- sd(proba_b) / mean(proba_b)
}

se_boot_manual <- sd(boot_manualny)
ci_boot_manual <- quantile(boot_manualny, c(0.025, 0.975))

print(paste("Manualny błąd standardowy (SE):", round(se_boot_manual, 4)))
[1] "Manualny błąd standardowy (SE): 0.1072"
library(boot)

v_func <- function(data, indices) {
  d <- data[indices]
  return(sd(d) / mean(d))
}


boot_systemowy <- boot(data = dane_analiza, statistic = v_func, R = B) 
se_boot_sys <- sd(boot_systemowy$t)

print(paste("Systemowy błąd standardowy (SE):", round(se_boot_sys, 4)))
[1] "Systemowy błąd standardowy (SE): 0.1052"
hist(boot_manualny, breaks = 35, probability = TRUE, col = "lightblue", 
     main = "Rozkład Bootstrapowy Współczynnika Zmienności",
     xlab = "Wartości V*", border = "gray")
lines(density(boot_manualny), col = "darkblue", lwd = 3)

2.5 Bootstrap dla Indeksów Agregatowych

W analizach ekonomiczno-społecznych techniki bootstrapowe znajdują szerokie zastosowanie przy szacowaniu precyzji indeksów agregatowych (np. indeksów cenowych Laspeyres’a, Paaschego czy Fishera). Ponieważ wskaźniki te stanowią nieliniowe relacje wielu zmiennych, wyznaczenie ich dokładnego rozkładu teoretycznego jest niezwykle skomplikowane.

Jako przykład rozpatrzmy agregatowy indeks cen typu Laspeyres’a:

$$I_L = \frac{\sum_{i=1}^n p_{1i} q_{0i}}{\sum_{i=1}^n p_{0i} q_{0i}} \times 100$$

Gdzie \(p_{1i}\) to cena w okresie badanym, \(p_{0i}\) cena w okresie bazowym, a \(q_{0i}\) to stała wagowa (ilość) z okresu bazowego. Próbkowanie bootstrapowe polega tutaj na losowaniu ze zwracaniem całych wektorów informacji o produktach.

set.seed(123)
n_produkty <- 50
p0 <- runif(n_produkty, 10, 100)
p1 <- p0 * rlnorm(n_produkty, meanlog = 0.04, sdlog = 0.08)
q0 <- round(runif(n_produkty, 5, 150))

calc_laspeyres <- function(p1_v, p0_v, q0_v) {
  sum(p1_v * q0_v) / sum(p0_v * q0_v) * 100
}

B_indeks <- 2000
boot_indeks <- numeric(B_indeks)

for(b in 1:B_indeks) {
  idx_b <- sample(1:n_produkty, size = n_produkty, replace = TRUE)
  boot_indeks[b] <- calc_laspeyres(p1[idx_b], p0[idx_b], q0[idx_b])
}

se_laspeyres <- sd(boot_indeks)
ci_laspeyres <- quantile(boot_indeks, c(0.025, 0.975))

cat("Wartość indeksu z próby głównej:", calc_laspeyres(p1, p0, q0), "\n")
Wartość indeksu z próby głównej: 105.3176 
cat("Błąd standardowy bootstrap:", se_laspeyres, "\n")
Błąd standardowy bootstrap: 1.250947

3 Głęboka Analiza: Metoda Jackknife

Metoda Jackknife (Quenouille, 1949; Tukey, 1958) to technika deterministyczna – nie występuje w niej stochastyczne losowanie. Jest historyczną poprzedniczką bootstrapu, zaprojektowaną głównie do redukcji obciążenia (bias reduction).

3.1 Intuicja i zastosowanie

W odróżnieniu od bootstrapu, Jackknife nie losuje niczego na ślepo. Działa metodycznie: usuwa pierwszą obserwację i przelicza wynik, potem ją przywraca, usuwa drugą i znowu przelicza wynik – i tak do końca próby. Pozwala to precyzyjnie ocenić stabilność modelu i sprawdzić, czy w próbie nie ma jednej „wystającej” obserwacji, która całkowicie fałszuje ostateczne wnioski.

Główną wadą metody jest to, że zupełnie nie radzi sobie ze statystykami „niegładkimi”, takimi jak mediana.

::: {.styled-note} Uwaga na ograniczenia: Ponieważ Jackknife opiera się na usuwaniu pojedynczych jednostek, całkowicie zawodzi w przypadku szacowania błędów dla kwantyli i mediany, ponieważ drobna zmiana jednej wartości nie zmienia wartości środkowej próby w sposób ciągły. :::

3.2 Pseudowartości Tukeya i Estymacja Wariancji

Algorytm opiera się na sekwencyjnym usuwaniu z próby dokładnie jednej obserwacji (metoda Leave-One-Out). Aby zredukować obciążenie estymatorów nieliniowych, Tukey wprowadził transformację w tzw. pseudowartości (\(\tilde{\theta}_i\)):

$$\tilde{\theta}_i = n\hat{\theta} - (n-1)\hat{\theta}_{(-i)}$$

Gdzie \(\hat{\theta}\) to ocena z pełnej próby, a \(\hat{\theta}_{(-i)}\) to ocena obliczona po wycięciu \(i\)-tej obserwacji. Ogólny błąd standardowy Jackknife dany jest wzorem:

$$\widehat{SE}_{JK}(\hat{\theta}) = \sqrt{\frac{1}{n(n-1)}\sum_{i=1}^n \left(\tilde{\theta}_i - \bar{\tilde{\theta}}\right)^2}$$

3.3 Implementacja i Symulacja w R

Zastosujemy obciążony estymator wariancji (dzielony przez \(n\)), aby sprawdzić sprawność algorytmu Jackknife w usuwaniu obciążenia.

obciazona_var <- function(x) { sum((x - mean(x))^2) / length(x) }
theta_full <- obciazona_var(dane_analiza)

n_jk <- length(dane_analiza)
jack_wyszukane <- numeric(n_jk)


for(i in 1:n_jk) {
  jack_wyszukane[i] <- obciazona_var(dane_analiza[-i])
}

pseudowartosci <- n_jk * theta_full - (n_jk - 1) * jack_wyszukane
se_jack_manual <- sqrt(var(pseudowartosci) / n_jk)

print(paste("Jackknife błąd standardowy (SE):", round(se_jack_manual, 4)))
[1] "Jackknife błąd standardowy (SE): 9.8073"
plot(jack_wyszukane, type = "h", col = "#1e40af", lwd = 2,
     main = "Wpływ usunięcia i-tej obserwacji na estymator Jackknife",
     xlab = "Indeks usuniętej obserwacji", ylab = "Wartość estymatora cząstkowego")
abline(h = theta_full, col = "#eab308", lty = 2, lwd = 2)

3.3.1 Wykres diagnostyczny wpływu obserwacji jednostkowych

Poniższy kod generuje wykres szpilkowy, który pozwala sprawdzić, jak usunięcie konkretnego elementu z bazy danych wpływa na wariancję, ułatwiając tym samym identyfikację obserwacji odstających.

4 Porównanie Właściwości: Bootstrap vs Jackknife

Poniższa tabela przedstawia zestawienie krytycznych różnic operacyjnych między obiema metodami

             Kryterium                Bootstrap                   Jackknife
1      Charakterystyka   Losowy (ze zwracaniem)      Deterministyczny (LOO)
2        Liczba kroków      Dowolna (B >= 1000)          Sztywna (zawsze n)
3           Główny cel     Błędy SE, przedziały     Redukcja obciążenia, SE
4 Podatność na medianę Niska (działa poprawnie) Wysoka (całkowicie zawodzi)

5 Metody Grup Losowych oraz BHS

W badaniach masowych o złożonych schematach losowania stosuje się podejścia strukturalne.

5.1 Metoda Niezależnych/Zależnych Grup Losowych - Opis

Metody te polegają na kontrolowanym dzieleniu wielkich zbiorów danych (np. danych ze spisów powszechnych GUS) na mniejsze, autonomiczne pakiety (podgrupy). Zamiast liczyć błąd dla wielomilionowej próby na raz, liczy się statystyki wewnątrz mniejszych grup, a następnie bada się rozbieżności pomiędzy nimi. Zapobiega to przeciążeniom pamięci komputera i uwzględnia strukturę warstwową badania.

5.2 Metoda Niezależnych Grup Losowych

Stosowana, gdy proces zbierania danych opiera się na \(R\) całkowicie niezależnych powtórzeniach identycznego schematu losowania:

$$\widehat{Var}(\bar{\theta}) = \frac{1}{R(R-1)}\sum_{r=1}^R (\hat{\theta}_r - \bar{\theta})^2$$

5.3 Metoda Zależnych Grup Losowych (Podział Post-Hoc)

Podział istniejącej próby na \(R\) podgrup za pomocą operatora modulo.

R_grupy <- 5
podzial_indeksow <- (1:length(dane_analiza)) %% R_grupy
list_grup <- split(dane_analiza, podzial_indeksow)

sapply(list_grup, mean)
       0        1        2        3        4 
2.359894 6.123962 4.064196 4.666025 6.869247

5.4 Metoda Zrównoważonych Półpróbek Powtarzanych (BHS) - Opis

Metoda BHS (Balanced Repeated Replication) to wysoce wyspecjalizowana technika używana w badaniach socjoekonomicznych, w których populacja jest podzielona na pary jednostek (np. dwa miasta z każdego województwa). Z każdej pary wybiera się losowo po jednym mieście, tworząc tzw. półpróbkę. Aby zachować matematyczną równowagę i ortogonalność wyborów, wykorzystuje się specjalne macierze kombinatoryczne – macierze Hadamarda.

5.5 Metoda Zrównoważonych Półpróbek Powtarzanych (BHS)

Wykorzystuje ortogonalne Macierze Hadamarda do szacowania wariancji w schematach warstwowych.

macierz_hadamarda <- matrix(c(1,  1,  1,  1, 
                              1, -1,  1, -1, 
                              1,  1, -1, -1, 
                              1, -1, -1,  1), nrow = 4, byrow = TRUE)
print(macierz_hadamarda)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1   -1    1   -1
[3,]    1    1   -1   -1
[4,]    1   -1   -1    1

6 Bibliografia i Literatura Źródłowa

  1. Efron, B. (1979) – Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, Vol. 7, No. 1, str. 2–8. Link do pełnego tekstu PDF

  2. Wasserman, L. (2004) – All of Statistics: A Concise Course in Statistical Inference. Carnegie Mellon University, Springer, str. 97-116. Link do darmowego podręcznika PDF

  3. Wolter, K. M. (2007) – Introduction to Variance Estimation. US Census Bureau, Springer, Second Edition, str. 22-28; 32-50; 151-169; 194-214 Pobierz podręcznik Woltera

  4. Lednicki, B., Wieczorkowski, R. (2020). Metody komputerowe szacowania precyzji zmiennych w badaniach reprezentacyjnych. Wiadomości Statystyczne. The Polish Statistician, nr 6, str. 23–31.

  5. Wywiał, J. L. (2014). Contributions to Testing Statistical Hypotheses in Auditing. Wydawnictwo Uniwersytetu Ekonomicznego w Katowicach.

    • Sekcja 3.3.3 (str. 34): Teoretyczne podstawy podejścia bootstrapowego (Bootstrap approach) w szacowaniu parametrów