1 Wprowadzenie

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.

1.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ść.

1.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).

b <- density(runif(10))
   f <- CDF(b)
   f(0.5)
## [1] 0.3749553
 plot(f)

1.4 Przeczytaj

Przeczytaj artykuł naukowy “Kernel-smoothed cumulative distribution function estimation with akdensity” autorstwa Philippe Van Kerm.

1.5 Zadanie

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).

dane_cale<-diagnoza %>% na.omit()

wykres1 <- ggplot(data=dane_cale, aes(x=gp64, fill=plec)) +
  labs(title="Dochód miesięczny netto", x = "Dochód", y = "Wynagrodzenia") +
  geom_density(alpha = 0.5) + 
  ggtitle("Wykres gęstości dochodu według płci i edukacji") +
  facet_wrap(~ eduk4_2013) +
  theme(
    text = element_text(color = "#003366"), 
    plot.title = element_text(size = 10),
    axis.title = element_text(size = 10, color = '#003366'),
    axis.title.y = element_text(vjust = .5, angle = 90),
    axis.title.x = element_text(hjust = .5),
    legend.position = "bottom",
    legend.title = element_text(size = 7)
  )

wykres1

Wnioski:

W miarę wzrostu poziomu wykształcenia widoczne jest przesunięcie dochodów w prawo, co sugeruje, że osoby z wyższym wykształceniem osiągają wyższe dochody. Osoby z wykształceniem wyższym i policealnym mają rozkłady gęstości skoncentrowane na wyższych dochodach niż osoby z niższymi poziomami edukacji. W każdej grupie wykształcenia większość osób ma dochody w niskich przedziałach, zwłaszcza wśród osób z wykształceniem podstawowym lub zasadniczym zawodowym. Wykształcenie wyższe i policealne wydaje się jednak być bardziej zróżnicowane pod względem dochodów, co sugeruje większy rozstrzał w wynagrodzeniach dla tej grupy.

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

dane <- data.frame(dochod=diagnoza$gp64, 
                   wojewodztwo=diagnoza$wojewodztwo,
                   plec=diagnoza$plec
                   )
dane <- dane %>% na.omit()

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

dochod<-dane$dochod
1.5.1 Rozkład gęstości

Typy jądra funkcji gęstości

par(mfrow = c(2, 1))
wykres1<-oszac_gestosci<-density(dochod, kernel="gaussian", bw=2)
plot(oszac_gestosci, main="Default(Gaussian)")
wykres2<-oszac_gestosci<-density(dochod, kernel="gaussian", bw=3)
plot(oszac_gestosci, main="Default(Gaussian)")

par(mfrow = c(2, 1))
wykres3<-oszac_gestosci<-density(dochod, kernel="rectangular", bw=1)
plot(oszac_gestosci, main="Rectangular")
wykres4<-oszac_gestosci<-density(dochod, kernel="rectangular", bw=2.5)
plot(oszac_gestosci, main="Rectangular")

par(mfrow = c(2, 1))
wykres5<-oszac_gestosci<-density(dochod, kernel="triangular", bw=2)
plot(oszac_gestosci, main="Triangular")

wykres6<-oszac_gestosci<-density(dochod, kernel="triangular", bw=2.5)
plot(oszac_gestosci, main="Triangular")

par(mfrow = c(2, 1))
wykres7<-oszac_gestosci<-density(dochod, kernel="epanechnikov", bw=3)
plot(oszac_gestosci, main="Epanechnikov")

wykres8<-oszac_gestosci<-density(dochod, kernel="epanechnikov", bw=0.5)
plot(oszac_gestosci, main="Epanechnikov")

Podział ze względu na płeć:

plot_gestosc_g <- ggplot(dane, aes(x = dochod, fill = plec)) +
  stat_density(alpha = 0.6, position = "identity") +  
  labs(title = "Estymacja gęstości dochodu (jądro Gaussowskie)", 
       x = "Dochód", 
       y = "Gęstość") +
  scale_fill_manual(values = c("#4C72B0", "#DD8452")) + 
  theme_minimal() + 
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, color = "#333333"),  # Stylizacja tytułu
    axis.title = element_text(size = 12, color = "#333333"), 
    axis.text = element_text(color = "#333333"),
    legend.position = "top",
    legend.title = element_blank(),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()
  ) +
  guides(fill = guide_legend(title = "Płeć")) 

plot_gestosc_g

plot_gestosc_r <- ggplot(dane, aes(x = dochod, fill = plec)) +
  stat_density(kernel = "rectangular", alpha = 0.6, position = "identity") +
  labs(
    title = "Estymacja gęstości dochodu z jądrem prostokątnym",
    x = "Dochód",
    y = "Gęstość"
  ) +
  scale_fill_manual(values = c("#4C72B0", "#DD8452")) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, color = "#333333"),
    axis.title = element_text(size = 12, color = "#333333"),
    axis.text = element_text(color = "#333333"),
    legend.position = "top",
    legend.title = element_blank(),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()  
  ) +
  guides(fill = guide_legend(title = "Płeć"))

plot_gestosc_r

plot_gestosc_t <- ggplot(dane, aes(x = dochod, fill = plec)) +
  stat_density(kernel = "triangular", alpha = 0.6, position = "identity") +
  labs(
    title = "Estymacja gęstości dochodu z jądrem trójkątnym",
    x = "Dochód",
    y = "Gęstość"
  ) +
  scale_fill_manual(values = c("#4C72B0", "#DD8452")) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, color = "#333333"),
    axis.title = element_text(size = 12, color = "#333333"),
    axis.text = element_text(color = "#333333"),
    legend.position = "top",
    legend.title = element_blank(),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()  
  ) +
  guides(fill = guide_legend(title = "Płeć"))

plot_gestosc_t

plot_gestosc_e <- ggplot(dane, aes(x = dochod, fill = plec)) +
  stat_density(kernel = "epanechnikov", alpha = 0.6, position = "identity") +
  labs(
    title = "Estymacja gęstości dochodu z jądrem epaschnikova",
    x = "Dochód",
    y = "Gęstość"
  ) +
  scale_fill_manual(values = c("#4C72B0", "#DD8452")) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5, color = "#333333"),
    axis.title = element_text(size = 12, color = "#333333"),
    axis.text = element_text(color = "#333333"),
    legend.position = "top",
    legend.title = element_blank(),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()  
  ) +
  guides(fill = guide_legend(title = "Płeć"))

plot_gestosc_e

Wnioski:

Na podstawie przedstawionych wykresów można zauważyć, że zmiana jądra (Gaussowskie, prostokątne, trójkątne, Epanechnikova) nie prowadzi do istotnych różnic w kształcie estymacji gęstości dochodu. Rozkład gęstości dochodów pozostaje praktycznie taki sam niezależnie od zastosowanego jądra. Wszystkie wykresy pokazują najwyższą gęstość w niskich wartościach dochodu, a następnie gwałtowny spadek. Najwyższe wartości gęstości są skupione blisko zera, co sugeruje, że większość obserwacji dotyczy osób z niskim dochodem. Ponieważ różnice między rozkładami są niewielkie, wybieramy standardowe jądro Gaussowskie

1.5.2 Selekcja pasma wygładzania (PDF)
den_nrd0 <- density(dane$dochod, bw = "nrd0")
plot(den_nrd0, main = "Domyślne pasmo nrd0", col="purple")

den_scotta <- density(dane$dochod, bw = "nrd")
plot(den_scotta, main = "Pasmo (Reguła Scotta)", col='green')

den_ucv <- density(dane$dochod, bw = "ucv")
## Warning in bw.ucv(x): odnotowano wartość minimalną na jednym z końców zakresu
# Kernel density plot
plot(den_ucv, main = "Pasmo UCV (nieobciążona cross-walidacja)", col='blue')

den_bcv <- density(dane$dochod, bw = "bcv")

# Kernel density plot
plot(den_bcv, main = "Pasmo BCV (obciążona cross-walidacja)", col='orange')

den_sj <- density(dane$dochod, bw = "SJ")

# Kernel density plot
plot(den_sj, main = "Pasmo (Metoda Sheather & Jones)", col='red')

Wnioski:

Metoda BCV przesadza z wygładzaniem danych (zbyt duża szerokość pasma), a metody SJ oraz UCV powodują zbyt mocne dopasowanie się do krzywej (zbyt mała szerokość pasma). Rozsądnym wyborem będzie zatem domyślne pasmo nrd0 lub reguła Scotta.

ggplot(dane, aes(x = dochod, fill = plec)) +
  stat_density(bw = "nrd0", alpha = 0.6, position = "identity") + 
  labs(
    title = "Estymacja gęstości dochodu z domyślnym pasmem (nrd0)",
    x = "Dochód",
    y = "Gęstość"
  ) +
  scale_fill_manual(values = c("#4C72B0", "#DD8452")) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5, color = "#333333"),
    axis.title = element_text(size = 12, color = "#333333"),
    axis.text = element_text(color = "#333333"),
    legend.position = "top",
    legend.title = element_blank(),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank()
  ) +
  guides(fill = guide_legend(title = "Płeć"))

1.5.3 Rozkład skumulowanej gęstości (dystrybuanty)

Typy jądra

gestosc_gaussian <- density(dane$dochod, kernel = "gaussian", bw = 4)
gestosc_rectangular <- density(dane$dochod, kernel = "rectangular", bw = 4)
gestosc_epanechnikov <- density(dane$dochod, kernel = "epanechnikov", bw = 4)

dystrybuanta_gaussian <- approxfun(gestosc_gaussian$x,
                          cumsum(gestosc_gaussian$y)/sum(gestosc_gaussian$y))
dystrybuanta_rectangular <- approxfun(gestosc_rectangular$x,
                          cumsum(gestosc_rectangular$y)/sum(gestosc_rectangular$y))
dystrybuanta_epanechnikov <- approxfun(gestosc_epanechnikov$x,
                          cumsum(gestosc_epanechnikov$y)/sum(gestosc_epanechnikov$y))

x_vals <- seq(min(dane$dochod), max(dane$dochod), length.out = 100)
par(mfrow = c(3, 1))

plot(x_vals, dystrybuanta_gaussian(x_vals), type = "l", col = "blue", lwd = 2,
     main = "CDF - jądro Gaussowskie", 
     xlab = "Dochód", ylab = "Dystrybuanta",
     ylim = c(0, 1)) 
grid(col = "gray80") 

plot(x_vals, dystrybuanta_rectangular(x_vals), type = "l", col = "darkgreen", lwd = 2,
     main = "CDF - jądro prostokątne", 
     xlab = "Dochód", ylab = "Dystrybuanta",
     ylim = c(0, 1))
grid(col = "gray80")

plot(x_vals, dystrybuanta_epanechnikov(x_vals), type = "l", col = "purple", lwd = 2,
     main = "CDF - jądro Epanechnikova", 
     xlab = "Dochód", ylab = "Dystrybuanta",
     ylim = c(0, 1))
grid(col = "gray80")

1.5.4 Selekcja pasma wygładzania (CDF)
dystrybuanta_nrd0 <- approxfun(den_nrd0$x,
                          cumsum(den_nrd0$y)/sum(den_nrd0$y))

dystrybuanta_scotta <- approxfun(den_scotta$x,
                          cumsum(den_scotta$y)/sum(den_scotta$y))

dystrybuanta_ucv <- approxfun(den_ucv$x,
                          cumsum(den_ucv$y)/sum(den_ucv$y))
dystrybuanta_bcv <- approxfun(den_bcv$x,
                          cumsum(den_bcv$y)/sum(den_bcv$y))

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

x_vals <- seq(min(dane$dochod), max(dane$dochod), length.out = 100)
plot(x_vals, dystrybuanta_nrd0(x_vals), pch=1, 
     main = "CDF - pasmo nrd0", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "purple", lwd = 2)

plot(x_vals, dystrybuanta_scotta(x_vals), pch=1, 
     main = "CDF - pasmo Scotta", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "green", lwd = 2)

plot(x_vals, dystrybuanta_ucv(x_vals), pch=1, 
     main = "CDF - pasmo UCV", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "blue", lwd = 2)

plot(x_vals, dystrybuanta_bcv(x_vals), pch=1, 
     main = "CDF - pasmo BCV", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "orange", lwd = 2)

plot(x_vals, dystrybuanta_SJ(x_vals), pch=1, 
     main = "CDF - pasmo Sheather & Jones", 
     xlab = "Dochód", ylab = "Dystrybuanta", 
     col = "red", lwd = 2)

Wnioski:

Estymacja gęstości (PDF) i skumulowanej gęstości (CDF) służą do analizy danych, do wizualizacji rozkładu i struktury danych w sposób intuicyjny i szczegółowy. Wybór odpowiedniego typu jądra i szerokości pasma jest kluczowy dla uzyskania wiarygodnych wyników, ponieważ może znacząco wpłynąć na wygląd i interpretację wykresu.

W przypadku estymacji gęstości (PDF), wybór szerokości pasma odgrywa szczególnie istotną rolę. Mniejsze pasmo umożliwia dokładniejsze odwzorowanie drobnych fluktuacji w danych, pokazując lokalne zmiany w zagęszczeniu, ale może prowadzić do bardziej “szorstkiego” wykresu, który uwidacznia nawet niewielkie zmiany, w tym wpływ wartości odstających. Z kolei większe pasmo wygładza wykres, ujmując ogólny kształt rozkładu, ale redukując widoczność mniejszych różnic.

Skumulowana funkcja gęstości (CDF) jest bardziej odporna na wybór szerokości pasma i typu jądra, ponieważ pokazuje skumulowane prawdopodobieństwo, które zawsze wzrasta do 1. W rezultacie, różnice w parametrach estymacji mają tu mniejszy wpływ na ogólny kształt wykresu. Wartości odstające, które mogą mieć duży wpływ na PDF, mają mniejsze znaczenie w CDF, ponieważ wpływają jedynie na końcowe wartości dystrybuanty, a proces sumowania prawdopodobieństw “wygładza” różnice, które mogłyby być bardziej widoczne na wykresie PDF.