Raport 1

Ć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?

Mniejsza szerokość pasma sprawia, że wykres pokazuje więcej szczegółów, ale może być „zbyt dopasowany” do danych, pokazując przypadkowe fluktuacje. Większa szerokość wygładza wykres, ale może ukrywać mniejsze różnice w danych. Różne typy jądra, jak Gaussian czy Epanechnikov, wpływają na to, jak kształtuje się gęstość – mogą bardziej lub mniej wygładzać dane w poszczególnych miejscach.

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

table(salaries_firefighters$Year)
## 
## 2014 
##  738

Wszystkie dane na temat wynagrodzeń strażaków pochodzą z roku 2014, co uniemożliwia sprawdzenie jak płace zmieniały się w czasie.

library(ggplot2)
ggplot(salaries_firefighters, aes(x = TotalPay)) +
  geom_density(fill = "cornflowerblue", alpha = 0.5) +
    labs(title = "Kernel Density Estimate for Firefighters' Total Pay",
       x = "Total Pay (in dollars)",
       y = "Density") +
  theme_minimal()

Szczyt rozkładu znajduje się w przedziale około 140,000 - 160,000 USD, co sugeruje, że większość strażaków może liczyć na zarobki w tym przedziale. Spadek po 160,000 USD oznacza, że wyższe wynagrodzenia są rzadsze, a sam wykres przypomina lekko rozkład normalny, w którym to większość pracowników zarabia kwoty równe medianie, a im wynagrodzenie bardziej różni się od mediany, tym mniej osób ma takie zarobki. Można jednak zauważyć, że nie jest to do końca rozkład normalny, swoim wyglądem przypomina również rozkład prawostronnie skośny, co co oznacza, że zauważalna jest obecność pracowników z bardzo wysokim wynagrodzeniem całkowitym, ale występuje to rzadziej niż niskie wynagrodzenie .

hist(salaries_firefighters$TotalPay,main="Total Pay", xlab="Pay (in dollars)")
abline(v = mean(salaries_firefighters$TotalPay),lty="dashed")
abline(v = median(salaries_firefighters$TotalPay))
legend("topright", legend=c("Mediana","Średnia"),lty=c("solid","dashed"))

hist(salaries_firefighters$TotalPay,main="Total Pay, Zoomed-in", xlab="Pay (in dollars)", xlim=c(60000,270000), breaks=100)

benf <- as.numeric(salaries_firefighters$Benefits)
hist(benf,main="Benefits", xlab="Pay (in dollars)")
abline(v = mean(benf),lty="dashed")
abline(v = median(benf))
legend("topright", legend=c("Mediana","Średnia"),lty=c("solid","dashed"))

Największa częstotliwość (około 500 osób) dotyczy przedziału wynagrodzeń bliskiego 40,000 dolarów, co sugeruje, że większość badanych osób otrzymuje benefity bliskie tej wartości. Histogram jest asymetryczny a średnia jest nieco mniejsza od mediany, z długim ogonem po lewej stronie, co jest typowe dla lewoskośnego rozkładu. Można więc przypuszczać, że wynagrodzenia są skupione wokół mediany, a lewoskośność wykresu wskazuje na mniejszą liczbę pracowników zarabiających wynagrodzenia znacznie poniżej średniej.

Raport 2

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

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

Rozkład gęstości

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

Typ jądra z rozróżnieniem płci

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

Wszystkie testowane jądra dały zbliżone wyniki, dlatego do dalszej analizy użyte zostanie jądro Gaussa

Selekcja pasma

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

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): minimum occurred at one end of the range
plot(d_ucv, main = "Pasmo UCV", col="cornflowerblue")

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

Metody UCV oraz Sheathera&Jonesa mają zbyt małe pasma, co powoduje powstanie nadmiernie dopasowanej krzywej. Zbyt duża szerokość pasma natomiast powoduje powstanie krzywej nadmiernie wygładzonej. Z tego powodu najlepsza byłaby metoda Scotta lub nrd0.

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

Oszacowanie skumulowanej gęstości

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

### 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 = "pink", 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 = "cornflowerblue", lwd = 2)

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