Statystyka i prawdopodobieństwo — od rzutów monetą do wnioskowania bayesowskiego

Prawdopodobieństwo jako ramy matematyczne

Statystyka i prawdopodobieństwo to nie ta sama dyscyplina. Prawdopodobieństwo to czysta matematyka: mając model, jakie wyniki są prawdopodobne? Statystyka odwraca to pytanie: mając wyniki, jaki model jest prawdopodobny? Symulacje w tym numerze znajdują się na granicy między tymi dwoma podejściami — każda pozwala zobaczyć, jak teoretyczny wynik wyłania się z powtarzanych losowych eksperymentów, czyniąc abstrakcję konkretną.

To, co je łączy, to losowość ze strukturą. Rzuty monetą są nieprzewidywalne pojedynczo, jednak 10 000 rzutów tworzy histogram nieodróżnialny od krzywej normalnej. Losowy spacer pijaka po mieście to chaos, jednak ułamek czasu spędzonego w każdej dzielnicy zbiega do dokładnej wartości. Losowość i porządek nie są przeciwieństwami — to to samo zjawisko widziane w różnych skalach czasu.

Warstwa 1: uniwersalna zbieżność — centralne twierdzenie graniczne

Symulacja centralnego twierdzenia granicznego

Centralne twierdzenie graniczne (CTG) jest prawdopodobnie najważniejszym wynikiem w całej statystyce. Mówi ono, że suma dużej liczby niezależnych zmiennych losowych — niezależnie od ich indywidualnych rozkładów — zbiega do rozkładu normalnego. Rzuć kostką raz, a otrzymasz płaski rozkład. Uśrednij dwa rzuty kostką, a pojawiają się piki. Uśrednij trzydzieści, a otrzymasz krzywą dzwonową tak gładką, że tablice rozkładu normalnego dają odpowiedzi dokładne do szóstego miejsca po przecinku.

Dlatego rozkład normalny pojawia się wszędzie w naturze: wzrost ludzi, błędy pomiarowe, stopy zwrotu z akcji, wyniki IQ i losowa prędkość cząsteczek gazu — wszystko to wyłania się z sumowania wielu małych niezależnych wkładów. CTG jest matematycznym powodem, dla którego rozkład Gaussa jest nie tylko wygodą, ale prawem wielkich liczb.

Klasyczne CTG (Lindeberg-Lévy):
  X₁, X₂, ..., Xₙ — n.i.o. o średniej μ, wariancji σ²
  S̄ₙ = (X₁ + X₂ + ... + Xₙ) / n  (średnia z próby)

  Gdy n → ∞:
    √n · (S̄ₙ − μ) / σ  →  N(0, 1)  (według rozkładu)

  Równoważnie:
    S̄ₙ ~ N(μ, σ²/n)  w przybliżeniu, dla dużych n

Błąd standardowy średniej:
  SE = σ / √n
  (zmniejszenie SE o połowę wymaga 4× więcej próbek)

Nierówność Berry'ego-Esseena (tempo zbieżności):
  sup_x |P(√n(S̄ₙ−μ)/σ ≤ x) − Φ(x)| ≤ C · ρ / (σ³ · √n)
  ρ = E[|X − μ|³]  (trzeci moment absolutny)
  C ≈ 0.4748  (dokładna stała, Shevtsova 2010)

Kluczowa obserwacja: zbieżność jest szybsza, gdy ρ/σ³ jest małe
  — rozkłady symetryczne zbiegają szybciej niż skośne
  — lżejsze ogony zbiegają szybciej niż ciężkie

Rozkłady źródłowe w symulacji:
  Jednostajny U(0,1):   μ=0.5, σ²=1/12,  szybkie CTG
  Wykładniczy(1):       μ=1,   σ²=1,      umiarkowane (skośny)
  Bimodalny:            μ=0.5, σ²≈0.25,   wolne (bimodalny)
  Poissona(λ=1):        μ=1,   σ²=1,      umiarkowane

Symulacja centralnego twierdzenia granicznego pozwala wybrać spośród pięciu rozkładów źródłowych — płaskiego, wykładniczego, bimodalnego, Poissona i prawoskośnego — i obserwować w czasie rzeczywistym, jak gromadzą się histogramy średnich z próby. Przeciągnij suwak rozmiaru próby od n = 1 do n = 200 i zobacz, jak rozgrywa się zbieżność: średnie z próby rozkładu bimodalnego wyglądają na rozłożone normalnie już przy n = 30, podczas gdy rozkład o ciężkich ogonach potrzebuje n = 100, aby zbiec się w pełni.

Warstwa 2: aktualizacja przekonań — wnioskowanie bayesowskie

Symulacja wnioskowania bayesowskiego

Statystyka klasyczna pyta: mając hipotezę, jakie jest prawdopodobieństwo zaobserwowanych danych? Statystyka bayesowska pyta odwrotnie: mając dane, jakie jest prawdopodobieństwo, że hipoteza jest prawdziwa? Brzmi to jak to samo pytanie, ale dają one różne odpowiedzi, a mylenie ich jest źródłem wielu błędów w literaturze naukowej.

Twierdzenie Bayesa dostarcza reguły aktualizacji: pomnóż prawdopodobieństwo a priori przez wiarygodność danych, a następnie znormalizuj. Każdy dowód przekształca rozkład prawdopodobieństwa nad hipotezami, przyciągając go do spójnych wyjaśnień. Zacznij od pełnej niepewności, a po wystarczającej liczbie obserwacji nawet niejasny rozkład a priori skoncentruje się ostro wokół prawdziwej wartości.

Twierdzenie Bayesa:
  P(H | E) = P(E | H) · P(H) / P(E)

  H = hipoteza (np. pacjent ma chorobę)
  E = dowód (np. wynik testu jest pozytywny)
  P(H)       = prawdopodobieństwo a priori H
  P(E | H)   = wiarygodność: P(test pozytywny | choroba)
  P(E)       = wiarygodność marginalna (stała normalizująca)
  P(H | E)   = prawdopodobieństwo a posteriori H przy danym E

Prawo prawdopodobieństwa całkowitego:
  P(E) = P(E|H)·P(H) + P(E|¬H)·P(¬H)

Przykład testu medycznego (w symulacji):
  Częstość choroby:        P(H)       = 0.01   (1 na 100)
  Czułość testu:            P(E|H)     = 0.99   (odsetek prawdziwie pozytywnych)
  Specyficzność testu:      P(¬E|¬H)   = 0.95   (odsetek prawdziwie negatywnych)
  → Odsetek fałszywie pozytywnych: P(E|¬H) = 0.05

  P(H|E) = 0.99 × 0.01 / [0.99×0.01 + 0.05×0.99]
         = 0.0099 / (0.0099 + 0.0495)
         = 0.0099 / 0.0594
         ≈ 0.167  ← tylko 16,7% szans na chorobę!

Dlaczego to nieintuicyjne:
  Fałszywie pozytywne wyniki dominują, gdy częstość choroby jest niska.
  Test "wygląda" na dokładny w 99%, ale prawdopodobieństwo a posteriori wynosi tylko 17%.
  Nazywa się to błędem zaniedbania częstości bazowej.

Symulacja wnioskowania bayesowskiego wizualizuje rozkłady a priori, wiarygodności i a posteriori jako krzywe, które można przeciągać. Scenariusz testu medycznego pokazuje, dlaczego nawet test dokładny w 99% ma niską wartość predykcyjną wyniku pozytywnego, gdy schorzenie jest rzadkie — wynik, który zaskakuje zarówno studentów, jak i klinicystów. Uruchom test z monetą, aby obserwować, jak rozkład a posteriori zaostrza się od niepewności 50/50 w kierunku prawdziwego obciążenia po każdym zaobserwowanym rzucie.

Warstwa 3: kontrintuicyjne prawdopodobieństwo — paradoks urodzin

Symulacja paradoksu urodzin

Ile osób potrzeba w pokoju, zanim dwie z nich będą miały ten sam dzień urodzin? Większość ludzi zgaduje wartość w okolicach 180 — mniej więcej połowę z 365. Poprawna odpowiedź brzmi 23. To drastyczne niedoszacowanie to właśnie paradoks urodzin i ujawnia ono, jak słabo ludzka intuicja radzi sobie z kombinatorycznym wybuchem porównań parami.

Przy 23 osobach istnieje C(23,2) = 253 różnych par. Każda para niezależnie ma szansę 1/365 na wspólny dzień urodzin. Prawdopodobieństwo, że żadna z 253 par nie ma wspólnych urodzin, wynosi w przybliżeniu e⁻²⁵³/³⁶⁵ ≈ 0,50. Dla 30 osób rośnie ono do 70%, dla 50 osób do 97%, a dla 70 osób do 99,9%.

Dokładne prawdopodobieństwo, że co najmniej dwie osoby mają wspólne urodziny:
  n = liczba osób,  d = 365 dni

  P(zgodność | n) = 1 − P(brak zgodności | n)
  P(brak zgodności | n) = d/d · (d−1)/d · (d−2)/d · ... · (d−n+1)/d
                        = d! / [(d−n)! · dⁿ]

Przybliżenie na podstawie wzoru Stirlinga (n ≪ d):
  P(brak zgodności) ≈ exp(−n(n−1) / (2d))
  P(zgodność)       ≈ 1 − exp(−n²/(2d))  dla n ≪ d

Kluczowe wartości:
  n = 23 → P(zgodność) ≈ 50,7%
  n = 30 → P(zgodność) ≈ 70,6%
  n = 50 → P(zgodność) ≈ 97,0%
  n = 70 → P(zgodność) ≈ 99,9%

Uogólnienie — znajdowanie kolizji w zbiorze d elementów:
  Oczekiwana liczba próbek przed pierwszą kolizją ≈ √(π·d/2)
  Dla d = 365: √(π·365/2) ≈ 23,9  ("granica urodzinowa")

Znaczenie kryptograficzne:
  Funkcja skrótu tworząca n-bitowe skróty ma d = 2ⁿ możliwych wyjść.
  Kolizji można się spodziewać po ≈ 2^(n/2) skrótach.
  SHA-256 (n=256): kolizja po ~2¹²⁸ skrótach — obliczeniowo niewykonalne.

Symulacja paradoksu urodzin uruchamia jednocześnie dokładną krzywą analityczną i żywy eksperyment Monte Carlo. Widok pokoju pokazuje 365 miejsc rozmieszczonych w okręgu, które zapełniają się w miarę wchodzenia ludzi i migają na czerwono po wykryciu kolizji. Uruchom 10 000 prób natychmiast przyciskiem trybu wsadowego i obserwuj, jak częstość empiryczna dokładnie pokrywa się z krzywą teoretyczną.

Warstwa 4: dopasowywanie linii — regresja liniowa

Symulacja metody najmniejszych kwadratów

Mając rozrzut punktów danych, regresja liniowa znajduje linię, która minimalizuje sumę kwadratów pionowych reszt — odległości od każdego punktu do dopasowanej linii. To kryterium zwykłych najmniejszych kwadratów (OLS) ma rozwiązanie w postaci zamkniętej, co czyni je najczęściej stosowaną metodą statystyczną w nauce, ekonomii i uczeniu maszynowym.

R², współczynnik determinacji, mierzy, jaka część wariancji y jest wyjaśniana przez x: R² = 1,0 oznacza idealne dopasowanie; R² = 0 oznacza, że linia regresji nie radzi sobie lepiej niż zgadywanie średniej. r Pearsona to pierwiastek kwadratowy z R² (ze znakiem), mierzący kierunek i siłę zależności liniowej. Wpływowy punkt odstający może dramatycznie przesunąć linię regresji — żywa demonstracja jest warta tysiąca słów.

Cel OLS:
  Zminimalizować SSE = Σᵢ (yᵢ − ŷᵢ)²  gdzie ŷᵢ = a + b·xᵢ

Rozwiązanie w postaci zamkniętej (prosta regresja liniowa):
  b = Σ(xᵢ − x̄)(yᵢ − ȳ) / Σ(xᵢ − x̄)²
    = Cov(x, y) / Var(x)

  a = ȳ − b·x̄  (przecięcie przez centroid)

  Równoważnie w postaci macierzowej (regresja wieloraka):
    β̂ = (XᵀX)⁻¹ Xᵀy
    X = [1, x₁; 1, x₂; ...; 1, xₙ]  (macierz projektowa)

Miary dobroci dopasowania:
  SST = Σ(yᵢ − ȳ)²           (całkowita suma kwadratów)
  SSR = Σ(ŷᵢ − ȳ)²           (suma kwadratów regresji)
  SSE = Σ(yᵢ − ŷᵢ)²          (suma kwadratów reszt)
  R²  = 1 − SSE/SST = SSR/SST  ∈ [0, 1]

Korelacja Pearsona:
  r = Cov(x,y) / [σₓ · σᵧ]  ∈ [−1, 1]
  b = r · σᵧ / σₓ            (zależność nachylenia i korelacji)

Założenia Gaussa-Markowa dla estymatorów BLUE:
  1. Liniowość: y = Xβ + ε
  2. Ścisła egzogeniczność: E[ε | X] = 0
  3. Homoskedastyczność: Var(ε | X) = σ²I
  4. Brak doskonałej współliniowości: rank(X) = p

Symulacja regresji liniowej zamienia twoją mysz w edytor danych: kliknij, aby dodać punkt, przeciągnij, aby go przesunąć, i obserwuj, jak linia OLS, wynik R², r Pearsona i kwadraty reszt aktualizują się w czasie rzeczywistym. Przeciągnij jeden punkt odstający daleko od skupiska i zobacz, jak R² się zawala — dobitnie pokazując kruchość korelacji wobec obserwacji wpływowych. Pięć gotowych ustawień demonstruje idealne dopasowanie, zerową korelację, wzorce kwadratowe (których OLS nie potrafi uchwycić) i klasyczne scenariusze z punktami odstającymi.

Warstwa 5: błądzenia losowe z pamięcią — łańcuchy Markowa

Symulacja łańcucha Markowa

Łańcuch Markowa to proces stochastyczny, w którym przyszły stan zależy wyłącznie od stanu obecnego — a nie od tego, jak się do niego doszło. Ta „bezpamięciowa” własność Markowa jest zaskakująco potężna: leży u podstaw algorytmu PageRank Google, wyceny opcji finansowych, rozpoznawania mowy, modeli zwijania białek i metod całkowania Monte Carlo.

Mając macierz przejść P, gdzie P[i][j] to prawdopodobieństwo przejścia ze stanu i do stanu j, łańcuch Markowa ma unikalny rozkład stacjonarny π spełniający πP = π. Zaczynając od dowolnego stanu początkowego, powtarzane przejścia zbiegają do π — matematyczna gwarancja o głębokich implikacjach dla symulacji, optymalizacji i analizy sieci.

Macierz przejść P (n×n, wierszowo-stochastyczna):
  P[i][j] ≥ 0,   Σⱼ P[i][j] = 1  dla wszystkich i

Wektor prawdopodobieństwa stanu w kroku t:
  π(t) = π(0) · P^t

Rozkład stacjonarny (ustalony) π*:
  π* · P = π*    (lewy wektor własny dla wartości własnej 1)
  Σᵢ π*ᵢ = 1

Zbieżność metody potęgowej (dla łańcuchów ergodycznych):
  π(t) → π*  gdy t → ∞
  Tempo zależy od drugiej co do wielkości wartości własnej λ₂:
    czas zbieżności ∝ 1 / (1 − |λ₂|)

Gotowe ustawienie "pogoda" (w symulacji):
  Stany: Słonecznie, Pochmurno, Deszczowo
  P = [0,70  0,20  0,10]   (ze Słonecznie)
      [0,30  0,40  0,30]   (z Pochmurno)
      [0,20  0,30  0,50]   (z Deszczowo)
  Stacjonarny: π* ≈ [0,46, 0,28, 0,26]

Gotowe ustawienie PageRank:
  d = 0,85 (współczynnik tłumienia)
  P_PR = d · A_norm + (1−d) · (1/n) · 1·1ᵀ
  Stacjonarny π* = wyniki ważności stron

Symulacja łańcuchów Markowa rysuje animowane strzałki przejść między stanami, których szerokość jest proporcjonalna do prawdopodobieństw przejścia. Animowany wędrowiec skacze między stanami z regulowaną prędkością, podczas gdy wykresy słupkowe porównują empiryczną częstość odwiedzin z teoretycznym rozkładem stacjonarnym — zbiegając do π* na twoich oczach. Pięć gotowych ustawień obejmuje przejścia pogodowe, PageRank, ruinę gracza, genetykę i losowy łańcuch czterostanowy.

Warstwa 6: fizyka prawdopodobieństwa — rozkład Maxwella-Boltzmanna

Symulacja rozkładu Maxwella-Boltzmanna

Mechanika statystyczna to miejsce, w którym spotykają się fizyka i prawdopodobieństwo. Pojemnik z gazem zawiera ~10²³ cząsteczek, z których każda porusza się z inną prędkością. Opis deterministyczny jest niewykonalny, ale opis statystyczny jest dokładny: ułamek cząsteczek o prędkościach w przedziale [v, v+dv] podlega rozkładowi Maxwella-Boltzmanna, wyprowadzonemu przez Maxwella (1860) przy użyciu zaledwie trzech założeń: izotropii, niezależności składowych prędkości i zachowania energii.

Wynik pokazuje, że rozkład prędkości cząsteczek nie ma szczytu w zerze — większość cząsteczek porusza się znacznie szybciej niż prędkość termiczna √(k_BT/m). W wyższych temperaturach szczyt przesuwa się w prawo, a rozkład się poszerza; w niższych — zwęża się. Ma to bezpośrednie konsekwencje dla szybkości reakcji chemicznych (Arrhenius), ucieczki atmosfery z planet i działania laserów gazowych.

Funkcja gęstości prawdopodobieństwa:
  f(v) = 4π · (m / 2πk_BT)^(3/2) · v² · exp(−mv² / 2k_BT)

  gdzie:
    m   = masa cząsteczki [kg]
    k_B = stała Boltzmanna = 1,380649 × 10⁻²³ J/K
    T   = temperatura bezwzględna [K]

Prędkości charakterystyczne:
  Prędkość najbardziej prawdopodobna:  v_p  = √(2k_BT / m)
  Prędkość średnia:                    ⟨v⟩  = √(8k_BT / πm)  = v_p · √(4/π)
  Prędkość RMS:                        v_rms = √(3k_BT / m)  = v_p · √(3/2)

Kolejność: v_p < ⟨v⟩ < v_rms  (stosunek ≈ 1 : 1,128 : 1,225)

Normalizacja (sprawdzenie):
  ∫₀^∞ f(v) dv = 1  ✓

Skalowanie z temperaturą:
  v_p ∝ √T   — podwojenie temperatury zwiększa prędkość najbardziej prawdopodobną o √2
  v_p ∝ 1/√m — azot (28 u) porusza się ~3,4× wolniej niż wodór (2 u) w tej samej T

Ucieczka atmosferyczna (ucieczka Jeansa):
  Prędkość ucieczki z Ziemi: v_esc ≈ 11,2 km/s
  Wysokoenergetyczny ogon f(v) determinuje tempo utraty
  Wodór (m=2u) ucieka; azot (m=28u) nie ucieka

Symulacja Maxwella-Boltzmanna modeluje gaz cząstek o kształcie twardych kul, zderzających się sprężyście w dwuwymiarowym pudełku. Histogram prędkości cząstek w czasie rzeczywistym aktualizuje się przy każdym zderzeniu, zbiegając do krzywej MB w miarę trwania symulacji. Dostosuj suwak temperatury i obserwuj, jak reagują zarówno krzywa teoretyczna, jak i prędkości cząstek — pokazując namacalnie, że temperatura jest po prostu średnią energią kinetyczną chaotycznego ruchu.

Warstwa 7: resampling bez założeń — bootstrap

Symulacja resamplingu bootstrapowego

Statystyka klasyczna zakłada, że znasz rozkład populacji — a następnie wyprowadza wzory na przedziały ufności i wartości p przy tym założeniu. Ale co, jeśli nie możesz założyć normalności? Co, jeśli masz tylko 15 punktów danych i chcesz uzyskać przedział ufności dla mediany lub złożonej statystyki, takiej jak współczynnik Giniego?

Resampling bootstrapowy, wprowadzony przez Bradleya Efrona w 1979 roku, całkowicie omija założenia dotyczące rozkładu. Pomysł jest prosty: potraktuj swoją próbę jako zastępstwo populacji i wylosuj z niej tysiące nowych prób ze zwracaniem. Rozkład twojej statystyki w tych ponownych próbach przybliża jej prawdziwy rozkład próbkowania — bez potrzeby żadnych wzorów.

Algorytm (bootstrap nieparametryczny):
  1. Mając dane oryginalne: x = {x₁, x₂, ..., xₙ}
  2. Dla b = 1 do B (typowo B = 10 000):
       Wylosuj x*_b = {x*₁, ..., x*ₙ} ze zwracaniem z x
       Oblicz statystykę θ̂*_b = f(x*_b)  (np. średnia, mediana, R²)
  3. Posortuj {θ̂*_1, ..., θ̂*_B}
  4. 95% PU = [θ̂*_{0,025·B}, θ̂*_{0,975·B}]  (metoda percentylowa)

Przedział skorygowany o obciążenie i przyspieszenie (BCa):
  Koryguje obciążenie i skośność w rozkładzie bootstrapowym
  Preferowany zamiast prostej metody percentylowej dla małych prób

Dlaczego to działa (spostrzeżenie Efrona):
  Rozkład empiryczny F̂ₙ zbiega do prawdziwego F, gdy n → ∞
  Świat bootstrapu:  F̂ₙ → F̂ₙ  (znany dokładnie)
  Świat rzeczywisty: F̂ₙ → F    (nieznany cel)
  Obie zmienności próbkowania są asymptotycznie równe

Kluczowe zalety:
  • Działa dla dowolnej statystyki (średnia, mediana, korelacja, Gini itd.)
  • Brak założenia normalności
  • Obsługuje złożone schematy próbkowania
  • Przedziały ufności dla wydajności modeli ML (walidacja krzyżowa)

Symulacja resamplingu bootstrapowego wizualizuje każdą z ponad 1 000 ponownych prób jako mały punkt na wykresie punktowym, budując na żywo rozkład bootstrapowy. Suwak wprowadzający wartość odstającą pozwala wstrzyknąć jedną skrajną wartość do oryginalnej próby i obserwować, jak przedziały ufności rozszerzają się w czasie rzeczywistym — demonstrując graficznie, jak odporne są przedziały ufności bootstrap w porównaniu z przedziałami t, gdy założenia są naruszone.

Wszystkie siedem symulacji na pierwszy rzut oka

Głębszy związek

Te siedem symulacji to nie odizolowane tematy — tworzą one spójny łańcuch intelektualny. CTG wyjaśnia, dlaczego rozkłady normalne pojawiają się wszędzie, co uzasadnia test t i regresję OLS. Wnioskowanie bayesowskie dostarcza ram do aktualizacji dowolnego oszacowania prawdopodobieństwa, w tym parametrów modelu regresji. Łańcuchy Markowa leżą u podstaw MCMC (Markov Chain Monte Carlo), techniki używanej do próbkowania ze złożonych rozkładów a posteriori bayesowskich. Paradoks urodzin jest szczególnym przypadkiem prawdopodobieństwa kolizji w tablicach mieszających i wydobywaniu kryptowalut. A rozkład Maxwella-Boltzmanna to CTG zastosowane do 10²³ cząsteczek gazu — ta sama matematyka, która wyjaśnia średnie z próby, rządzi także termodynamiką.

Statystyka to nie zbiór wzorów do zapamiętania, lecz sposób myślenia o niepewności, dowodach i wnioskowaniu. Każda symulacja tutaj buduje jeden pokój w tej mentalnej architekturze — a razem tworzą dom dla myślenia ilościowego.