Test pustych cel

Author

Julia Szymczyk

Published

June 9, 2026

1. Wstęp i idea testu

Test pustych cel (znany też jako test Davida-Hellwiga) to nieparametryczny test statystyczny, służący do weryfikacji hipotezy, czy badana próba losowa pochodzi z określonego rozkładu (zazwyczaj bada się w ten sposób rozkład jednostajny).

W praktyce

Zmieniamy naszą próbę tak, aby (teoretycznie) pochodziła z rozkładu jednostajnego na przedziale \([0, 1]\). Następnie ten przedział dzielimy na \(m\) równych i rozłącznych “cel” (przedziałów). Wrzucamy do nich nasze \(n\) obserwacji z próby i zliczamy, ile cel pozostało całkowicie pustych. Statystyką testową jest \(K\) – liczba pustych cel.

Hipoteza zerowa (\(H_0\)): Próba pochodzi z testowanego rozkładu (obserwacje wpadają do cel losowo, rozpraszając się w miarę równomiernie).
Hipoteza alternatywna (\(H_1\)): Próba nie pochodzi z testowanego rozkładu (np. obserwacje zbijają się w klastry, zostawiając nienaturalnie dużo pustych cel).

2. Rozkład statystyki, parametry i liczba cel

Ile cel (\(m\)) należy przyjąć?

Z reguły w klasycznym podejściu przyjmuje się, że liczba cel jest równa liczbie obserwacji, czyli \(m = n\).

Rozkład statystyki \(K\)

Przy założeniu prawdziwości \(H_0\), dokładne prawdopodobieństwo tego, że liczba pustych cel wyniesie \(k\), wyraża się skomplikowanym wzorem kombinatorycznym:

\[P(K=k) = \binom{m}{k} \sum_{r=0}^{m-k} (-1)^r \binom{m-k}{r} \left(1 - \frac{k+r}{m}\right)^n\]

Dla testowania hipotez łatwiej jest oprzeć się na wartości oczekiwanej i wariancji tego rozkładu:

Wartość oczekiwana: \(E(K) = m \left(1 - \frac{1}{m}\right)^n\)
Wariancja: \(Var(K) = m(m-1)\left(1 - \frac{2}{m}\right)^n + m\left(1 - \frac{1}{m}\right)^n - m^2\left(1 - \frac{1}{m}\right)^{2n}\)

Przybliżenie dla dużych prób

Ponieważ liczenie silni z powyższego wzoru dla dużych prób “zapycha” pamięć komputerów, stosuje się przybliżenia:

Dla bardzo dużych \(n\) i \(m\), rozkład zbiega do rozkładu Poissona.
Asymptotycznie (gdy \(n\) i \(m\) rosną w stałej proporcji), rozkład statystyki \(K\) dąży do rozkładu normalnego.

3. P-value i tabele wartości krytycznych

Wartości krytyczne

W teście tym hipotezę zerową odrzucamy najczęściej wtedy, gdy pustych cel jest zbyt dużo (wartość \(K \ge k_{kryt}\)). Oznacza to, że dane nie są losowe, lecz skupiają się w jednym miejscu.Wartości krytyczne dla małych prób (np. do \(n = 30\)) odczytuje się ze specjalnych tablic statystycznych opracowanych przez Z. Hellwiga.

Wyznaczanie p-value

Dla większych prób (gdzie działa aproksymacja normalna), p-value wyznacza się standaryzując wynik. Ze względu na to, że statystyka \(K\) jest dyskretna (przyjmuje tylko wartości całkowite), stosuje się statystykę \(Z\) z tzw. poprawką na ciągłość:

\[Z = \frac{K - 0.5 - E(K)}{\sqrt{Var(K)}}\]

P-value odczytujemy następnie z tablic (lub funkcji) standardowego rozkładu normalnego jako prawdopodobieństwo prawostronne (\(P(Z > z_{obs})\)).

4. Prompty do generowania w języku R

Poniższe prompty zostały wysłane do Gemini, a następnie opracowane w celu wyrzucenia nieistotnych fragmentów.

Porównanie wyników testu pustych cel dla dwóch różnych zbiorów danych

Porównanie próby losowej (rozkład jednostajny) z próbą zgrupowaną (klastry)

“Napisz kod w R, który przeprowadza test pustych cel. Napisz funkcję liczącą ten test i wyznaczającą p-value. Następnie wygeneruj dwie sztuczne próby danych (jedną losową, a drugą zgrupowaną), podziel je na cele i zademonstruj działanie swojej funkcji”

# 1. Funkcja pomocnicza: Obliczanie dokładnego prawdopodobieństwa P(S=k)
# Oparta na zasadzie włączeń i wyłączeń.
prob_empty_cells <- function(k, n, m) {
  if (k < 0 || k > m - 1) return(0)
  
  sum_val <- sum(sapply(0:(m - k), function(v) {
    (-1)^v * choose(m - k, v) * (1 - (k + v) / m)^n
  }))
  
  return(choose(m, k) * sum_val)
}

# 2. Główna funkcja: Test pustych cel
# data    - wektor danych liczbowych
# m       - liczba cel (przedziałów)
# min_val - dolna granica zakresu (domyślnie minimum z próby)
# max_val - górna granica zakresu (domyślnie maksimum z próby)
empty_cells_test <- function(data, m, min_val = min(data), max_val = max(data)) {
  
  n <- length(data)
  
  # Tworzenie m przedziałów o równej szerokości
  breaks <- seq(min_val, max_val, length.out = m + 1)
  
  # Podział danych na cele
  podzial <- cut(data, breaks = breaks, include.lowest = TRUE)
  czestosci <- table(podzial)
  
  # Zliczenie zaobserwowanej liczby pustych cel
  S_obs <- sum(czestosci == 0)
  
  # Obliczenie p-value (prawdopodobieństwo, że liczba pustych cel >= S_obs)
  p_value <- sum(sapply(S_obs:(m - 1), prob_empty_cells, n = n, m = m))
  
  # Zwracamy listę z wynikami
  return(list(
    Metoda = "Test pustych cel",
    n = n,
    m = m,
    S_obs = S_obs,
    p_value = p_value,
    Czestosci_cel = czestosci
  ))
}

# DEMONSTRACJA DZIAŁANIA

# Ustawienie ziarna dla powtarzalności wyników
set.seed(42)

# Parametry globalne dla demonstracji
n_elementow <- 100
m_cel <- 20
zakres_min <- 0
zakres_max <- 100

# Próba 1: Rozkład całkowicie losowy (jednostajny U(0, 100))
# Elementy powinny być w miarę równomiernie rozłożone we wszystkich celach.
proba_losowa <- runif(n_elementow, min = zakres_min, max = zakres_max)

# Próba 2: Rozkład zgrupowany (klastry)
# Tworzymy dwa silne skupiska danych wokół wartości 20 i 80.
proba_zgrupowana <- c(
  rnorm(n_elementow / 2, mean = 20, sd = 4),
  rnorm(n_elementow / 2, mean = 80, sd = 4)
)
# Ograniczenie wygenerowanych wartości do zakresu [0, 100]
proba_zgrupowana <- pmin(pmax(proba_zgrupowana, zakres_min), zakres_max)


# PRZEPROWADZENIE TESTÓW

cat("TEST 1: PRÓBA LOSOWA (Rozkład Jednostajny)\n")

TEST 1: PRÓBA LOSOWA (Rozkład Jednostajny)

wynik_losowa <- empty_cells_test(proba_losowa, m = m_cel, min_val = zakres_min, max_val = zakres_max)

cat(sprintf("Liczba obserwacji (n): %d\n", wynik_losowa$n))

Liczba obserwacji (n): 100

cat(sprintf("Liczba cel (m): %d\n", wynik_losowa$m))

Liczba cel (m): 20

cat(sprintf("Liczba pustych cel (S_obs): %d\n", wynik_losowa$S_obs))

Liczba pustych cel (S_obs): 0

cat(sprintf("Wartość p-value: %.5f\n", wynik_losowa$p_value))

Wartość p-value: 1.00000

if(wynik_losowa$p_value < 0.05) {
  cat("Wniosek: ODRZUCAMY hipotezę zerową (wykryto klasteryzację).\n")
} else {
  cat("Wniosek: BRAK PODSTAW do odrzucenia hipotezy zerowej (dane wyglądają na losowe).\n")
}

Wniosek: BRAK PODSTAW do odrzucenia hipotezy zerowej (dane wyglądają na losowe).

cat("TEST 2: PRÓBA ZGRUPOWANA (Klastry)\n")

TEST 2: PRÓBA ZGRUPOWANA (Klastry)

wynik_zgrupowana <- empty_cells_test(proba_zgrupowana, m = m_cel, min_val = zakres_min, max_val = zakres_max)

cat(sprintf("Liczba obserwacji (n): %d\n", wynik_zgrupowana$n))

Liczba obserwacji (n): 100

cat(sprintf("Liczba cel (m): %d\n", wynik_zgrupowana$m))

Liczba cel (m): 20

cat(sprintf("Liczba pustych cel (S_obs): %d\n", wynik_zgrupowana$S_obs))

Liczba pustych cel (S_obs): 10

cat(sprintf("Wartość p-value: %.10f\n", wynik_zgrupowana$p_value))

Wartość p-value: 0.0000000000

if(wynik_zgrupowana$p_value < 0.05) {
  cat("Wniosek: ODRZUCAMY hipotezę zerową (wykryto znaczną klasteryzację!).\n")
} else {
  cat("Wniosek: BRAK PODSTAW do odrzucenia hipotezy zerowej.\n")
}

Wniosek: ODRZUCAMY hipotezę zerową (wykryto znaczną klasteryzację!).

Wniosek: * Próba losowa (p-value = 1.0, 0 pustych cel):Dane są idealnie rozproszone w przestrzeni. Nie ma w nich żadnych ukrytych skupisk, a ich układ jest w pełni równomierny i przypadkowy.

Próba zgrupowana (p-value = 0.0, 10 pustych cel): Dane silnie zbijają się w grupy (klastry). Fakt, że aż połowa przedziałów pozostała pusta, jednoznacznie wyklucza przypadek.

Teoretyczny rozkład i symulacja testu pustych cel

“Napisz kod w R, który pokaże jak wygląda rozkład testu pustych cel na wybranym przykładzie, wypisze wartość p-value i stworzy tabele wartości krytycznych”

# 1. Parametry naszego przykładu
n <- 60 # liczba elementów (np. rzutów / losowań)
m <- 20 # liczba kategorii (cel)

# 2. Funkcja obliczająca dokładne prawdopodobieństwo P(S=k)
prob_empty_cells <- function(k, n, m) {
  # Zabezpieczenie przed niemożliwymi wartościami
  if (k < 0 || k > m - 1) return(0)
  
  # Implementacja wzoru opartego na zasadzie włączeń i wyłączeń
  sum_val <- sum(sapply(0:(m - k), function(v) {
    (-1)^v * choose(m - k, v) * (1 - (k + v) / m)^n
  }))
  
  return(choose(m, k) * sum_val)
}

# Generujemy wektor wszystkich możliwych wartości pustych cel (od 0 do m-1)
k_vals <- 0:(m - 1)

# Obliczamy prawdopodobieństwa dla każdego k
probs <- sapply(k_vals, prob_empty_cells, n = n, m = m)

# 3. Symulacja danych i statystyka testowa
set.seed(42) # Ustawienie ziarna dla powtarzalności wyników
# Losujemy n elementów do m cel (rozkład jednostajny)
obserwacje <- sample(1:m, n, replace = TRUE)

# Zliczamy ile unikalnych cel zostało zajętych i ile pozostało pustych
zajete_cele <- length(unique(obserwacje))
S_obs <- m - zajete_cele

# 4. Obliczenie p-value
# Standardowo w teście pustych cel testuje się nadmierne skupienie danych,
# co objawia się WIĘKSZĄ niż oczekiwana liczbą pustych cel (prawostronny obszar krytyczny).
p_value <- sum(probs[k_vals >= S_obs])

# 5. Wizualizacja
# Tworzymy wykres słupkowy i zapisujemy współrzędne słupków (bp) do narysowania linii
bp <- barplot(probs, names.arg = k_vals, col = "steelblue", 
              main = sprintf("Rozkład liczby pustych cel (n=%d, m=%d)", n, m),
              xlab = "Liczba pustych cel (k)", ylab = "Prawdopodobieństwo P(S=k)",
              border = "white", ylim = c(0, max(probs) * 1.2))

# Zaznaczenie obserwowanej wartości na wykresie
idx <- which(k_vals == S_obs)
if(length(idx) > 0) {
  abline(v = bp[idx], col = "red", lwd = 2, lty = 2)
  legend("topright", legend = paste("S_obs =", S_obs), col = "red", lty = 2, lwd = 2)
}

# 6. Tabele wyników i wartości krytycznych
# Prawdopodobieństwo skumulowane prawego ogona: P(S >= k)
p_prawy_ogon <- rev(cumsum(rev(probs)))

tabela_rozkladu <- data.frame(
  k = k_vals,
  P_S_eq_k = round(probs, 6),   # Prawdopodobieństwo punktowe
  P_S_ge_k = round(p_prawy_ogon, 6) # P-value dla danej wartości (prawy ogon)
)

# Szukamy wartości krytycznych dla standardowych poziomów istotności
alpha_levels <- c(0.10, 0.05, 0.01, 0.001)
wartosci_krytyczne <- sapply(alpha_levels, function(a) {
  # Znajdź najmniejsze k, dla którego prawdopodobieństwo znalezienia się w ogonie <= alpha
  k_crit <- k_vals[p_prawy_ogon <= a]
  if(length(k_crit) > 0) min(k_crit) else NA
})

tabela_krytyczna <- data.frame(
  Poziom_Istotnosci = alpha_levels,
  Wartosc_Krytyczna_k = wartosci_krytyczne
)

# 7. Wypisanie wyników do konsoli
cat("=== WYNIKI TESTU PUSTYCH CEL ===\n")

=== WYNIKI TESTU PUSTYCH CEL ===

cat(sprintf("Liczba kul (n): %d\n", n))

Liczba kul (n): 60

cat(sprintf("Liczba cel (m): %d\n", m))

Liczba cel (m): 20

cat(sprintf("Zaobserwowana liczba pustych cel (S_obs): %d\n", S_obs))

Zaobserwowana liczba pustych cel (S_obs): 2

cat(sprintf("P-value P(S >= %d): %.5f\n\n", S_obs, p_value))

P-value P(S >= 2): 0.22906

cat("=== TABELA WARTOŚCI KRYTYCZNYCH (Prawy ogon) ===\n")

=== TABELA WARTOŚCI KRYTYCZNYCH (Prawy ogon) ===

cat("Odrzucamy hipotezę zerową o losowości (jednostajności) na danym poziomie\n")

Odrzucamy hipotezę zerową o losowości (jednostajności) na danym poziomie

cat("istotności, jeżeli zaobserwowane S >= Wartość Krytyczna.\n")

istotności, jeżeli zaobserwowane S >= Wartość Krytyczna.

print(tabela_krytyczna)

  Poziom_Istotnosci Wartosc_Krytyczna_k
1             0.100                   3
2             0.050                   3
3             0.010                   4
4             0.001                   5

cat("\n=== FRAGMENT ROZKŁADU ===\n")


=== FRAGMENT ROZKŁADU ===

# Wypisujemy tylko wartości mające statystyczne znaczenie (żeby nie zaśmiecać konsoli zerami)
print(tabela_rozkladu[tabela_rozkladu$P_S_eq_k > 0.0001, ], row.names = FALSE)

 k P_S_eq_k P_S_ge_k
 0 0.360605 1.000000
 1 0.410331 0.639395
 2 0.182182 0.229064
 3 0.041247 0.046882
 4 0.005232 0.005635
 5 0.000386 0.000403

Wniosek: Ponieważ zasymulowane w kodzie dane pochodzą z rozkładu jednostajnego (sample z opcją replace = TRUE), zaobserwowana liczba pustych cel wpisuje się w środek wyliczonego rozkładu teoretycznego. Wynikowe p-value jest wysoki, jednoznacznie pokazując, że w danych nie ma ukrytych skupisk, a ich układ jest w pełni równomierny i przypadkowy.

Porównanie metody symulacyjnej, rzeczywistych obliczeń oraz używanych w teorii przybliżeń

Porównanie wraz z wizualizacją wyników

“Napisz kod w R, który pokaże porównanie metody symulacyjnej, rzeczywistych obliczeń, przybliżenia rozkładem Poissona lub asymptotycznie rozkładem normalnym dla testu pustych cel”

# Wczytanie bibliotek
library(ggplot2)
library(tidyr)

# ==========================================
# 1. PARAMETRY TESTU
# ==========================================
set.seed(123)
n <- 80      # liczba elementów (losowań)
m <- 40      # liczba cel (kategorii)
n_sim <- 10000 # liczba symulacji Monte Carlo

# Ustalenie zakresu wartości k (liczba pustych cel), które chcemy analizować
# (Dla n=80 i m=40 rzadko zobaczymy mniej niż 2 lub więcej niż 12 pustych cel)
k_vals <- 0:15 

# ==========================================
# 2. DEFINICJE OBLICZEŃ
# ==========================================

# A. METODA DOKŁADNA (Exact)
prob_exact <- function(k, n, m) {
  if (k < 0 || k > m - 1) return(0)
  sum_val <- sum(sapply(0:(m - k), function(v) {
    (-1)^v * choose(m - k, v) * (1 - (k + v) / m)^n
  }))
  return(choose(m, k) * sum_val)
}
p_exact <- sapply(k_vals, prob_exact, n = n, m = m)

# B. METODA SYMULACYJNA (Monte Carlo)
# Wielokrotne losowanie i sprawdzanie, ile cel zostało pustych
sim_results <- replicate(n_sim, {
  losowania <- sample(1:m, n, replace = TRUE)
  zajete <- length(unique(losowania))
  return(m - zajete) # liczba pustych
})
# Obliczamy proporcje (prawdopodobieństwo empiryczne)
p_sim <- as.numeric(table(factor(sim_results, levels = k_vals))) / n_sim

# Parametry analityczne dla przybliżeń asymptotycznych
# Wartość oczekiwana E(S)
mu <- m * (1 - 1/m)^n
# Wariancja Var(S)
var_S <- m * (m - 1) * (1 - 2/m)^n + m * (1 - 1/m)^n - m^2 * (1 - 1/m)^(2 * n)
sigma <- sqrt(var_S)

# C. PRZYBLIŻENIE POISSONA
# Zakładamy rozkład Poissona z parametrem lambda = E(S)
p_pois <- dpois(k_vals, lambda = mu)

# D. PRZYBLIŻENIE NORMALNE
# Ze względu na to, że k to zmienna dyskretna, stosujemy poprawkę na ciągłość: P(k-0.5 < X < k+0.5)
p_norm <- pnorm(k_vals + 0.5, mean = mu, sd = sigma) - pnorm(k_vals - 0.5, mean = mu, sd = sigma)

# ==========================================
# 3. ZESTAWIENIE WYNIKÓW W TABELI
# ==========================================
wyniki_df <- data.frame(
  k = k_vals,
  Symulacja = round(p_sim, 5),
  Dokladne = round(p_exact, 5),
  Poisson = round(p_pois, 5),
  Normalny = round(p_norm, 5)
)

cat("=== PORÓWNANIE PRAWDOPODOBIEŃSTW P(S=k) ===\n")

=== PORÓWNANIE PRAWDOPODOBIEŃSTW P(S=k) ===

cat(sprintf("Parametry: n = %d, m = %d, E(S) = %.2f, Var(S) = %.2f\n\n", n, m, mu, var_S))

Parametry: n = 80, m = 40, E(S) = 5.28, Var(S) = 3.19

print(wyniki_df)

    k Symulacja Dokladne Poisson Normalny
1   0    0.0006  0.00100 0.00511  0.00313
2   1    0.0104  0.00961 0.02694  0.01347
3   2    0.0407  0.04145 0.07109  0.04274
4   3    0.1052  0.10635 0.12507  0.09985
5   4    0.1814  0.18205 0.16501  0.17186
6   5    0.2220  0.22106 0.17417  0.21792
7   6    0.1970  0.19762 0.15320  0.20361
8   7    0.1345  0.13323 0.11550  0.14016
9   8    0.0706  0.06885 0.07619  0.07108
10  9    0.0283  0.02757 0.04468  0.02655
11 10    0.0072  0.00862 0.02358  0.00730
12 11    0.0020  0.00211 0.01131  0.00148
13 12    0.0001  0.00041 0.00498  0.00022
14 13    0.0000  0.00006 0.00202  0.00002
15 14    0.0000  0.00001 0.00076  0.00000
16 15    0.0000  0.00000 0.00027  0.00000

# ==========================================
# 4. WIZUALIZACJA W GGPLOT2
# ==========================================
# Przekształcenie danych do formatu "długiego" (long format), wymaganego przez ggplot2
df_long <- pivot_longer(wyniki_df, cols = c(Dokladne, Poisson, Normalny), 
                        names_to = "Metoda", values_to = "Prawdopodobienstwo")

wykres <- ggplot() +
  # Słupki dla symulacji empirycznej (nasz "ground truth" i tło)
  geom_col(data = wyniki_df, aes(x = k, y = Symulacja, fill = "Symulacja"), 
           alpha = 0.3, color = "black") +
  # Linie i punkty dla metod analitycznych
  geom_line(data = df_long, aes(x = k, y = Prawdopodobienstwo, color = Metoda), 
            linewidth = 1, linetype = "dashed") +
  geom_point(data = df_long, aes(x = k, y = Prawdopodobienstwo, color = Metoda), 
             size = 3) +
  scale_fill_manual(name = "Empiryczna", values = c("Symulacja" = "grey50")) +
  scale_color_manual(name = "Aproksymacje", 
                     values = c("Dokladne" = "black", "Poisson" = "#E69F00", "Normalny" = "#56B4E9")) +
  scale_x_continuous(breaks = k_vals) +
  theme_minimal() +
  labs(
    title = "Test pustych cel: Porównanie rozkładów",
    subtitle = sprintf("Liczba kul: %d | Liczba cel: %d | Symulacje: %d", n, m, n_sim),
    x = "Liczba pustych cel (k)",
    y = "Prawdopodobieństwo P(S=k)"
  ) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold", size = 14)
  )

# Wyświetlenie wykresu
print(wykres)

Wniosek: * Pełna zgodność: Metoda dokładna (wzór matematyczny) i symulacja Monte Carlo pokrywają się w 100%. Wykazują, że w losowym układzie najczęściej zostaje około 5 pustych cel.

Ograniczenia przybliżeń: Rozkłady Poissona i Normalny pokazują dobry ogólny trend, ale tracą precyzję na krańcach wykresu (w ogonach rozkładu).
Werdykt: Ponieważ obliczenia dla n=80 i m=40 trwają ułamek sekundy, nie warto iść na skróty. Do rzetelnej analizy i testowania hipotez stosuj wyłącznie metodę dokładną lub symulację, aby uniknąć błędów aproksymacji.

Bibliografia

Hellwig Z. (1965), Test zgodności dla małej próbki, Przegląd Statystyczny, s. 99-112.
Domański C. (2011), Własności testu Davida-Hellwiga, Prace Naukowe Uniwersytetu Ekonomicznego we Wrocławiu, nr 165, s. 41-49.
Csorgo M., Guttman I. (1962), On the Empty Cell Test, Technometrics, Vol. 4, No. 2, s. 235-247.
Kończak G. (2005), On the Modification of David-Hellwig Test, s.3-4.