Problem

Oszacować średnią na podstawie próby

Problem 1

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

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

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

Liczba rugbystów: 623.

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 turnieju RWC’2015 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 2 zawodników pobranych losowo

Powtarzamy eksperyment 1000 razy (dwóch bo dla jednego nie obliczmy wariancji)

s02 <- mks(2, wN)

summary(s02)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    75.0    96.0   103.0   102.6   109.0   126.5
sd(s02)
## [1] 9.334801

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. 
##    89.8   100.1   102.8   102.7   105.4   115.7
sd(s10)
## [1] 4.041939

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

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

s40 <- mks(40, wN)
summary(s40)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   95.33  101.36  102.81  102.75  104.13  110.41
sd(s40)
## [1] 2.028034

Wykres

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

error10 <- true.mean.w  * 0.1

p1 <- all.samples %>%
  pivot_longer(cols = c(s02, 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

Wnioski z eksperymentu

Wszystkie średnie są zbliżone do wartości prawdziwej (to się nazywa nieobciążoność); jeżeli będziemy oceniać wartość prawdziwej średniej na podstawie próby, a naszą ocenę powtórzymy wielokrotnie, to średnia będzie zbliżona do wartości prawdziwej (a nie np. niższa czy wyższa) Ta cecha jest niezależna od wielkości próby.

Jeżeli rośnie liczebność próby to prawdopodobieństwo, że wartość oceniona na podstawie średniej z próby będzie zbliżona do wartości szacowanego parametru rośnie (to się nazywa zgodność)

Jeżeli mamy dwie metody oszacowania parametru obie nieobciążone oraz zgodne, to którą wybrać? Tę która ma mniejszą wariancję. Taką metodę nazywa się bardziej efektywną.

Metodę o której mowa, formalnie funkcję elementów z próby, nazywa się w statystyce estymatorem

Estymator zatem powinien być nieobciążony, zgodny oraz efektywny (czyli mieć małą wariancję). Można matematycznie udowodnić, że jakiś estymator ma tak małą wariancję, że niemożliwe jest wynalezienie czegoś jeszcze bardziej efektywnego. Takim estymatorem średniej w populacji jest średnia z próby…

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   34.00   44.00   45.77   58.00   82.00
sd(k1)
## [1] 14.58705

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. 
##   32.30   42.80   45.70   45.97   49.30   57.90
sd(k10)
## [1] 4.529905

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. 
##   38.75   44.55   46.27   46.19   47.91   52.88
sd(k40)
## [1] 2.320317
## 46 / 2.37
diffMx(k40, true.mean.w)
## [1] 38

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.37   45.15   46.27   46.28   47.47   50.86
sd(k70)
## [1] 1.728121
##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.2965  5.8853  9.5380 11.0716 14.3450 54.8529
sd(k1)
## [1] 7.317069
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. 
##   5.384   9.577  10.842  11.138  12.525  20.804
sd(k10)
## [1] 2.212277

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

k100 <- mks(100, wN)
summary(k100)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    8.81   10.65   11.16   11.17   11.66   13.56
sd(k100)
## [1] 0.7340225

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

k1000 <- mks(1000, wN)
summary(k1000)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.28   10.99   11.15   11.15   11.30   11.78
sd(k1000)
## [1] 0.2247542

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