Problem

Oszacować średnią na podstawie próby

Problem 1

Plik rwc2015.csv zawiera dane (w tym dotyczące wagi) dla 623 rugbystów uczestniczących w turnieju o Puchar Świata w 2016 roku.

r <- read.csv("rwc2015.csv", sep = ';', dec = ",",  header=T, na.string="NA")

w <- as.vector(na.omit(r$weight))
wN <- length(w)

Obliczamy (prawdziwą) średnią, odchylenie standardowe i współczynnik zmienności:

summary(w)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    65.0    93.0   103.0   102.8   113.0   145.0
## 102.8 kg
sd(w)
## [1] 12.92366
## 12.9 kg
true.mean.w <- mean(w)
## 46.24 lat
sd(w, na.rm = T)
## [1] 12.92366
max.err <- 0.1 * true.mean.w 
max.err.sd <- sd(w, na.rm = T) 

max.err.sd / true.mean.w * 100
## [1] 12.57213

Czyli średnio rugbysta na RWC2015 ważył 102.7961477 kg a odchylenie standardowe wyniosło 12.9236635 kg.

Wykres (rozkład jest dwumodalny; bo w rugby są dwie grupy zawodników, wcale nie wszyscy > 110 kg):

q1 <- ggplot(r, aes(x=weight)) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_histogram(binwidth=2, alpha=.5, fill="steelblue") +
  ggtitle("Rozkład wagi zawodników")
q1
## Warning: Removed 1 rows containing non-finite values (stat_bin).

Szacujemy średnią na podstawie 1 zawodnika pobranego losowo

Powtarzamy eksperyment 1000 razy

s1 <- mks(1, wN)

summary(s1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    65.0    92.0   102.0   102.2   113.0   145.0       4
sd(s1)
## [1] NA

szacujemy średnią na podstawie 10 zawodników pobranych losowo

Powtarzamy eksperyment 1000 razy

s10 <- mks(10, wN)

summary(s10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    84.5   100.3   103.0   102.8   105.5   114.0
sd(s10)
## [1] 4.06763

szacujemy średnią na podstawie 40 zawodników pobranych losowo

Uwaga: 40 zawodników to około 10% całego zbioru. Powtarzamy eksperyment 1000 razy

s40 <- mks(40, wN)
summary(s40)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   96.53  101.22  102.80  102.71  104.12  108.72
sd(s40)
## [1] 2.067104

Wykres

all.samples <- data.frame(s1, s10, s40)

error10 <- true.mean.w  * 0.1

p1 <- all.samples %>%
  pivot_longer(cols = c(s1, s10, s40), names_to = 'k', values_to = 'v') %>%
  ggplot(aes(x=v)) +
  facet_wrap(~ k) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_vline(xintercept = true.mean.w - max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w + max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w - error10, colour="orange", size=.4) +
  geom_vline(xintercept = true.mean.w + error10, colour="orange", size=.4) +
  geom_histogram(binwidth=1, alpha=.5, fill="steelblue") +
  ggtitle("rozkład średniej wagi rugbystów w zależności od wielkości próby")
p1
## Warning: Removed 4 rows containing non-finite values (stat_bin).

Problem 2

Plik kandydaci_ws_2018_3.csv zawiera dane (w tym dotyczące wieku) dla ponad 7000 kandydatów do Sejmików Wojewódzkich w wyborach samorządowych w 2018 roku.

r <- read.csv("kandydaci_ws_2018_3.csv", sep = ';', dec = ",",  header=T, na.string="NA")

w <- as.vector(na.omit(r$wiek))
wN <- length(w)

Dokładniej to kandydatów jest 7076.

Obliczamy (prawdziwą) średnią, odchylenie standardowe i współczynnik zmienności:

summary(w)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   34.00   46.00   46.24   58.00   91.00
true.mean.w <- mean(w)
sd(w)
## [1] 14.61312
max.err <- 0.1 * true.mean.w 
max.err.sd <- sd(w) 
max.err.sd / true.mean.w * 100
## [1] 31.60347

Czyli średnio kandydat miał 46.2389768 lat a odchylenie standardowe wieku wyniosło 14.613121 lat.

Wykres (rozkład znowu jest dwumodalny z jakiś powodów):

q1 <- ggplot(r, aes(x=wiek)) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_histogram(binwidth=2, alpha=.5, fill="steelblue") +
  ggtitle("Rozkład wieku kandydatów")
q1

Szacujemy średnią na podstawie 1 kandydata pobranego losowo

Powtarzamy eksperyment 1000 razy

k1 <- mks(1, wN)

summary(k1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.00   35.00   46.00   46.31   58.00   83.00
sd(k1)
## [1] 14.52132

Szacujemy średnią na podstawie 10 kandydatów pobranych losowo

Powtarzamy eksperyment 1000 razy.

k10 <- mks(10, wN)
summary(k10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   33.40   42.90   46.10   46.15   49.30   61.80
sd(k10)
## [1] 4.576852

Szacujemy średnią na podstawie 40 kandydatów pobranych losowo

Uwaga: 40 kandydatów to ok 0.6% całości. Powtarzamy eksperyment 1000 razy.

40/wN * 100
## [1] 0.5652911
k40 <- mks(40, wN)

summary(k40)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   39.20   44.64   46.23   46.18   47.73   53.35
sd(k40)
## [1] 2.312611
## 46 / 2.37
diffMx(k40, true.mean.w)
## [1] 46

Szacujemy średnią na podstawie 70 kandydatów pobranych losowo

Uwaga: 70 kandydatów to około ok 1% całości (1000 powtórzeń)

70/wN * 100
## [1] 0.9892595
k70 <- mks(70, wN)
summary(k70)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   41.51   45.03   46.34   46.24   47.39   52.34
sd(k70)
## [1] 1.762296
##diffMx(k70, true.mean.w)

Wykres

all.samples <- data.frame(k1, k10, k40, k70)

error10 <- true.mean.w  * 0.1

p1 <- all.samples %>%
  pivot_longer(cols = c(k1, k10, k40, k70), names_to = 'k', values_to = 'v') %>%
ggplot(aes(x=v)) +
  facet_wrap(~ k) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_vline(xintercept = true.mean.w - max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w + max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w - error10, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w + error10, colour="red", size=.4) +
  geom_histogram(binwidth=1, alpha=.5, fill="steelblue") +
  ggtitle("rozkład średniej wieku kandydatów w zależności od wielkości próby")
p1

Problem 3

Plik airline_delay_causes.csv (pobrany z http://www.transtats.bts.gov/OT_Delay/OT_DelayCause1.asp) zawiera dane dotyczące opóźnień w ruchu lotniczym w USA w latach 2003–2021. Są to wskaźniki miesięczne dla każdego przewoźnika/lotniska docelowego. Wykorzystujemy zmienne arr_flights: Number of flights which arrived at the airport; arr_del15: Number of flights delayed (>= 15minutes late); arr_delay: Total time (minutes) of delayed flights

dc <- read.csv("airline_delay_causes.csv", 
              sep = ',', dec = ".",  header=T, na.string="NA")

Obliczamy wskaźnik: średnie opóźnienie miesięczne (pomijamy lotniska o liczbie lotów mniejszej niż 30/miesiąc czyli 1/dzień):

dc0 <- dc %>% select (arr_flights, arr_del15, arr_delay) %>%
  filter (arr_flights > 30) %>%
  mutate (del1 = arr_delay/arr_flights, 
          del15p = arr_del15/ arr_flights)

Liczymy prawdziwe średnie itp:

nr0 <- nrow(dc0)

w <- as.vector(na.omit(dc0$del1))

wN <- length(w)
summary(w)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   6.062   9.724  11.149  14.648 172.194
true.mean.w <- mean(w)
sd(w)
## [1] 7.22357
max.err <- 0.1 * true.mean.w 
max.err.sd <- sd(w) 
max.err.sd / true.mean.w * 100
## [1] 64.79395

średnia wartość wskaźnika linia/lotnisko wynosi zatem 11.1485249 (odchylenie standardowe 7.22357). Wszytkich wskaźników jest 272324

Rozkład wskaźnika jest cośkolwiek skośny:

q1 <- ggplot(dc0, aes(x=del1)) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_histogram(binwidth=1, alpha=.5, fill="steelblue") +
  ggtitle("Miesięczny wskaźnik opóźnień przewoźnik/lotnisko")
q1

Szacujemy średnią na podstawie 1 wartości pobranej losowo

Wykonujemy jak zwykle 1000 powtórzeń

k1 <- mks(1, wN)
summary(k1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   6.262   9.855  11.417  15.168  59.683
sd(k1)
## [1] 7.375452
sample.m1 <- mean(k1)

szacujemy średnią na podstawie 10 wartości pobranych losowo

k10 <- mks(10, wN)
summary(k10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.899   9.609  10.880  11.160  12.495  22.128
sd(k10)
## [1] 2.318091

szacujemy średnią na podstawie 100 wartości pobranych losowo

k100 <- mks(100, wN)
summary(k100)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   9.051  10.706  11.168  11.177  11.623  14.967
sd(k100)
## [1] 0.7101456

szacujemy średnią na podstawie 1000 wartości pobranych losowo

k1000 <- mks(1000, wN)
summary(k1000)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.51   11.00   11.15   11.15   11.31   12.03
sd(k1000)
## [1] 0.2283832

Wykres

error10 <- true.mean.w  * 0.1

all.samples <- data.frame(k1, k10, k100, k1000)
p1 <- all.samples %>%
  pivot_longer(cols = c(k1, k10, k100, k1000), names_to = 'k', values_to = 'v') %>%
  ggplot(aes(x=v)) +
  facet_wrap(~ k) +
  geom_vline(xintercept = true.mean.w, colour="forestgreen", size=.4) +
  geom_vline(xintercept = true.mean.w - max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w + max.err.sd, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w - error10, colour="red", size=.4) +
  geom_vline(xintercept = true.mean.w + error10, colour="red", size=.4) +
  geom_histogram(binwidth=.5, alpha=.5, fill="steelblue") +
  ggtitle("rozkład miesięcznego wskaźnika opóźnień linia-lotnisko (minuty)")
p1

Wniosek: precyzja wnioskowania zwiększa się wraz z liczebnością próby; tym szybciej im rozproszenie w populacji generalnej jest mniejsze. Żeby z dużą dokładnością wnioskować o średniej dla dużej populacji wcale nie trzeba pobierać dużej próby. W ostatnim przykładzie wystarczyło 100/1000 (czyli 0,04–0,4% liczebności populacji).

Żeby było ciekawiej istnieje rozkład teoretyczny, który z dobrą dokładnością opisuje rozkłady empiryczne będące wynikiem powyższej zabawy. Ten rozkład (zwany normalnym) zależy tylko od dwóch parametrów: średniej i rozproszenia, gdzie średnia będzie równa (prawdziwej) średniej w populacji a rozproszenie wartości rozproszenia w populacji podzielonej przez pierwiastek z wielkości próby.

Dla próby 100-elementowej wygląda to tak:

this.sample.size <- 100
k1000m <- true.mean.w
k1000sd <- max.err.sd / sqrt(this.sample.size)
this.binwd <- .1
p1b <- data.frame(k100) %>%
  ggplot(aes(x=k100)) +
  geom_histogram(binwidth= this.binwd, alpha=.5, fill="steelblue") +
  stat_function(fun = function(x) 
  {dnorm(x, mean = k1000m, sd = k1000sd) * sample.size * this.binwd},
  color="red")
p1b

dla próby 1000-elementowej:

this.sample.size <- 1000
k1000m <- true.mean.w
k1000sd <- max.err.sd / sqrt(this.sample.size)
this.binwd <- .1
p1b <- data.frame(k1000) %>%
  ggplot(aes(x=k1000)) +
  geom_histogram(binwidth= this.binwd, alpha=.5, fill="steelblue") +
  stat_function(fun = function(x) 
  {dnorm(x, mean = k1000m, sd = k1000sd) * sample.size * this.binwd},
  color="red")
p1b

Prawda, że wynik jest całkiem dobry? Teoretyczność czerwonej krzywej polega na tym, że ona zawsze będzie identyczna, podczas gdy histogram będzie różny. Gdybyśmy powtórzyli nasz eksperyment (generowania 1000 losowych prób przypominam), to zapewne trochę by się różnił, bo byśmy wylosowali inne wartości do prób. Ta teoretyczna abstrakcja nazywa się prawdopodobieństwem. Rzucając monetą 1000 razy spodziewamy się po 500 orłów i reszek, co w modelu matematycznym będzie opisane jak: prawdopodobieństwo wyrzucenia orła wynosi 0,5. Rzucanie monetą to bardzo prosty eksperyment; nasz z liczeniem średniej wskaźnika opóźnień jest dużo bardziej skomplikowany więc miło jest się dowiedzieć, że używając czerwonej krzywej można łatwo obliczyć jak bardzo prawdopodobne jest na przykład popełnienie błędu większego niż 10%, albo większego niż 5 minut. Albo jak duża powinna być próba żeby ten błąd był nie większy niż 5 minut.

Na tym polega wnioskowanie statystyczne: można określić z jakim prawdopodobieństem oszacujemy coś-tam (w tym przypadku wartość średnią) na podstawie próby z dokładnością ileś-tam. Dokładność i prawdopodobieństwo jest założone z góry. Chcemy dokładniej lub z większą szansą, że się nie rąbniemy, to musimy zwiększyć wielkość próby. Idea zatem jest całkiem prosta, ale matematyczna teoria stojąca za tym już niestety nie. Więc na tym zakończymy :-)