Nieklasyczne metody statystyki

Nieparametryczna estymacja dystrybuanty

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

1.3 Przykład 1.

   b <- density(runif(10))
   f <- CDF(b)
   f(0.5)

## [1] 0.603946

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

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 wg płci

ggplot(dane, aes(x = dochod, fill = plec)) +
  geom_histogram(position="dodge", alpha = 0.5)+
  labs(title = "Rozkład dochodów miesięcznych netto według płci")+
  facet_wrap(~ plec)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.5.1 Rozkład gęstości

Estymacja gęstości, w przeciwieństwie do histogramów, tworzy bardziej płynny obraz rozkładu danych, nie dzieląc ich na sztywne kategorie. Dodatkowo uwypuklają nam się różnice w kształcie rozkładu między grupami.

1.5.1.1 Typ jądra

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")
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")
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")
plot_t

Wniosek: Różne typu jądra dają podobne wyniki, zatem do dalszej analizy wystarczy użyć jądra domyślnego - Gausowskiego.

1.5.1.2 Selekcja pasma

d_0 <- density(dane$dochod, bw = "nrd0")
plot(d_0, main = "Pasmo nrd0", col="yellowgreen")

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

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

## Warning in bw.ucv(x): odnotowano wartość minimalną na jednym z końców zakresu

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

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

Wniosek: Szerokość pasma wpływa na wyniki oszacowania gestości. Metody SJ i UCV mają zbyt małą szerokośc pasma, zatem są nadmiernie dopasopwuje sie do krzywej.

Optymalnym wyborem może być reguła Scotta lub domyślna.

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ć")

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

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

1.5.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 = "yellowgreen", lwd = 2)

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

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

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

Wniosek: Szerokość pasma znikomo wpływa na wyniki oszacowania dystrybuanty.

Estymacja gęstości i skumulowanej gęstości są potężnymi narzędziami do analizy danych. Różne typy jąder i metody doboru pasma mogą prowadzić do różnych rezultatów, co podkreśla znaczenie starannego wyboru metod analizy w badaniach statystycznych. Natomiast analizowane dane mają wartości odstające, które znacznie wpłynęły na analizę gęstości i dystrybuanty. Spowodowały one, że różnice pomiędzy różnymi pasmami czy jądrami, nie wpływaja istonie na analizę ( z wyjątkiem doboru pasma dla gęstości ).