Nieklasyczne metody statystyki

Nieparametryczna estymacja gęstości/dystrybuanty - Raport 1 + 2

1 Krótkie wprowadzenie teoretyczne

Kernel Density Estimation (KDE), czyli estymacja jądrowa gęstości, jest nieparametryczną metodą oszacowania funkcji gęstości rozkładu zmiennej losowej na podstawie dostępnych danych. W przeciwieństwie do metod parametrycznych, takich jak estymacja rozkładu normalnego, KDE nie wymaga przyjęcia konkretnej formy rozkładu danych. Zamiast tego, KDE buduje estymację funkcji gęstości w sposób gładki, opierając się na rzeczywistych obserwacjach z próby. Dla każdego punktu danych umieszczana jest funkcja jądrowa (kernel), najczęściej o kształcie zbliżonym do krzywej dzwonowej, np. rozkładu Gaussa, której zadaniem jest przypisanie “wagi” obserwacjom w pobliżu danego punktu. Szerokość pasma (bandwidth) kontroluje stopień wygładzenia całkowitej funkcji gęstości.

2 Ćwiczenie 1.

Uruchom demo estymatora funkcji gęstości kernel. Zmieniaj zarówno dane wejściowe, jak i opcje estymatora - szerokość pasma oraz rodzaj funkcji jądrowej. Czy widzisz istotne różnice w oszacowaniu?

#install.packages("remotes") #tylko raz! potem #
#remotes::install_github("hericks/KDE") #tylko raz! potem #
#library(KDE)
#shiny_kde() 

Kilka wniosków:

-Małe pasmo (wąska szerokość): Gdy szerokość pasma jest bardzo mała, estymacja gęstości staje się bardzo szczegółowa i dokładnie dopasowuje się do poszczególnych obserwacji, co powoduje, że estymator jest podatny na lokalne wahania i “szumy” w danych.

-Duże pasmo (szeroka szerokość): Przy dużym paśmie estymator gęstości jest wygładzony, co sprawia, że lokalne szczegóły są mniej widoczne, ale za to lepiej oddaje ogólny kształt rozkładu danych. Przy zbyt dużym paśmie może jednak dojść do nadmiernego wygładzenia (underfitting), gdzie szczegóły rozkładu są zupełnie zgubione, a estymator przestaje odzwierciedlać realne zmiany w danych.

-Przy mniejszej liczbie danych estymator KDE jest bardziej wrażliwy na dobór szerokości pasma i rodzaju jądra. Przy małych próbach każde jądro ma duży wpływ na wynik, a małe pasmo generuje wyraźny szum.

-Przy większej liczbie obserwacji KDE działa bardziej stabilnie. Większe próby lepiej odzwierciedlają rzeczywisty rozkład danych, a estymacja staje się bardziej precyzyjna. Dobór szerokości pasma i jądra ma mniejszy wpływ na wynik przy dużej liczbie danych.

-Użycie walidacji krzyżowej (cross-validation) do doboru szerokości pasma okazało się skuteczne w uzyskaniu dobrze zbalansowanego estymatora, który zachowuje kształt rozkładu bez nadmiernego wygładzania.

-Jądro Gaussowskie daje chyba najbardziej intuicyjne i gładkie estymacje.

3 Ćwiczenie 2.

Wykorzystując dowolną funkcję R do estymacji funkcji gęstości oszacuj jej przebieg dla wynagrodzeń (zbiór danych salaries) strażaków w San Francisco. Wykorzystaj metody graficzne dostępne w pakiecie ggplot2. Mile widziane przekroje oraz odpowiedzi na pytania badawcze zadane na wstępie.

salaries <- read.csv("https://github.com/kflisikowski/ds/raw/master/Salaries.csv")

3.1 Histogram

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")

hist(firefighter_salaries$TotalPay, 
     main = "Total Pay for Firefighters", 
     xlab = "Pay (in dollars)", 
     col = "lightblue", 
     border = "black")

abline(v = mean(firefighter_salaries$TotalPay, na.rm = TRUE), col = "red", lty = "dashed")
abline(v = median(firefighter_salaries$TotalPay, na.rm = TRUE), col = "blue", lty = "solid")
legend("topright", 
       legend = c("Mediana", "Średnia"), 
       col = c("blue", "red"), 
       lty = c("solid", "dashed"))

Większość strażaków zarabia w przedziale około 100,000 - 175,000 dolarów rocznie, co stanowi największe skupisko wartości na wykresie. Rozkład ten jest względnie symetryczny, co sugeruje, że wartości wynagrodzeń są równomiernie rozłożone wokół centralnej wartości. Z tego wykresu wynika, że mediana i średnia są bardzo zbliżone do siebie, co sugeruje, że w rozkładzie nie występują duże wartości odstające, które mogłyby znacząco zniekształcić średnią.

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")
par(mfrow = c(2, 2))
hist(firefighter_salaries$TotalPay, 
     main = "Total Pay for Firefighters, default breaks", 
     xlab = "Pay (in dollars)", 
     col = "lightblue", 
     border = "black")
hist(firefighter_salaries$TotalPay, 
     main = "Total Pay for Firefighters, breaks=100", 
     xlab = "Pay (in dollars)", 
     breaks = 100, 
     col = "lightblue", 
     border = "black")
hist(firefighter_salaries$TotalPay, 
     main = "Total Pay for Firefighters, breaks=1000", 
     xlab = "Pay (in dollars)", 
     breaks = 1000, 
     col = "lightblue", 
     border = "black")

3.2 Wykres pudełkowy

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")

boxplot(firefighter_salaries$TotalPay, 
        main = "Total Pay for Firefighters", 
        ylab = "Pay (in dollars)", 
        col = "lightblue", 
        border = "black")

abline(h = median(firefighter_salaries$TotalPay, na.rm = TRUE), col = "blue", lty = "dashed")
abline(h = mean(firefighter_salaries$TotalPay, na.rm = TRUE), col = "red", lty = "dotted")
legend("topright", 
       legend = c("Mediana", "Średnia"), 
       col = c("blue", "red"), 
       lty = c("dashed", "dotted"))

3.3 Wykres ggplot

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")

ggplot(firefighter_salaries, aes(x = TotalPay)) +
  geom_density(fill = "skyblue", alpha = 0.5, color = "darkblue") +
  geom_vline(xintercept = median(firefighter_salaries$TotalPay, na.rm = TRUE), color = "blue", linetype = "dashed", linewidth = 1) +
  geom_vline(xintercept = mean(firefighter_salaries$TotalPay, na.rm = TRUE), color = "red", linetype = "dotted", linewidth = 1) +
  labs(title = "Distribution of Total Pay for Firefighters",
       x = "Pay (in dollars)", y = "Density") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  ) 

Rozkład wynagrodzeń strażaków jest względnie symetryczny, z głównym skupiskiem w okolicach 150,000 dolarów. Dodatkowy szczyt przy niższych wartościach może wynikać z różnych poziomów doświadczenia lub stanowisk w grupie strażaków.

3.4 Estymacja nieliniowej zależności

Poniżej wykorzystano metodę locpoly() z pakietu KernSmooth do estymacji nieliniowej zależności.

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")

fit <- locpoly(x = firefighter_salaries$TotalPay, 
               degree = 0, 
               bandwidth = 10000) %>% as_tibble()

ggplot(firefighter_salaries, aes(x = TotalPay)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", color = "black", alpha = 0.5) +
  geom_line(data = fit, aes(x = x, y = y), color = "purple", linewidth = 1) +
  labs(
    title = "Kernel Density Estimation of Total Pay for Firefighters",
    x = "Total Pay (in dollars)", 
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

3.5 Modyfikacje - estymacja nieliniowej zależności

firefighter_salaries <- subset(salaries, JobTitle == "Firefighter")

fit_1 <- locpoly(firefighter_salaries$TotalPay, bandwidth = 5000, degree = 0) %>% as_tibble()
fit_2 <- locpoly(firefighter_salaries$TotalPay, bandwidth = 10000, degree = 1) %>% as_tibble()
fit_3 <- locpoly(firefighter_salaries$TotalPay, bandwidth = 20000, degree = 2) %>% as_tibble()

ggplot(firefighter_salaries, aes(x = TotalPay)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", color = "black", alpha = 0.5) +
  geom_line(data = fit_1, aes(x = x, y = y), color = "purple", linewidth = 1, linetype = "solid", show.legend = TRUE) +
  geom_line(data = fit_2, aes(x = x, y = y), color = "red", linewidth = 1, linetype = "dashed", show.legend = TRUE) +
  geom_line(data = fit_3, aes(x = x, y = y), color = "blue", linewidth = 1, linetype = "dotdash", show.legend = TRUE) +
  labs(
    title = "Comparison of Kernel Density Estimation",
    x = "Total Pay (in dollars)", 
    y = "Density"
  ) +
  scale_color_manual(
    name = "Bandwidth and Degree",
    values = c("purple", "red", "blue"),
    labels = c("Bandwidth=5000, Degree=0", "Bandwidth=10000, Degree=1", "Bandwidth=20000, Degree=2")
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12)
  )

4 Wprowadzenie - raport 2

Cel: wyznaczenie obszaru ufności dla dystrybuanty nieznanego rozkładu, a nie tylko oszacowania parametrów, od jakich zależą jej wartości.

• Brzegi tego obszaru są wykresami funkcji „przedziałami stałych” (funkcji schodkowych).

• Jeżeli przy wyznaczaniu pasma ufności dla dystrybuanty otrzymamy lewy kraniec przedziału będący liczbą ujemną, to zastępujemy ją przez zero.

• Jeżeli otrzymamy prawy kraniec przedziału większy od jedności, to przyjmujemy, że jest on równy jeden.

• Określenie obszaru ufności dla dystrybuanty w przedstawiony sposób polega na wyznaczeniu przedziałowego oszacowania dla każdej wartości dystrybuanty.

4.1 Funkcja CDF

Funkcja w programie R odpowiedzialna za estymację to np. CDF z pakietu spatstat. CDF jest metodą ogólną, z metodą dla klasy “gęstość”.

Oblicza ona skumulowaną funkcję rozkładu, której gęstość prawdopodobieństwa została oszacowana i zapisana w obiekcie f. Obiekt f musi należeć do klasy “gęstość” i zazwyczaj zostałby uzyskany z wywołania funkcji gęstość.

4.2 Funkcja kde

Pakiet R o nazwie snpar zawiera kilka uzupełniających metod statystyki nieparametrycznej, w tym test kwantylowy, test trendu Coxa-Stuarta, test przebiegów, test normalnego wyniku, estymację jądra PDF i CDF, estymację regresji jądra i test jądra Kołmogorowa-Smirnowa.

Funkcja kde zawiera obliczanie zarówno nieparametrycznego estymatora jądra funkcji gęstości prawdopodobieństwa (PDF) jak i funkcji rozkładu skumulowanego (CDF).

4.3 Zadanie 1

Posłużymy się zbiorem danych diagnozy społecznej.

Na jego podstawie Twoim zadaniem jest oszacowanie rozkładu “p64 Pana/Pani wlasny (osobisty) dochod miesieczny netto (na reke)” według województw/płci. Postaraj się oszacować zarówno rozkład gęstości jak i skumulowanej gęstości (dystrybuanty).

data("diagnoza")
data("diagnozaDict")

dane <- data.frame(dochod=diagnoza$gp64, 
                   plec=diagnoza$plec, 
                   woj=diagnoza$wojewodztwo)

dane <- dane %>% 
  drop_na(dochod, woj, plec)

dane$dochod <- as.numeric(dane$dochod)

Histogram rozkładu dochodu według płci

ggplot(dane, aes(x = dochod, fill = plec)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 30) + 
  labs(
    title = "Rozkład dochodów miesięcznych netto według płci",
    x = "Dochód netto (PLN)",
    y = "Liczba osób"
  ) +
  facet_wrap(~ plec) + 
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "none" 
  )

4.3.1 Rozkład gęstości

4.3.1.1 **Typ jądra

Poniżej wykonano estymację jądrową gęstości (Kernel Density Estimation - KDE) dla zmiennej dochod z wykorzystaniem różnych funkcji jądrowych (kernels) oraz tej samej szerokości pasma (bandwidth = 3). Wykorzystano trzy różne typy jąder: Gaussa (domyślne), Epanechnikova oraz prostokątne (rectangular). Wybór jądra wpływa na gładkość i dokładność estymacji. Jądro Gaussa daje najbardziej gładką estymację, podczas gdy jądro prostokątne jest bardziej sztywne i mniej wygładza dane.

d_gaussian <- density(dane$dochod, kernel = "gaussian", bw = 3)
d_epanechnikov <- density(dane$dochod, kernel = "epanechnikov", bw = 3)
d_rectangular <- density(dane$dochod, kernel = "rectangular", bw = 3)

par(mfrow = c(3, 1))
plot(d_gaussian,main = "Default")
plot(d_epanechnikov, main = "Epanechnikov")
plot(d_rectangular, main = "Rectangular")

Z rozróżnieniem na płeć:

plot_g <- ggplot(dane, aes(x=dochod, fill=plec))+
  stat_density()+
  labs(title="Estymacja gęstości z jądrem Gaussa")+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "none" 
  )
plot_g

plot_e <- ggplot(dane, aes(x=dochod, fill=plec))+
  stat_density(kernel = "epanechnikov")+
  labs(title="Estymacja gęstości z jądrem Epanechnikova")+
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "none" 
  )
plot_e

plot_t <- ggplot(dane, aes(x=dochod, fill=plec))+
  stat_density(kernel = "rectangular")+
  labs(title="Estymacja gęstości z jądrem prostokątnym")+
    theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    legend.position = "none" 
  )
plot_t

Wniosek: Ponieważ różne typy jąder dają podobne wyniki, do dalszej analizy można użyć domyślnego jądra Gaussowskiego. Wybór ten upraszcza obliczenia, jednocześnie zapewniając stabilne i dokładne estymacje rozkładów dochodów według województw i płci.

par(mfrow = c(2, 2))
d_0 <- density(dane$dochod, bw = "nrd0")
plot(d_0, main = "Pasmo nrd0", col="red")

d_scott <- density(dane$dochod, bw = "nrd")
plot(d_scott, main = "Pasmo Scotta", col='yellow')

d_ucv <- density(dane$dochod, bw = "ucv")

## Warning in bw.ucv(x): minimum occurred at one end of the range

plot(d_ucv, main = "Pasmo UCV", col="green")

d_SJ <- density(dane$dochod, bw = "SJ")
plot(d_SJ, main = "Pasmo Sheathera & Jonesa", col="blue")

Wnioski:

Pasmo UCV i Pasmo Sheathera & Jonesa (SJ) ujawniają najwięcej szczegółów i lokalnych wariacji w rozkładzie dochodów, ale mogą wprowadzać „szum”, który utrudnia ogólną interpretację danych. Są one odpowiednie, jeśli potrzebujemy szczegółowej analizy lokalnych zmian w rozkładzie.

Pasmo nrd0 i Pasmo Scotta zapewniają gładszy obraz rozkładu, który jest bardziej odporny na szum i fluktuacje lokalne, co jest korzystne, gdy interesuje nas bardziej ogólny kształt rozkładu. Pasmo Scotta daje najbardziej wygładzony wynik, idealny dla dużych zestawów danych.

ggplot(dane, aes(x=dochod, fill=plec))+
  stat_density(bw = "nrd")+
  labs(title="Estymacja gęstości z pasmem Scotta z podziałem na płeć")

4.3.2 Oszacowanie skumulowanej gęstości

Dystrybuanta to funkcja, która opisuje prawdopodobieństwo, że zmienna losowa przyjmie wartość mniejszą lub równą danej wartości x.

4.3.2.1 Typ jądra

dystrybuanta_g <- approxfun(d_gaussian$x,
                          cumsum(d_gaussian$y)/sum(d_gaussian$y))

dystrybuanta_e <- approxfun(d_epanechnikov$x,
                          cumsum(d_epanechnikov$y)/sum(d_epanechnikov$y))

dystrybuanta_r <- approxfun(d_rectangular$x,
                          cumsum(d_rectangular$y)/sum(d_rectangular$y))

x_vals <- seq(min(dane$dochod), max(dane$dochod), length.out = 100)

par(mfrow = c(3, 1))
plot(x_vals, dystrybuanta_g(x_vals), pch=1,
     main = "CDF dla jądra Gaussowskiego", 
     xlab = "Dochód", ylab = "Dystrybuanta")
plot(x_vals, dystrybuanta_e(x_vals), pch=1,
     main = "CDF dla jądra Epanechnikova", 
     xlab = "Dochód", ylab = "Dystrybuanta")
plot(x_vals, dystrybuanta_r(x_vals), pch=1,
     main = "CDF dla jądra prostokątnego", 
     xlab = "Dochód", ylab = "Dystrybuanta")

Wszystkie trzy dystrybuanty (dla różnych jąder) mają bardzo podobny kształt. Każdy z wykresów pokazuje, że większość wartości dochodów znajduje się w dolnym zakresie, co powoduje szybki wzrost dystrybuanty na początku, a następnie wypłaszczenie. Oznacza to, że dla większości populacji dochody są stosunkowo niskie.

4.3.2.2 Selekcja pasma

dystrybuanta_0 <- approxfun(d_0$x,
                          cumsum(d_0$y)/sum(d_0$y))

dystrybuanta_s <- approxfun(d_scott$x,
                          cumsum(d_scott$y)/sum(d_scott$y))

dystrybuanta_ucv <- approxfun(d_ucv$x,
                          cumsum(d_ucv$y)/sum(d_ucv$y))

dystrybuanta_SJ <- approxfun(d_SJ$x,
                          cumsum(d_SJ$y)/sum(d_SJ$y))

x_vals <- seq(min(dane$dochod), max(dane$dochod), length.out = 100)

plot(x_vals, dystrybuanta_0(x_vals), pch=1, 
     main = "CDF dla pasma nrd0", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "red", lwd = 2)

Szerokość pasma ma ograniczony wpływ na wyniki estymacji dystrybuanty, co obserwujemy na wykresach, gdzie krzywe CDF dla różnych pasm i jąder są do siebie zbliżone. Zarówno estymacja gęstości, jak i skumulowanej gęstości są wartościowymi narzędziami w analizie rozkładu danych, pozwalając na dokładniejsze zrozumienie struktury danych. Różne jądra i metody wyboru pasma mogą jednak przynieść subtelne różnice w wynikach, co podkreśla znaczenie przemyślanego doboru tych parametrów, szczególnie w analizach statystycznych, gdzie precyzja ma kluczowe znaczenie.