Statystyka - matematyka, która odróżnia trop od przypadku

Ostatni dzień sezonu, @name. Masz w notesie liczby zwierząt z dziesięciu wodopojów, a dyrekcja rezerwatu prosi o jedną liczbę: ile zwierząt przypada na jeden wodopój. Dziesięć obserwacji ma się skurczyć do jednej cyfry w raporcie, a od tego, którą wybierzesz, zależy, czy raport mówi prawdę.

Drugie pytanie jest jeszcze gorsze. Ważyliście lwy na północy i na południu rezerwatu, po osiem sztuk. Północ wyszła cięższa o dwadzieścia jeden kilogramów. Czy północne lwy naprawdę są większe, czy po prostu traficie akurat na osiem grubszych osobników? Oko tego nie rozstrzygnie. Ja przez trzy sezony byłem przekonany, że gepardy z zachodniej doliny biegają szybciej, i dopiero policzenie tego na papierze pokazało, że różnica mieściła się w zwykłym rozrzucie pomiarów.

Statystyka to zestaw narzędzi do tych dwóch pytań: jak streścić zbiór liczb do kilku sensownych wartości i czy to, co widzisz, nie jest przypadkiem. W tej lekcji przejdziemy oba, od średniej w notesie po prostą, która przewiduje przyszły sezon.

Skąd bierze się

stats

Do liczenia potrzebne są dwie biblioteki. Pierwszą już znasz: NumPy trzyma pomiary w tablicach i liczy na nich szybko, a przyjęło się skracać go do

np

stats

scipy

pyplot

Skoro

stats

from ... import ...

from

scipy

import

stats

from

scipy

import

stats

SyntaxError

Istnieje druga poprawna droga:

import scipy.stats as stats

as

from scipy import stats

1import numpy as np
2from scipy import stats

Po tych dwóch liniach nic się jeszcze nie policzyło i nic nie wypisało. W pamięci wiszą tylko dwie etykiety:

np

stats

Jedna liczba zamiast dziesięciu

Zacznijmy od pytania dyrekcji. Najbardziej oczywiste streszczenie zbioru liczb to średnia arytmetyczna: suma wszystkich pomiarów podzielona przez ich liczbę. Liczy ją funkcja

mean

np.array

populations

1populations = np.array([120, 450, 85, 200, 330, 150, 280, 95, 400, 175])
2
3print(np.mean(populations))
4# 228.5

Wyszło 228.5 zwierzęcia na wodopój i ta liczba jest policzona bezbłędnie. Problem w tym, że nie opisuje rezerwatu. Sprawdź sam: poniżej 228.5 leży sześć wodopojów z dziesięciu, powyżej tylko cztery. Średnia wypadła wyżej niż większość obserwacji, bo dwa wodopoje z tłumami 400 i 450 zwierząt pociągnęły ją w górę. Średnia bierze każdą liczbę jednakowo poważnie, więc każda skrajność ma na nią realny wpływ.

Środek, którego nie da się przekupić

Skoro średnia daje się ciągnąć, potrzebujemy miary, która patrzy nie na wielkość liczb, ale na ich pozycję. Tym jest mediana: ustawiasz pomiary od najmniejszego do największego i bierzesz ten ze środka. Gdy pomiarów jest parzyście, jak u nas dziesięć, mediana to średnia dwóch środkowych. Liczy ją funkcja

median

mean

1print(np.median(populations))
2# 187.5

Posortowany rezerwat wygląda tak: 85, 95, 120, 150, 175, a potem 200, 280, 330, 400, 450. Dwie środkowe wartości to 175 i 200, więc mediana wynosi 187.5. Różnica wobec średniej to czterdzieści jeden zwierząt - i to jest już informacja o kształcie danych, a nie usterka. Gdy średnia jest wyraźnie wyższa od mediany, znaczy to, że w zbiorze siedzi kilka bardzo dużych wartości.

Kiedy przez wodopój przechodzi stado gnu

Teraz zróbmy to, co robi natura raz na sezon. Przy wodopoju, gdzie zwykle było 450 zwierząt, przechodzi migrujące stado gnu i liczymy 4200 sztuk. Jeden pomiar, dziewięć pozostałych bez zmian. Tak wygląda wartość odstająca, po angielsku outlier: obserwacja, która nie jest błędem, tylko rzadkim zdarzeniem, a mimo to potrafi rozsadzić statystykę. Policzmy na nowo oba środki.

1with_herd = np.array([120, 4200, 85, 200, 330, 150, 280, 95, 400, 175])
2
3print(np.mean(with_herd))
4# 603.5
5print(np.median(with_herd))
6# 187.5

Zatrzymaj się przy tych dwóch liniach, bo to najważniejszy moment lekcji. Średnia skoczyła z 228.5 na 603.5, czyli prawie potroiła się przez jedną obserwację. Mediana nie drgnęła ani o jotę: 187.5 przed stadem i 187.5 po stadzie. Powód jest prosty - dla mediany liczy się tylko to, że 4200 leży gdzieś na prawym końcu kolejki. Czy jest to 4200, czy 450, czy 4 miliony, środkowa para wciąż zostaje ta sama. Dlatego to właśnie mediana jest miarą tendencji centralnej odporną na wartości odstające, i to ona trafia do raportu, gdy dane bywają dziwne.

Warto od razu wyprostować trzy pomyłki, które w tym miejscu wracają najczęściej. Średnia arytmetyczna jest miarą tendencji centralnej, ale odporna nie jest - właśnie zobaczyłeś, jak jeden pomiar przesunął ją o 375 jednostek. Wariancja i odchylenie standardowe w ogóle nie są miarami tendencji centralnej: nie mówią, gdzie leży środek, tylko jak szeroko dane są rozrzucone. Co gorsza, obie budują się na kwadratach odległości od średniej, więc na wartości odstające reagują jeszcze gwałtowniej niż sama średnia. Zaraz to policzymy i zobaczysz skalę.

Trzeci środek: dominanta

Miary tendencji centralnej są trzy, więc domknijmy komplet. Dominanta (moda) to wartość, która występuje najczęściej. Ma sens tam, gdzie pomiary się powtarzają: liczba patroli dziennie, liczba młodych w miocie, gatunek najczęściej widziany przy wodopoju. Liczy ją funkcja

mode

stats

mode

count

1patrols = np.array([3, 5, 4, 5, 2, 5, 4])
2
3result = stats.mode(patrols)
4print(result.mode)
5# 5
6print(result.count)
7# 3

Najczęściej wyjeżdżaliście na pięć patroli dziennie i zdarzyło się to trzy razy. Zwróć uwagę, że nie policzyliśmy dominanty na tablicy

populations

Jak bardzo dane są rozrzucone

Dwa rezerwaty mogą mieć identyczną średnią i nie mieć ze sobą nic wspólnego. W pierwszym każdy wodopój zbiera około dwustu zwierząt, w drugim połowa świeci pustkami, a druga połowa pęka w szwach. Do opisania tej różnicy służą miary rozproszenia. Najważniejsza z nich to odchylenie standardowe: mniej więcej typowa odległość pojedynczego pomiaru od średniej. Liczy je

np.std

np.var

1print(np.std(populations))
2# 123.06603918222119
3print(np.var(populations))
4# 15145.25

Przy średniej 228.5 typowy wodopój odchyla się od niej o jakieś 123 zwierzęta - rozrzut jest ogromny i to jest prawdziwy portret rezerwatu. Wariancja 15145.25 niesie tę samą informację, ale w jednostkach podniesionych do kwadratu, czyli w "zwierzętach kwadratowych", których nikt nie umie sobie wyobrazić. Dlatego w raportach podaje się odchylenie standardowe, a wariancję zostawia się wzorom w środku obliczeń. Zapamiętaj też, że

np.std

ddof=1

Obiecałem pokazać, jak te dwie miary reagują na stado gnu. Odchylenie standardowe dla

with_herd

Gdzie mieszka połowa stada

Jest jeszcze trzecie spojrzenie na te same dane, odporne jak mediana, ale bogatsze. Percentyl to wartość, poniżej której leży zadany procent obserwacji: percentyl 25 odcina najniższą ćwiartkę pomiarów, percentyl 95 mówi, poniżej czego mieści się prawie wszystko. Trzy percentyle mają własne nazwy i zwie się je kwartylami: Q1 to percentyl 25, Q2 to percentyl 50 (czyli po prostu mediana), Q3 to percentyl 75. Odstęp między Q3 a Q1 nazywa się rozstępem międzykwartylowym, w skrócie IQR, i opisuje przedział, w którym mieszka środkowa połowa danych. Liczy go gotowa funkcja

iqr

stats

np.ptp

1print(np.percentile(populations, 25))
2# 127.5
3print(np.percentile(populations, 50))
4# 187.5
5print(np.percentile(populations, 75))
6# 317.5
7print(stats.iqr(populations))
8# 190.0
9print(np.ptp(populations))
10# 365

Środkowa połowa wodopojów mieści się między 127.5 a 317.5 zwierzęcia, a Q2 wyszedł identyczny jak wcześniejsza mediana, bo to dwie nazwy tej samej rzeczy. Najciekawsze dzieje się jednak, gdy te same cztery linie puścisz na tablicy

with_herd

Kształt, czyli w którą stronę przechyla się rozkład

Wiesz już, gdzie jest środek i jak szeroki jest rozrzut. Brakuje trzeciej informacji: czy dane rozkładają się symetrycznie wokół środka, czy przechylają się w jedną stronę. Mierzy to skośność, po angielsku skewness, licząca się funkcją

skew

1print(stats.skew(populations))
2# 0.5269796953903054
3print(stats.skew(with_herd))
4# 2.6335940542443765

Zwykły rezerwat ma skośność 0.53, czyli lekki przechył w prawo - i właśnie dlatego średnia 228.5 leżała wyżej niż mediana 187.5. Po dodaniu stada gnu skośność skacze do 2.63, co jest już krzykiem "w tych danych siedzi coś skrajnego". Skośność jest więc szybkim czujnikiem wartości odstających, zanim jeszcze zaczniesz ich szukać po jednej. Istnieje jeszcze siostrzana miara, kurtoza z funkcji

stats.kurtosis

Rozkład normalny, czyli skąd biorą się dane do ćwiczeń

Pomiary w przyrodzie bardzo często układają się w ten sam kształt: najwięcej wyników skupia się w okolicy środka, a im dalej od niego, tym rzadziej coś się trafia. Ten kształt nazywa się rozkładem normalnym i wygląda jak dzwon. Opisują go dwie liczby: środek dzwonu, w SciPy i NumPy nazywany

loc

scale

Nie zawsze masz pod ręką prawdziwe dane, więc przydaje się umiejętność ich wylosowania. Robi to

np.random.normal

loc

scale

size

np.random.seed

round

1np.random.seed(7)
2heights = np.random.normal(loc=500, scale=40, size=200)
3
4print(round(np.mean(heights), 2))
5# 498.88
6print(round(np.std(heights), 2))
7# 39.62

Prosiliśmy o średnią 500 i odchylenie 40, a dostaliśmy 498.88 i 39.62. Ta drobna rozbieżność nie jest błędem - to normalny szum losowania, ten sam, przez który dwie ekipy liczące ten sam rezerwat nigdy nie wracają z identycznym wynikiem. Przy dwustu żyrafach jest niewielki, przy dziesięciu byłby dużo większy. Ten sam mechanizm ma braci od innych kształtów:

np.random.poisson

np.random.exponential

Czy dane naprawdę są normalne

Wiele metod statystycznych zakłada, że dane pochodzą z rozkładu normalnego, więc wypada to sprawdzić, a nie zgadywać z wykresu. Służy do tego test Shapiro-Wilka, wywoływany jako

stats.shapiro

stats

shapiro

1stat, p_value = stats.shapiro(heights)
2
3print(round(stat, 4))
4# 0.9912
5print(round(p_value, 4))
6# 0.2614

P-value wynosi 0.2614, czyli spokojnie powyżej progu 0.05, więc żyrafy przechodzą test - czego się spodziewaliśmy, bo sami wylosowaliśmy je z rozkładu normalnego. Zwróć uwagę na ostrożne sformułowanie: test nie udowodnił normalności, tylko nie znalazł powodu, by ją odrzucić. To subtelna, ale ważna różnica w całej statystyce.

Zobaczmy teraz, jak wygląda wynik, gdy dane naprawdę są krzywe. Puśćmy ten sam test na rezerwacie ze stadem gnu, a próg sprawdźmy tym razem wprost porównaniem

p_value < 0.05

True

False

1stat, p_value = stats.shapiro(with_herd)
2
3print(round(stat, 4))
4# 0.44
5print(p_value < 0.05)
6# True

Statystyka spadła z 0.99 do 0.44, a p-value zeszło poniżej jednej milionowej, więc warunek

p_value < 0.05

Czy różnica dwóch stad jest prawdziwa

Wracamy do drugiego pytania z początku lekcji. Zważyliście po osiem lwów w dwóch częściach rezerwatu i chcecie wiedzieć, czy północ naprawdę hoduje cięższe zwierzęta. Zacznijmy od zapisania pomiarów i policzenia obu średnich narzędziem, które już znasz.

1north = np.array([190, 205, 178, 212, 195, 188, 201, 199])
2south = np.array([172, 168, 185, 176, 180, 165, 178, 174])
3
4print(np.mean(north))
5# 196.0
6print(np.mean(south))
7# 174.75

Różnica wynosi 21.25 kilograma na korzyść północy. Sama w sobie nie znaczy jeszcze nic, bo pomiarów jest tylko osiem na grupę, a lwy potrafią różnić się między sobą o kilkanaście kilogramów. Pytanie brzmi: czy taka różnica mogłaby wyskoczyć z czystego przypadku, gdyby oba stada były w rzeczywistości identyczne?

Odpowiada na to t-test, po polsku test t-Studenta, dostępny jako

stats.ttest_ind

1t_stat, p_value = stats.ttest_ind(north, south)
2
3print(round(t_stat, 4))
4# 4.8177
5print(round(p_value, 6))
6# 0.000273

P-value wynosi 0.000273, czyli 0.03 procenta. Gdyby oba stada były takie same, taka różnica trafiałaby się rzadziej niż raz na trzy tysiące prób - więc uznajemy różnicę za rzeczywistą i możemy wpisać ją do raportu. Zauważ, że średnie same się nie zmieniły; test niczego nie przelicza w danych, tylko dokłada do nich ocenę wiarygodności.

Skoro t-test jest narzędziem często mylonym, ustalmy wprost, do czego nie służy. Nie normalizuje danych - przeskalowaniem pomiarów zajmuje się osobna operacja, w której odejmujesz średnią i dzielisz przez odchylenie standardowe, i poznasz ją przy inżynierii cech. Nie liczy korelacji - od tego są

stats.pearsonr

stats.spearmanr

np.random

np.random.normal

Gdy dane nie są liczbami, tylko kategoriami

T-test potrzebuje liczb, bo liczy średnie. Co jednak, gdy obie zmienne to kategorie? Notujecie porę wyprawy (rano albo wieczorem) i to, czy widziano lwa (tak albo nie), i chcecie wiedzieć, czy pora dnia ma z tym cokolwiek wspólnego. Takie dane układa się w tabelę kontyngencji: wiersze to jedna cecha, kolumny druga, a w środku liczba przypadków. Bada ją test chi-kwadrat, wywoływany jako

stats.chi2_contingency

dof

expected

1observed = np.array([[50, 30],
2                     [20, 40]])
3
4chi2, p_value, dof, expected = stats.chi2_contingency(observed)
5
6print(round(chi2, 4))
7# 10.5292
8print(round(p_value, 4))
9# 0.0012
10print(dof)
11# 1

Rano widziano lwa 50 razy przy 30 nieudanych wyprawach, wieczorem odwrotnie: 20 sukcesów na 40 prób. P-value 0.0012 jest znacznie poniżej 0.05, więc pora dnia i spotkanie lwa nie są niezależne - poranne wyprawy naprawdę mają wyższą skuteczność. W tabeli

expected

Czy dwie liczby idą w parze

Kolejne pytanie z terenu brzmi inaczej niż poprzednie: nie porównujemy dwóch grup, tylko sprawdzamy, czy dwie wielkości zmieniają się razem. Im dłużej siedzimy w ukryciu, tym więcej zwierząt widzimy - ale jak mocny jest ten związek? Mierzy to współczynnik korelacji Pearsona, liczony funkcją

stats.pearsonr

O współczynniku musisz wiedzieć trzy rzeczy. Mieści się zawsze w przedziale od -1 do 1. Jego znak mówi o kierunku: dodatni znaczy "rośnie jedno, rośnie drugie", ujemny "rośnie jedno, maleje drugie". Jego wartość bezwzględna mówi o sile: blisko 1 to bardzo silny związek, blisko 0 to brak związku. I najważniejsze zastrzeżenie: Pearson mierzy wyłącznie zależność liniową, czyli to, jak dobrze punkty układają się wzdłuż prostej.

1hours = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
2seen = np.array([4, 3, 9, 12, 10, 18, 15, 21, 19, 26])
3
4corr, p_value = stats.pearsonr(hours, seen)
5
6print(round(corr, 4))
7# 0.9526
8print(round(p_value, 6))
9# 2.1e-05

Wyszło 0.9526, w zaokrągleniu 0.95. Czytamy to jednym zdaniem: silna dodatnia korelacja liniowa - im dłuższa obserwacja, tym więcej zwierząt, a punkty leżą blisko prostej. Gdyby wyszło -0.95, siła byłaby taka sama, ale kierunek odwrotny: dłuższa obserwacja oznaczałaby mniej zwierząt. Brak korelacji to okolice zera, powiedzmy od -0.1 do 0.1, a nie 0.95. I ostatnia pułapka: 0.95 nie mówi nic o zależności nieliniowej, bo Pearson takich zależności po prostu nie mierzy - wysoka wartość świadczy właśnie o tym, że związek jest bliski prostej.

Kiedy Pearson nie widzi oczywistego związku

Zobaczmy tę granicę na żywo. Populacja pewnego gatunku owadów mniej więcej podwaja się z każdym tygodniem, więc związek między numerem tygodnia a liczebnością jest doskonały - tyle że wykładniczy, a nie prostoliniowy. Do takich przypadków służy korelacja Spearmana, wywoływana jako

stats.spearmanr

1week = np.array([1, 2, 3, 4, 5, 6, 7, 8])
2insects = np.array([2, 4, 9, 20, 45, 100, 230, 520])
3
4corr, p_value = stats.pearsonr(week, insects)
5print(round(corr, 4))
6# 0.8171
7
8rho, p_value = stats.spearmanr(week, insects)
9print(round(rho, 4))
10# 1.0

Te same dane, dwie różne odpowiedzi. Pearson mówi 0.82, bo krzywa wykładnicza kiepsko udaje prostą. Spearman mówi 1.0, czyli związek jest doskonały - każdy kolejny tydzień daje wynik wyższy niż poprzedni, bez jednego wyjątku. Nie znaczy to, że jeden ze współczynników kłamie: one mierzą co innego. Pearson pyta "czy to prosta", Spearman pyta "czy to rośnie". Gdy podejrzewasz zależność nieliniową, sięgaj po Spearmana.

Ile warta jest Twoja średnia

Każda średnia policzona z próbki jest tylko oszacowaniem. Zmierzyłeś dwieście żyraf, a nie wszystkie na świecie, więc prawdziwa średnia gatunku leży gdzieś w pobliżu Twojej, ale niekoniecznie dokładnie na niej. Przedział ufności podaje ten "gdzieś w pobliżu" w liczbach: przedział 95-procentowy to zakres, o którym możemy z 95-procentową pewnością powiedzieć, że mieści prawdziwą średnią. Potrzebne są do tego dwa narzędzia:

stats.sem

stats.t.interval

len

loc

scale

1mean = np.mean(heights)
2sem = stats.sem(heights)
3ci_95 = stats.t.interval(0.95, len(heights) - 1, loc=mean, scale=sem)
4
5print(round(sem, 3))
6# 2.808
7print(round(ci_95[0], 2), round(ci_95[1], 2))
8# 493.34 504.42

Zamiast suchego "średnia wysokość żyrafy to 498.88 cm" możesz teraz napisać zdanie mocniejsze: prawdziwa średnia mieści się z 95-procentową pewnością między 493.34 a 504.42 cm. Sama średnia się nie zmieniła - zyskała tylko uczciwą informację o swojej niepewności. Gdybyś zmierzył nie dwieście, a dziesięć żyraf, przedział byłby wielokrotnie szerszy, bo mniejsza próbka to większa niepewność. To najprostszy sposób, by raport przestał udawać większą precyzję, niż faktycznie ma.

Regresja liniowa w czterech krokach

Zostało ostatnie i najbardziej praktyczne pytanie: skoro stado rośnie, ile zwierząt będzie w przyszłym sezonie? Korelacja mówi tylko "rośnie razem", a Ty potrzebujesz konkretnej liczby. Daje ją regresja liniowa: dopasowanie do punktów prostej o równaniu

y = a * x + b

a

slope

y

x

b

intercept

y

x

stats.linregress

Krok pierwszy to przygotowanie danych

x

y

1seasons = np.array([1, 2, 3, 4, 5, 6])
2herd = np.array([34, 41, 45, 52, 58, 63])

Nic się jeszcze nie policzyło, ale dane mają teraz właściwy kształt: sześć sezonów i sześć odpowiadających im pomiarów, ustawione parami w tej samej kolejności. To jedyny moment, w którym decydujesz, co jest przyczyną, a co skutkiem -

x

y

Krok drugi to wywołanie

stats.linregress(x, y)

1result = stats.linregress(seasons, herd)

I znów cisza - nic się nie wypisało, bo wynik wylądował w zmiennej

result

slope

intercept

rvalue

pvalue

stderr

result

result.slope

Krok trzeci to odczytanie

slope

intercept

1print(round(result.slope, 2))
2# 5.8
3print(round(result.intercept, 2))
4# 28.53

Prosta ma równanie

y = 5.8 * x + 28.53

Krok czwarty to sprawdzenie R-squared, bo prostą da się dopasować nawet do chmury zupełnie przypadkowych punktów i nikt Cię przed tym nie ostrzeże.

rvalue

**

y

pvalue

1print(round(result.rvalue ** 2, 4))
2# 0.9964
3print(round(result.pvalue, 6))
4# 5e-06

R-kwadrat wynosi 0.9964, czyli prosta wyjaśnia ponad 99 procent zmienności liczebności stada - przewidywaniu 69 zwierząt można zaufać. Gdyby wyszło 0.24, jak w danych skaczących w górę i w dół bez trendu, wynik

slope

slope

intercept

Co zabierasz z obozu

• Miary tendencji centralnej to średnia, mediana i dominanta. Odporna na wartości odstające jest mediana - średnią jeden pomiar potrafi przesunąć o setki jednostek, a wariancja i odchylenie standardowe w ogóle nie opisują środka, tylko rozproszenie, i reagują na skrajności jeszcze mocniej.
• Rozproszenie opisujesz przez

  np.std

  i

  np.var

  , a odporniej przez kwartyle i

  stats.iqr

  . Zwykły rozstęp

  np.ptp

  ginie od jednej wartości odstającej.
• Moduł statystyczny importujesz w kolejności

  from

  ,

  scipy

  ,

  import

  ,

  stats

  . Zapis w drugą stronę to

  SyntaxError

  .
• T-test (

  stats.ttest_ind

  ) służy do porównania średnich dwóch grup. Nie normalizuje danych, nie liczy korelacji i nie generuje liczb losowych.
• Test chi-kwadrat sprawdza zależność między cechami kategorycznymi z tabeli kontyngencji.
• Korelacja Pearsona równa 0.95 oznacza silną dodatnią korelację liniową. Silna ujemna to -0.95, brak korelacji to okolice zera, a zależności nieliniowej Pearson w ogóle nie mierzy - od tego jest Spearman.
• Test normalności wywołujesz z sześciu części:

  stats

  ,

  .

  ,

  shapiro

  ,

  (

  ,

  data

  ,

  )

  , a wynik czytasz z p-value wobec progu 0.05.
• Przedział ufności z

  stats.sem

  i

  stats.t.interval

  dokłada do średniej uczciwą miarę niepewności.
• Regresja liniowa ma cztery kroki w stałej kolejności: przygotuj dane

  x

  i

  y

  , wywołaj

  stats.linregress(x, y)

  , odczytaj

  slope

  i

  intercept

  , sprawdź R-kwadrat.

Statystyka nie dokłada do notesu ani jednej nowej obserwacji, @name - daje Ci tylko prawo powiedzieć na odprawie, że trop jest prawdziwy, a nie że tak Ci się w terenie wydawało.