Statystyka w ujęciu Bayesowskim

Z Skrypty dla studentów Ekonofizyki UPGOW

(Różnice między wersjami)

Aktualna wersja na dzień 20:12, 7 cze 2011

Spis treści

1 Statystyka jako wnioskowanie w warunkach niepewności
2 Problemy fundamentalne: jak rozumieć prawdopodobieństwo?
- 2.1 Podstawowe reguły operacyjne rachunku prawdopodobieństwa
3 Wnioskowanie Bayesowskie
4 Bayesowska teoria estymacji parametrów modelu
- 4.1 Marginalizacja
5 Prawdopodobieństwa aprioryczne
6 Wiarygodność Bayesowska
7 Metodyka uzyskiwania rozkładów a posteriori
- 7.1 Metody próbkowania
- 7.2 Metoda Łańcuchów Markowa.
8 Selekcja modeli
9 Literatura

Statystyka jako wnioskowanie w warunkach niepewności

We współczesnym ujęciu statystyka rozumiana jest jako dyscyplina umożliwiająca wnioskowanie w warunkach niepewności.

Rola wnioskowania dedukcyjnego w nauce jest ugruntowana od starożytności. Opiera się ono na logice, która dostarcza reguł wnioskowania. Wywodzą się one zazwyczaj z tautologii. W matematyce budując jakiś określony dział, np. geometrię czy algebrę liniową, punktem wyjścia są aksjomaty – zdania przyjmowane jako oczywiste. Aksjomaty wprowadzają podstawowe elementy (pojęcia) danej teorii, definiują jej zakres, a reguły wnioskowania służą do formułowania twierdzeń (wyprowadzania ich z aksjomatów). Z kolei w naukach empirycznych zajmujących się badaniem jakiegoś obszaru rzeczywistości, odkrywana tam struktura ujmowana jest pod postacią praw, których słuszności staramy się dociec formułując testowalne predykcje (przewidywania). Niestety, testy eksperymentalne przewidywań teoretycznych nie dostarczają prostych zero-jedynkowych odpowiedzi „prawda” lub „fałsz”. Nasz stan wiedzy zawsze jest niekompletny, zawsze możliwy jest do pomyślenia bardziej wyrafinowany eksperyment, a przyrządy pomiarowe zawsze mają skończoną dokładność. Wnioskowanie o słuszności teorii czy modelu opisującego jakiś obszar rzeczywistości w warunkach niepełnej wiedzy (tj. częściowej ignorancji, niepewności) jest przedmiotem statystyki.

Statystyka, w swej koncepcji, oparta jest na rachunku prawdopodobieństwa, w teorii prawdopodobieństwa zaś – od momentu jej powstania – istnieją dwie szkoły pojmowania, czym operacyjnie (w praktyce) jest prawdopodobieństwo. Pierwsza z nich tzw. częstościowa – pochodząca od Abrahama de Moivre’a – twierdzi iż prawdopodobieństwo wystąpienia zdarzenia losowego w praktyce jest reprezentowane przez częstość występowania tego zdarzenia w bardzo dużej liczbie identycznych prób. Ale co począć, gdy rozważane zdarzenie, zjawisko, sytuacja są z natury niepowtarzalne? A niepowtarzalny jest świat w którym żyjemy, lub czasy w których żyjemy, każdy człowiek jest niepowtarzalny i w wielu przypadkach zdroworozsądkowego używania terminu „prawdopodobieństwo” trudno byłoby uzasadnić jego poprawność na gruncie częstościowej definicji. Czy, na przykład stwierdzenie typu: „prawdopodobieństwo awarii współczesnej elektrowni jądrowej wynosi jeden na milion” ma sens? W rozumieniu częstościowym – nie ! W historii ludzkości nie zbudowano jeszcze miliona elektrowni jądrowych, a spośród tych które zbudowano kilka uległo mniej lub bardziej poważnym awariom. Jednak czujemy sensowność przytoczonych słów jako określenie miary niezawodności technologii konstrukcji elektrowni jądrowych.

Tu pojawia się druga szkoła pojmowania prawdopodobieństwa, której autorem jest Thomas Bayes. Zgodnie z tą szkołą, prawdopodobieństwo (bezwarunkowe, tzw. a priori) wystąpienia zdarzenia losowego jest niczym innym jak miarą racjonalnego przekonania, że dane zdarzenie wystąpi. Chcąc zmienić (zmodyfikować, wzbogacić) nasze przekonania wykonujemy eksperymenty (obserwacje) dotyczące interesującego nas zdarzenia. Wyniki badań przekształcają prawdopodobieństwo a priori (wstępne oczekiwania) w tzw. prawdopodobieństwo a posteriori (prawdopodobieństwo wynikowe, miara racjonalnego oczekiwania wystąpienia zdarzenia po uzyskaniu wyników badań). Generalnie w taki sposób odbywa się nasze poznawanie świata. Zwolennikami filozofii Bayesa byli tacy wielcy matematycy jak P.S. Laplace czy H. Poincare lub wybitny ekonomista John Keynes, ale dopiero obecnie rozwój technik numerycznych umożliwił przekształcenie się Bayesowskiego podejścia w poważny, awangardowy nurt współczesnej statystyki, umożliwiający stawianie i rozwiązywanie problemów niedostępnych dla klasycznej statystyki częstościowej.

Problemy fundamentalne: jak rozumieć prawdopodobieństwo?

Czym jest prawdopodobieństwo? W klasycznym ujęciu tzw. „częstościowym” prawdopodobieństwo zdarzenia A jest częstością występowania tego zdarzenia w długiej serii identycznych prób, czyli stosunkiem

$$P(A) = \frac{n_A}{n}$$

gdzie $n_A$ reprezentuje liczbę zdarzeń $A$ w $n$ – próbach, gdzie „$n$ jest duże” tzn. „$n$ dąży do nieskończoności”. Definicja taka jest niesatysfakcjonująca z kilku powodów. Po pierwsze, co to znaczy że próby są identyczne? To znaczy, że szansa (czyli potocznie rozumiane „prawdopodobieństwo”) wystąpienia zdarzenia $A$ w każdej z prób jest taka sama – zatem definiujemy pojęcie prawdopodobieństwa w istocie odwołując się do niego samego. W logice jest to fundamentalny błąd zwany petitio principium. Po drugie definicja taka nie ma zastosowania do zdarzeń niepowtarzalnych. Na przykład stwierdzenie historyka, typu: „ na podstawie zapisków historycznych, z dużym prawdopodobieństwem można stwierdzić iż podczas koronacji Władysława Jagiełły nie padał deszcz” w myśl częstościowej koncepcji prawdopodobieństwa jest kompletnie pozbawiona sensu. Z drugiej jednak strony nie widać w tej wypowiedzi absurdu, tym bardziej jeśli wynika ona z analizy konkretnych danych historycznych i znając ów szerszy kontekst można taką tezę uznać nie tylko za sensowną, ale wręcz za prawdziwą. Po trzecie, definicja częstościowa dotyczy pewnego zachowania granicznego, i to w dodatku „granicznego” w specyficzny, słabo określony sposób. Mianowicie nie jest to granica w zwykłym matematycznym sensie znanym z analizy (wręcz można pokazać, że taka „zwykła” granica nie istnieje) i nie jest wcale oczywistym kiedy $n$ będzie dostatecznie duże aby mierzona częstość $n_A/n$ faktycznie reprezentowała prawdopodobieństwo. Znacznie szerszą, nieuwikłaną w wyżej wymienione problemy, perspektywą jest Bayesowskie rozumienie prawdopodobieństwa jako miary racjonalnego zaufania w prawdziwość danej tezy, zaufania uwarunkowanego posiadaną informacją. W tym kontekście np. przestajemy mieć problem z wypowiedzią o stanie pogody podczas koronacji Jagiełły.

Zalety Bayesowskiego rozumienia prawdopodobieństwa:

1. Ma zastosowanie zarówno do zjawisk powtarzalnych jak i niepowtarzalnych, tzn. pytanie: „Z jakim prawdopodobieństwem będzie jutro padał deszcz?” jest w pełni sensownym pytaniem probabilistycznym.

2. Prawdopodobieństwo staje się teraz terminem opisu niepewności niezależnie od jej pochodzenia. Na równych prawach traktowane są tzw. „błędy statystyczne” tzn. niepewność wynikająca z błędów pomiarowych czyli ze skończonej dokładności przyrządów oraz „błędy systematyczne” czyli niepewność związana z naszą niewiedzą odnośnie skądinąd czysto deterministycznych procesów.

3. Losowość, czy przypadkowość zjawisk rozumiana jest tu jako wyraz niepełnej informacji jaką posiadamy.

4. W odróżnieniu od podejścia częstościowego, Bayesowskie rozumienie prawdopodobieństwa nie odnosi się do granicznych własności estymatorów czy statystyk. Ma ono zastosowanie zarówno w przypadkach gdy podejście częstościowe ma zastosowanie, jak również w przypadkach gdzie nie ma ono sensu.

5. Bayesowskie podejście automatycznie radzi sobie z nieistotnymi parametrami (ang. nuisance parameters) w modelach statystycznych

6. W podejściu Bayesowskim istotna jest informacja aprioryczna, np. ważne przesłanki że masa, czy objętość są zawsze dodatnie. W rzeczy samej, ignorowanie lub pomijanie takich apriorycznych faktów prowadzić może do fałszywych wniosków.

7. Statystyka w ujęciu Bayesowskim zawsze (tj. w obliczeniach i interpretacji) odnosi się tylko do danych, które faktycznie zostały otrzymane, podczas gdy statystyka częstościowa odnosi się (w swej interpretacji) do rozkładu wyników, które są potencjalnie możliwe, lecz de facto nie zostały zaobserwowane.

Różnicę pomiędzy podejściem częstościowym i Bayesowskim można syntetycznie opisać następująco:

W podejściu częstościowym, gdy piszemy prawdopodobieństwo $p(X)$, to $X$ oznacza zmienną losową czyli taką, która może przyjmować różne wartości w nieskończonym zespole (wyimaginowanym ansamblu) identycznych eksperymentów. Eksperymenty te są potencjalne, wirtualne, wyimaginowane. $X$ ma w tym ansamblu rozkład prawdopodobieństwa $p$, czyli częstość znalezienia wartości $x$ w przedziale $[x, x+dx]$ wynosi $p(x)dx$.

W podejściu Bayesowskim $ X$ ma pewną ustaloną, konkretną – lecz obarczoną niepewnością – wartość, natomiast $p(x)dx$ opisuje (warunkowy) rozkład miary racjonalnego zaufania co do wartości $x$. Rozkład ten jest warunkowy, gdyż uwarunkowany jest posiadaną informacją – ogólnym kontekstem problemu oraz danymi z eksperymentu. W nauce dostępna informacja jest zawsze niepełna, więc nasza wiedza o prawach natury jest probabilistyczna.

We wnioskowaniu Bayesowskim, funkcja rozkładu prawdopodobieństwa jest sposobem „zakodowania” niepewności odnośnie pewnych parametrów modelu lub konkurujących teorii czy modeli, przy danym stanie wiedzy reprezentowanym przez informację $I$.

Podstawowe reguły operacyjne rachunku prawdopodobieństwa

Rachunek prawdopodobieństwa można operacyjnie sprowadzić do następujących reguł:

0. Każdemu zdarzeniu losowemu $A$ (w rozumieniu Bayesowskim: zdarzeniu lub wypowiedzi obarczonymi niepewnością) można przypisać liczbę $ p(A|I) ;\,\, 0 \leq p(A|I) \leq 1 $ zwaną prawdopodobieństwem

1. Reguła sumy: dla zdarzenia A i jego dopełnienia zachodzi $ p(A|I) + p(\sim A|I) = 1 $

2. Reguła iloczynu: prawdopodobieństwo zajścia łącznego zdarzeń wynosi: \[ p(A,B|I) = p(A|B,I) P(B|I) = p(B|A,I) p(A|I) \]

Powyższe reguły implikują tzw. regułę marginalizacji \[ p(B|I) = \sum_A p(A,B|I) \]

Mówi ona, że jeśli znamy prawdopodobieństwo łączne zdarzeń $A$ i $B$, to prawdopodobieństwo samego zdarzenia $B$ – niezależnie od zdarzenia $A$ uzyskamy sumując prawdopodobieństwo łączne po wszystkich możliwych wartościach $A$.

Z reguły iloczynu wynika Twierdzenie Bayesa:

\[ p(A|B,I) = \frac{p(B|A,I) p(A|I)}{p(B|I)} \]

Wnioskowanie Bayesowskie

Przy odpowiedniej interpretacji, twierdzenie Bayesa pokazuje w jaki sposób nowe dane eksperymentalne/obserwacyjne modyfikują wstępne oczekiwania tj. prawdopodobieństwo a priori.

Oznaczając: $H_i$ = zdanie deklarujące słuszność i-tej hipotezy $I$ = zdanie reprezentujące informacje a priori $D$ = zdanie reprezentujące dane (konkretne uzyskane dane!) przepisujemy wzór Bayesa w postaci

\[ p(H_i|D,I) = \frac{p(D|H_i,I) p(H_i|I)}{p(D|I)} \]

gdzie: $p(D|H_i, I)$ = prawdopodobieństwo otrzymania danych $D$, pod warunkiem słuszności hipotezy $H_i$ (tzw. Bayesowska funkcja wiarygodności )

$p(H_i|I)$ = prawdopodobieństwo hipotezy a priori, tzn. co wiemy o hipotezie $H_i$ zanim zobaczymy dane

$p(H_i|D, I)$ = prawdopodobieństwo hipotezy a posteriori, w świetle uzyskanych danych $D$, czyli stan naszej wiedzy o $H_i$ po otrzymaniu danych

$p(D|I)$ = czynnik normalizacyjny – niezależny od hipotez $H_i$ , tzw. wiarygodność globalna, czy „ewidencja” (ang. evidence)

\[ p(D|I) = \sum_i p(H_i|I) p(D|H_i,I) \]

Istotnie powyższy wzór zapewnia normalizację prawdopodobieństw a posteriori:

\[ \sum_i p(H_i|D,I) = 1 \]

Należy podkreślić, że Bayesowskie podejście nie dyskredytuje podejścia częstościowego, argumenty częstościowe mogą być pomocne w ustalaniu wiarygodności Bayesowskiej lub prawdopodobieństw a apriorycznych (wyjściowych).

Bayesowskie podejście do wnioskowania statystycznego streścić można w następujący sposób:

• zawsze posiadamy a priori pewną wiedzę (informację $I$) odnośnie badanego zjawiska i na jej podstawie formułujemy hipotezy $H_i$. Wiedzę tą wyraża prawdopodobieństwo a priori $p(H_i|I)$ – prawdopodobieństwo jest tu miarą naszych wstępnych oczekiwań co do słuszności $H_i$.

• w celu lepszego poznania badanego zjawiska, wykonuje się odpowiednio zaprojektowane badania odzwierciedlające nasze dotychczasowe jego zrozumienie i wynikające z niego przewidywania.

• czy i w jakim stopniu nowe badanie zwiększy nasze poznanie, jest zdeterminowane przez tzw. bayesowską funkcję wiarygodności $p(D| H_i)$.

• po wykonaniu badań ich wyniki (dane $D$) modyfikują nasze oczekiwania wiodąc do tzw. prawdopodobieństwa a posteriori $p(H_i|D)$, które wzmacnia lub osłabia nasze oczekiwania (w zależności czy jest większe czy mniejsze od a priorycznego).

Prawdopodobieństwo warunkowe $p(H_i|D)$ mierzy stopień poparcia dawanego hipotezie $H_i$ przez dane doświadczalne $D$, tzn. mierzy w jakim stopniu dane $D$ czynią rozróżnienie pomiędzy hipotezą $H_i$ i innymi alternatywnymi hipotezami. Wzór Bayesa natomiast pokazuje w jaki sposób prawdopodobieństwo a priori $p(H_i)$ wiąże się z prawdopodobieństwem a posteriori $p(H_i |D)$:

$ p(H_i |D) = p(H_i) p(D | H_i) / p(D) $

Rysunek 3.1. Zbiorczy wykres rozkładu a priori, wiarygodności bayesowskiej oraz rozkładu a posteriori jako ilustracja wnioskowania bayesowskiego.

Przykłady wnioskowania Bayesowskiego z dyskretną przestrzenią hipotez

Przykład 3.1.1

W problemie rzutu monetą, stwierdzenie iż prawdopodobieństwo otrzymania orła (O) wynosi 1/2 jest w istocie wiarygodnością Bayesowską: $p(O | moneta\; rzetelna) = 1/2$ tzn. wyrażeniem typu p(D|M) (dane: $D$ = O, model: $M$ = rzetelna moneta).

Przy innym rozważanym modelu – nazwijmy go M’ (tj. monety nierzetelnej – takiej która preferuje jeden z wyników np. O ) mielibyśmy: p(O|M’) = $ \pi $ oraz p(R | M’) = 1 – $ \pi $ Oczywiście wartość $ \pi $ byłaby tu identyfikatorem modelu M’ .

Przypuśćmy, że ktoś rzucił monetą 4 razy uzyskując za każdym razem orła $ D={4 \times O} $. Rozumując częstościowo: Co możemy powiedzieć o rzetelności monety? Z jakim prawdopodobieństwem eksperyment ten dyskredytuje model $M$ monety rzetelnej? Czy można wskazać model $ M’$ tzn. identyfikujące go prawdopodobieństwo $ \pi $ ?

W podejściu Bayesowskim rozważymy ten przykład później.

Przykład 3.1.2

Załóżmy, że mamy urnę w której jest 6 kul, które mogą być białe (B) lub czarne (Cz) wykluczamy sytuację jednakowego koloru. Nie wiemy natomiast jaki jest skład zawartości urny. Poznanie składu tj. proporcji kul Cz i B jest naszym celem badawczym. Potencjalne możliwości tzn. możliwe do pomyślenia modele tej sytuacji są następujące:

$ M1 = 5 \times B + 1 \times Cz $

$ M2 = 4 \times B + 2 \times Cz $

$ M3 = 3 \times B + 3 \times Cz $

$ M4 = 2 \times B + 4 \times Cz $

$ M5 = 1 \times B + 5 \times Cz $

Każdy z modeli jest równoprawny, więc przypisujemy im takie samo prawdopodobieństwo aprioryczne:

$ P(M1) = P(M2) = … = P(M5) = 1/5 = 0.20 $

Aby przekonać się, który model jest słuszny, losujemy 2 kule i powiedzmy że wynikiem jest $D = {2 x B}$. Jaki jest stopień poparcia dany poszczególnym modelom przez ten wynik?

Zdroworozsądkowo czujemy że preferowany jest model $M1$, a model $M5$ jest przez ten wynik sfalsyfikowany. Dojdziemy do tego systematycznie poprzez zastosowanie wzoru Bayesa:

$ p(M_i|D) = \frac{ p(D|M_i)}{ p(D)}\,p(M_i) $

Potrzebne nam będą wiarygodności Bayesowskie: $p(2 \;B | M1) = 5/6 \times 4/5 = 0.667 $

$p(2 \; B | M2) = 4/6 \times 3/5 = 0.400 $

$p(2 \; B | M3) = 3/6 \times 2/5 = 0.200$

$p(2 \; B | M4) = 2/6 \times 1/5 = 0.067$

$p(2 \; B | M5) = 1/6 \times 0/5 = 0. $

Widzimy tu formalnie, że model $M5$ nie przewiduje wylosowania 2 białych kul, czyli w świetle uzyskanych danych, jego wiarygodność Bayesowska, a co za tym idzie także i prawdopodobieństwo a posteriori są zerowe.

Czynnik normalizacyjny:

$ p(D) = p( 2\;B) = \sum_i p(2\;B|M_i) p(M_i) = 0.667 \cdot 0.20 + 0.40 \cdot 0.20 + 0.067 \cdot 0.20 = 0.2668 $

Prawdopodobieństwa a posteriori:

$ p(M_1|2\;B) = \frac{0.667 \cdot 0.20}{0.2668} = 0.502 $

$ p(M_2|2\;B) = \frac{0.400 \cdot 0.20}{0.2668} = 0.300 $

$ p(M_3|2\;B) = \frac{0.20 \cdot 0.20}{0.2668} = 0.150 $

$ p(M_4|2\;B) = \frac{0.067 \cdot 0.20}{0.2668} = 0.050 $

$ p(M_5|2\;B) = \frac{0. \cdot 0.20}{0.2668} = 0. $

Widzimy zatem jak rozłożyła się alokacja poparcia dawanego różnym modelom przez uzyskane dane.

Procedury diagnostyczne w terminach Bayesowskich

W praktyce nierzadko mamy do czynienia z testami diagnostycznymi, tj. procedurami których wynik interpretowany jest jako „dodatni” $T+ $ lub „ujemny” $T- $ . Wynikiem procedury diagnostycznej może być liczba, sygnał lub zmiana barwy wskaźnika, ważne jest aby istniała jednoznaczna interpretacja wyniku $T+$ lub $T-$ .

Procedury diagnostyczne stosowane są w celu rozróżnienia pomiędzy 2 alternatywnymi (komplementarnymi w sensie logicznym) hipotezami $H+$ oraz $H-$.

Typowym przykładem testów diagnostycznych są testy stosowane w medycynie, np. test nosicielstwa wirusa HCV, próba wysiłkowa w diagnostyce choroby wieńcowej itp. Wówczas hipotezy $H+$ i $H-$ mają typowo znaczenie: $H+ $ = choroba / nosicielstwo wirusa itp. $H-$ = brak choroby / brak nosicielstwa wirusa itp.

Procedura diagnostyczna scharakteryzowana jest przez zespół prawdopodobieństw warunkowych przytoczonych w tabeli poniżej:

*Prawdopodobieństwa warunkowe charakteryzujące procedurę* diagnostyczną
	Choroba
Test		H +	H -
T +	$ p(T+ \| H+) $ czułość testu; częstość wyników prawdziwie dodatnich	$ p(T+ \| H-) $ 1 - swoistość; częstość wyników fałszywie dodatnich
T -	$ p(T - \| H+) $ 1 - czułość; częstość wyników fałszywie ujemnych	$ p(T -\| H -) $ swoistość testu; częstość wyników prawdziwie ujemnych

Jak widać z czterech przytoczonych wielkości tylko dwie są niezależne, zatem test diagnostyczny scharakteryzowany jest w pełni przez swoistość i czułość. Lub też alternatywnie poprzez częstość wyników fałszywie ujemnych oraz fałszywie dodatnich. Ten drugi alternatywny sposób charakteryzowania testu jest typowy w obszarze testowania hipotez statystycznych.

Czułość i swoistość są wewnętrznymi cechami procedury diagnostycznej, natomiast w praktyce interesuje nas interpretacja wyników dodatnich $T+$ oraz ujemnych $T-$, czyli wartość dodana wyników $T+$ oraz $T-$.

Owa wartość dodana reprezentowana jest łącznie przez dwa prawdopodobieństwa a posteriori, odpowiednio: $p(H+ | T+)$ oraz $p(H- | T-)$. W biostatystyce noszą one miano, odpowiednio: wartości predyktywnej wyników dodatnich oraz wartości predyktywnej wyników ujemnych.

Są one dane wzorami:

\[ p(H+|T+) = \frac{p(T+|H+) p(H+)}{p(T+|H+) p(H+)+ p(T+|H-) p(H-)} \]

\[ p(H-|T-) = \frac{p(T-|H-) p(H-)}{p(T-|H-) p(H-)+ p(T-|H+) p(H+)} \]

Do zastosowania powyższych wzorów, oprócz czułości i swoistości testu należy znać prawdopodobieństwo a priori $p(H+)$ o którym wiedzę czerpiemy z częstości występowania danej choroby (H+) czy nosicielstwa wirusa w populacji generalnej.

Uwaga

Weryfikacja hipotez statystycznych w klasycznej statystyce odbywa się w pełnej analogii do testów diagnostycznych. Podana wyżej tabela zachowuje ważność, przy zmianie oznaczeń $(H+,H-)$ na, odpowiednio $(H_0, H_1)$. Klasycznie $H_0$ jest symbolem tzw. hipotezy zerowej, $H_1$ tzw. hipotezy alternatywnej. Tu także procedura wnioskowania tj. test statystyczny scharakteryzowana jest przez dwie liczby (prawdopodobieństwa warunkowe):

1. Prawdopodobieństwo wyników fałszywie dodatnich; nazywa się ono historycznie prawdopodobieństwem błędu I rodzaju (oznaczane $\alpha$)

2. Prawdopodobieństwo wyników fałszywie ujemnych; historycznie nazywane prawdopodobieństwem błędu II rodzaju (oznaczane $\beta$)

Ich pochodne wielkości to: czułość – zwana mocą testu $(1 – \beta)$ oraz swoistość $(1 – \alpha)$ Różnicą jest jedynie to, że interpretacja testu jako dodatniego czy ujemnego związana jest z założonym prawdopodobieństwem błędu I rodzaju, zwanym także poziomem istotności testu.

Przykład Bayesowskiego wnioskowania w procesie diagnostycznym

Do lekarza zgłasza się mężczyzna uskarżający się na nadmierną męczliwość, w wywiadzie podaje chorobę wrzodową oraz epizody kamicy nerkowej. Lekarz rozważa jako jedną z możliwości nadczynność przytarczyc, pomimo iż nie jest to częsta choroba. Jej częstość występowania w populacji generalnej, na podstawie literatury, może być oceniona na 2%. Jest to zarazem nasze prawdopodobieństwo a priori $p(H+)$. Zleca oznaczenie poziomu wapnia w surowicy (aby wykluczyć wstępne rozpoznanie). Wynik badania jest jednak pozytywny tzn. poziom wapnia jest podniesiony ($T+$).

Jak wynik badania zmodyfikował prawdopodobieństwo a priori tzn. stopień poparcia dawanego wyjściowej diagnozie? Załóżmy, że test podniesionego poziomu wapnia ma czułość 90% (czyli p(T+|H+) = 0.90 ) i swoistość 95% (czyli p(T-|H-) = 0.95).

Wnioskujemy stąd prawdopodobieństwo wyników fałszywie dodatnich $p(T+ | H- ) = 0.05$ i posiadamy wszelkie dane aby zastosować wzór Bayesa

$$ p(H+|T+) = \frac{p(T+|H+) p(H+)}{p(T+|H+) p(H+)+ p(T+|H-) p(H-)} = \frac{ 0.90 \times 0.02} {0.90 \times 0.02 + 0.05 \times 0.98} = 0.269 = 27 \% $$

Widzimy jak wzrósł stopień poparcia dla rozpoznania pierwotnego po uzyskaniu dodatniego wyniku testu. Niskie prawdopodobieństwo a priori 2% przeobraziło się w znacznie (dziesięciokrotnie) wyższe prawdopodobieństwo a posteriori 27% .

Kolejny krok w procesie diagnostycznym

Mimo, że prawdopodobieństwo a posteriori jest znacznie większe od początkowego – a apriorycznego, daleko mu jeszcze do pewności (100%). Zatem diagnosta zleca kolejne, badanie tym razem już swoiste dla wstępnej diagnozy. Badaniem tym jest oznaczenie poziomu parathormonu (PTH) – hormonu produkowanego przez przytarczyce. Gdyby prawdziwa była diagnoza nadczynności tego gruczołu, spodziewamy się podniesionego poziomu PTH (T+), w przeciwnym przypadku uznajemy test za ujemny (T-).

Załóżmy, że test PTH ma swoistość 98% (tzn. $p(T- | H-) = 0.98$) i czułość 95% (tzn.$p(T+ | H+) = 0.95$), prawdopodobieństwo wyników fałszywie dodatnich wynosi więc $p(T+ | H-) = 0.02 $ natomiast jako prawdopodobieństwo a aprioryczne przyjmujemy, zgodnie z Bayesowską filozofią wnioskowania, poprzednie prawdopodobieństwo końcowe $p(H+) = 27 \% = 0.27$.

W jakim stopniu dodatni wynik testu zmodyfikuje stopień poparcia dawanego diagnozie? Odpowiedź tkwi w formule Bayesa

$ p(H+|T+) = \frac{p(T+|H+) p(H+)}{p(T+|H+) p(H+)+ p(T+|H-) p(H-)} = \frac{ 0.95 \times 0.269} {0.95 \times 0.269 + 0.02 \times 0.731} = 0.9449 = 94.5 \% $

Silny wzrost prawdopodobieństwa a posteriori obrazuje rozstrzygające znaczenie wykrycia podwyższonego poziomu parathormonu dla ostatecznej diagnozy.

Przykład - dlaczego przy dwóch sprzecznych wynikach testu trzeci jest rozstrzygający

Kobieta posiadająca brata chorego na hemofilię spodziewa się dziecka. Przypomnijmy że hemofilia jest poważnym zaburzeniem krzepliwości krwi uwarunkowanym genetycznie tzn. sprzężonym z chromosomem X. Zatem może potencjalnie być ona nosicielką hemofilii ($H+$). Szanse tego (czyli prawdopodobieństwo a priori) wynoszą $p(H+) = 0.5$ (albo jest albo nie jest nosicielką). Aby się o tym przekonać wykonuje ona test DNA, który jest swoisty na poziomie 90% i czuły na poziomie 80%.

Załóżmy, że wynik testu jest ujemny ($ T- $). Jak wpływa to na stopień jej oczekiwań, że może być nosicielką. Innymi słowy, jakie jest prawdopodobieństwo a posteriori p(H+|T-) ?

Odpowiedź daje wzór Bayesa

$ p(H+|T-) = \frac{p(T-|H+) p(H+)}{p(T-|H+) p(H+)+ p(T-|H-) p(H-)} = \frac{ 0.1 \times 0.5} {0.1 \times 0.5 + 0.8 \times 0.5} = 0.11 = 11 \% $

Czyli prawdopodobieństwo nosicielstwa – innymi słowy ryzyko zmalało prawie pięciokrotnie. Mimo, że wynik ten znacznie ją uspokoił, postanowiła wykonać kolejny test i złóżmy, że teraz wypadł pozytywnie $T+$. Jak zmieniło się teraz prawdopodobieństwo a posteriori, czyli jakie jest obecnie ryzyko bycia nosicielem? Znów stosujemy wzór Bayes, pamiętając jednak, że rolę prawdopodobieństwa a priori P(H+) pełni teraz wynik poprzedni, czyli poprzednie prawdopodobieństwo wynikowe (a posteriori). Jest to w pełni zgodne z filozofią wnioskowania Bayesowskiego: na wstępne subiektywne oczekiwania tzn. prawdopodobieństwa a aprioryczne składa się posiadana wyjściowo wiedza. Zauważmy, że poprzedni poziom oczekiwań wynosił 0.5 nie tylko dlatego, że nacechowany był niewiedzą o stanie faktycznym ale także z powodu już istniejącej wiedzy na temat natury hemofilii.

Zatem stosujemy po raz kolejny wzór Bayesa, obecnie w postaci

$ p(H+|T+) = \frac{p(T+|H+) p(H+)}{p(T+|H+) p(H+)+ p(T+|H-) p(H-)} = \frac{ 0.90 \times 0.11} {0.90 \times 0.11 + 0.20 \times 0.89} = 0.357 = 35.7 \% $

Sens wyniku jest zrozumiały: oczywiście ryzyko wzrosło (3.5 krotnie) ale nadal jest mniejsze od wstępnych oczekiwań, bo przecież wśród wykonanych badań mieliśmy wynik ujemny!

Zdrowy rozsądek podpowiada, że dwa sprzeczne wyniki (T+ oraz T-) są niekonkluzywne, więc w takich sytuacjach wykonuje się trzecie – rozstrzygające badanie. Powiedzmy, że wypada ono negatywnie (T-). Jaki jest teraz poziom ryzyka bycia nosicielem?

Stosujemy po raz kolejny wzór Bayesa, w celu uaktualnienia racjonalnego zaufania

$ p(H+|T-) = \frac{p(T-|H+) p(H+)}{p(T-|H+) p(H+)+ p(T-|H-) p(H-)} = \frac{ 0.1 \times 0.357} {0.1 \times 0.357 + 0.8 \times 0.643} = 0.065 = 6.5 \% $

Teraz widać, że wstępne a aprioryczne prawdopodobieństwo (ryzyko) uległo w sumie ok. 10 krotnej redukcji, co powinno w pełni uspokoić pacjentkę.

Jako ćwiczenie dobrze jest przeanalizować w analogiczny sposób inne konfiguracje dwóch negatywnych i jednego pozytywnego wyniku tj. (T+,T-,T-) oraz (T-, T- , T+).

Przykład ten ilustruje także, następujące fakty. Po pierwsze zdroworozsądkową regułę iż przy dwóch sprzecznych wynikach testu trzeci jest decydujący. Po drugi pokazuje, że stopień zaufania w wynik rozstrzygający jest różny od naiwnego wniosku 2/3. Po trzecie, że ów stopień zaufania zależy od jakości testu, tzn. jego czułości i swoistości, które to wielkości były uwzględniane w naszych obliczeniach jako $p(T+|H+)$ oraz $p(T-|H-)$.

Bayesowska teoria estymacji parametrów modelu

Często mamy do czynienia z problemem, w którym pewien model jest z założenia słuszny, a przestrzeń hipotez dotyczy wartości parametrów tego modelu. Na przykład w problemie regresji liniowej zakładamy liniową zależność pomiędzy zmienną objaśnianą $Y$ i zmienną objaśniającą $X$. Zatem $ M: \;\; Y = a X + b $ jest modelem, jego parametrami są tu $a$ i $ b$ czyli współczynnik kierunkowy (ang. slope) oraz wyraz wolny (ang. intercept). Jeśli parametry modelu nie zostaną podane explicite, będziemy je oznaczać ogólnie symbolem $ \theta $. Punktem wyjścia jest dla nas, jak zwykle wzór Bayesa

\[ p(H_i|D,I) = \frac{p(D|H_i,I)}{p(D|I)} p(H_i|I) \]

Z warunkiem normalizacyjnym określającym wyrażenie p(D|I)

\[p(D|I) = \sum_i p(D|H_i,I) p(H_i|I) \]

Z wykorzystaniem wzoru Bayesa w przypadku zespołu dyskretnych hipotez, którym przypisać można dyskretne liczbowe wartości prawdopodobieństwa rozumianego jako miarę racjonalnego poparcia, spotkaliśmy się już wcześniej w przykładach. Dział klasycznej statystyki znany jako teoria estymacji parametrów w języku Bayesowskim stanowi uogólnienie wzoru Bayesa na przypadek ciągłej przestrzeni hipotez. Hipotezą jest wówczas wartość estymowanego parametru. Wówczas prawdopodobieństwo hipotezy $H_0 $ czyli $p(H_0 |D,I)$ rozumiane jest podobnie jak gęstość rozkładu prawdopodobieństwa

\[ p(H_0|D) = \lim_{\Delta h \rightarrow 0} \frac{ p(h \leq H_0 \leq h + \Delta h) }{\Delta h} \]

Niech $M$ oznacza aprioryczną informację $I$, że model parametryzowany parametrem $\theta $ jest słuszny, wówczas $ p(\theta|M)$ jest aprioryczną gęstością prawdopodobieństwa, tzn. $p(\theta|M)\,d \theta $ jest prawdopodobieństwem iż prawdziwa wartość $\theta $ leży w przedziale $[\theta ,\; \theta + d \theta]$. W tym przypadku globalna wiarygodność (zwana w tym kontekście wiarygodnością modelu) jest uogólnieniem wzoru opisującego czynnik normalizacyjny na przypadek ciągłej przestrzeni hipotez

\[ p(D|M) = \int p(D| \theta ,M) p(\theta |M)\;d \theta = L(M) \]

Globalna wiarygodność modelu jest więc równa średniej ważonej wiarygodności jego parametrów. W podejściu Bayesowskim estymacja parametrów nie polega na konstrukcji estymatorów, tzn. znalezieniu punktu w przestrzeni parametrów. Odpowiedzią jest tu raczej pełna gęstość rozkładu prawdopodobieństwa a posteriori $p(\theta|D,M)$, czyli po uzyskaniu danych $D$. \[ p(\theta |D,M) \sim p(D| \theta , M) \cdot p(\theta| M) \]

Zwróćmy uwagę na fakt iż w problemie estymacji parametrów rola wiarygodności globalnej modelu tzn. $ p(D|M) = L(M)$ sprowadza się do czynnika normalizującego i jest zatem nieistotna. Jeśli posiadamy wiedzę o rozkładzie a posteriori, może być ona podsumowana jakąś charakterystyką skupienia rozkładu a posteriori: np. wartością modalną, medianą lub średnią (tj. wartością oczekiwaną)

\[ \langle \theta \rangle = \int \theta \; p( \theta| D,M) \; d\theta \]

Gdy moda i średnia znacząco się różnią, co oznacza, że rozkład jest silnie skośny, stosujemy estymację przedziałową. W statystyce klasycznej estymacja przedziałowa polega na konstrukcji przedziałów ufności CI (ang. confidence intervals). W podejściu Bayesowskim analogicznym tworem jest obszar ufności CR (ang. credible region). Jest to przedział na osi (ogólniej obszar w przestrzeni) parametrów, w którym mieści się określona (zazwyczaj pokaźna, np. $C=0.95$ tj. 95%) część C funkcji rozkładu – ściślej pola pod krzywą funkcji rozkładu – prawdopodobieństwa a posteriori

\[ \int_{CR} p(\theta | D,M) \; d \theta = C \]

Sytuacja jest tu podobna do klasycznego pojęcia przedziału ufności. Jednak filozofia tych dwóch pojęć CI oraz CR jest zasadniczo różna.

Marginalizacja

Często model parametryczny posiada kilka parametrów, a nas interesuje jedynie ich podzbiór. Na przykład mając sygnał okresowy może nas interesować jedynie częstość, niezależnie od wartości amplitudy i fazy sygnału. Analiza Bayesowska daje nam pełną funkcję rozkładu a posteriori dla dowolnej wartości wszystkich parametrów. Informację dotyczącą jedynie interesujących nas parametrów można uzyskać całkowaniem po pozostałych parametrach. Na przykład, jeśli model $ M$ zależy do kilku parametrów (np. $ \theta$ i $\phi$) $M =M(\theta, \phi)$ wtedy

\[ p(\theta | D,M) = \int p(\theta , \phi |D, M) \; d \phi \]

Jest to tzw. prawdopodobieństwo marginalne. Bayesowska marginalizacja rozkładu a posteriori jest uogólnieniem klasycznej teorii propagacji błędów.

Podsumowanie.

Bayesowska procedura estymacji parametrów ${\theta,\phi,\psi} $ modelu $M = M(\theta,\phi,\psi)$, koncepcyjnie sprowadza się do alokacji stopnia racjonalnego zaufania w obrębie ciągłej przestrzeni hipotez i polega na wyznaczeniu łącznego prawdopodobieństwa a posteriori dla parametrów $p(\theta,\phi,\psi|D,M)$ po otrzymaniu danych $D$ z eksperymentu (obserwacji). Wiedza ta jest – w rozważanym kontekście – najpełniejszą z możliwych i dalej może być podsumowana w postaci:

statystyki opisowej – jako np. średnia $\langle \theta \rangle, \langle \phi \rangle,\langle \psi \rangle $ , moda lub mediana (analogicznie do estymacji punktowej w statystyce klasycznej);

jako obszary wiarygodności $CR(\theta), CR(\phi), CR(\psi)$ dla parametrów, czyli najkrótsze przedziały zawierające np. 68% lub 95% rozkładu a posteriori dla $\theta, \phi $ lub $\psi$ (w analogii do estymacji przedziałowej w statystyce klasycznej);

wykresów prawdopodobieństw marginalnych typu $ p(\theta|D,M)$ otrzymanych drogą wycałkowania prawdopodobieństwa łącznego po pozostałych parametrach.

Prawdopodobieństwa aprioryczne

Prawdopodobieństwo aprioryczne modelu czy hipotezy $p(H_i|I)$ jest ważnym elementem wnioskowania Bayesowskiego. Często podejście Bayesowskie bywa krytykowane z powodu stosowania założeń apriorycznych. Podnoszone argumenty są typu: prawdopodobieństwa aprioryczne są subiektywne, teoria – tzn. nauka szczegółowa, w ramach której stosujemy statystykę – nie precyzuje postaci rozkładu apriorycznego, różne założenia a priori prowadzą do różnych wyników a posteriori czyli do różnych wniosków. Krytycyzm ten nie jest tak silny jak mogłoby się wydawać. Prawdopodobieństwa a aprioryczne mają na celu uchwycenie stanu wiedzy o testowanej hipotezie czy modelu. W końcu nie wysuwamy hipotez ani nie budujemy modeli „na oślep”, bezmyślnie czy automatycznie, lecz bazujemy na konkretnej wiedzy i przemyśleniach co składa się na element aprioryczny poznania. Na przykład w problemie estymacji parametrów, wszelkie wcześniejsze próby, badania czy eksperymenty – wszelkie wcześniejsze doświadczenie powinny znaleźć odbicie w postaci rozkładu prawdopodobieństwa a priori. W szczególności fakt, że pewne fizyczne wielkości jak masa, gęstość, objętość, energia itd. są nieujemne musi być należycie uwzględniony. Wnioskowanie Bayesowskie nie jest „czarną skrzynką”. Poza tym potrzeba specyfikacji założeń apriorycznych wymusza ich wysłowienie explicite. Jest to ze wszech miar pożądana cecha. I tak w konkretnych problemach badacz przyjmuje pewne założenia wstępne (czyli de facto aprioryczne) i najczęściej nie są one jawnie wypowiadane – procedura Bayesowska wymusza niejako ich ujawnienie. To, że różny wstępny stan wiedzy (różne prawdopodobieństwa aprioryczne) prowadzi do różnych konkluzji, nie powinno być zaskoczeniem. Natomiast silna zależność wyniku (tzn. prawdopodobieństwa a posteriori) od założeń apriorycznych jest sygnałem słabości procedury badawczej (podsumowanej jako wiarygodność Bayesowska). Innymi słowy w takim przypadku wyniki badań (eksperymentu czy obserwacji) niewiele wnoszą nowego. Przyjęcie określonej postaci rozkładu apriorycznego leży po stronie badacza prowadzącego analizę statystyczną. Jakimi przesłankami powinien się kierować?

Wybór rozkładu prawdopodobieństwa apriorycznego - ogólne wskazówki

W poniższych rozważaniach z uwagi na oszczędność notacji będziemy opuszczać symbol apriorycznej informacji $I$ , która nadaje kontekst badanemu problemowi i od której zawsze uwarunkowane są wszelkie występujące tu prawdopodobieństwa. Należy dążyć do tego aby rozkład prawdopodobieństwa apriorycznego był unormowany tzn. \[ \int p( \theta |M) \; d \theta = 1 \]

Rozkłady takie nazywamy właściwymi (ang. proper prior), czasem jednak konieczne staje się przyjęcie rozkładu nieunormowanego, wówczas mówimy o rozkładzie apriorycznym niewłaściwym (ang. improper prior). Należy następnie pamiętać, że rozkłady aprioryczne nie są niezmiennicze względem nieliniowych transformacji parametrów $ \psi = \psi(\theta)$

\[ p( \psi |M) = p(\theta |M) \left| \frac{d \theta}{d \psi} \right| \] w przypadkach wielowymiarowych moduł pochodnej zastąpić należy jakobianem przekształcenia.

Płaskie rozkłady aprioryczne

W przypadku dyskretnej przestrzeni hipotez, nie posiadając żadnej wstępnej wiedzy preferującej jedne z hipotez nad pozostałe, stosujemy podejście równego rozdziału wstępnego poparcia na wszystkie rozważane możliwości. Na przykład, gdy $H_i \in {H_1, H_2 , \dots , H_m} $ można przyjąć \[p(H_i|M) = 1/m .\] Podobnie w przypadku estymacji parametru, gdy wiemy skądinąd iż jego wartości leżą w pewnym przedziale $ \theta \in [a, b] $ można tą wiedzę aprioryczną opisać rozkładem jednorodnym \[p(\theta|M) = 1 / (b – a).\] Z rozkładami jednorodnymi należy jednak postępować ostrożnie, gdyż jeśli granice przedziału $a$ i $b$ zmierzać będą do nieskończoności rozkład taki stanie się niewłaściwym.

Rozkłady aprioryczne sprzężone

Ważną klasą rozkładów apriorycznych są tzw. rozkłady sprzężone (ang. conjugate priors). Posiadają one tą własność, że rozkłady a posteriori należą do tej samej rodziny funkcji rozkładu co rozkład aprioryczny. Oczywiście konkretne parametry tych rozkładów (apriorycznego i końcowego – a posteriori) będą różne. Rolą Bayesowskiej funkcji wiarygodności jest tu jedynie „uaktualnienie” apriorycznych parametrów modelu, przy zachowaniu jego funkcjonalnej postaci. Z przykładami takich rozkładów sprzężonych spotkamy się poniżej. Rodziny rozkładów sprzężonych były szczególnie ważnymi klasami funkcji w dobie małej mocy obliczeniowej. Obecnie ich znaczenie praktyczne jest mniejsze.

Rozkład Gaussa

Rodzina rozkładów normalnych Gaussa jest przykładem rodziny rozkładów sprzężonych. Załóżmy, że zmienna $X$ ma rozkład normalny $N(\theta,\sigma^2)$ – to jest nasz model statystyczny $M$. Klasycznie powiedzielibyśmy iż tego rozkładu wylosowana została obserwowana wartość zmiennej $X$, tj. dane $D = {x}$ . W języku Bayesowskim powiemy, że \[p(D|M) \equiv p(x|\theta) \sim N(\theta, \sigma^2).\]

Załóżmy też, że rozkład aprioryczny dla parametru $\theta$ jest także rozkładem Gaussa: \[p(\theta) \sim N(\mu, \tau^2).\] Rozkład a posteriori ma zatem postać: \[ p(\theta | D,M) \sim p(D | \theta, M)\cdot p(\theta |M) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{(x- \theta)^2}{2 \sigma^2}}\; \frac{1}{\sqrt{2 \pi \tau^2}} e^{- \frac{(\theta- \mu)^2}{2 \tau^2}} \]

licząc wykładnik eksponenty, po elementarnych, choć nużących przekształceniach dostajemy: \[ - \frac{1}{2 \sigma^2} (x - \theta)^2 - \frac{1}{2 \tau^2} (\theta - \mu)^2 = - \frac{1}{2 \rho^2} [ \theta - \rho^2 (\frac{x}{\sigma^2} + \frac{\mu}{\tau^2}) ]^2 - \frac{1}{2(\sigma^2 + \tau^2)} (x - \mu)^2 \]

gdzie oznaczyliśmy: \[ \rho^2 := \frac{\sigma^2 \tau^2 }{ \sigma^2 + \tau^2}. \] Tak więc rozkład a posteriori jest także rozkładem normalnym \[p(\theta|x,M) \sim N(\nu,\rho^2)\] o średniej \[ \nu = \rho^2 ( \frac{x}{\sigma^2} + \frac{\mu}{\tau^2} ) = \frac{\tau^2}{\sigma^2 + \tau^2} x + \frac{\sigma^2}{\sigma^2 + \tau^2} \mu \] i wariancji \[ \rho^2 = \frac{ \sigma^2 \tau^2}{\sigma^2 + \tau^2}. \]

Z obliczeń powyższych łatwo jest także zobaczyć, że globalna wiarygodność modelu, czyli czynnik normalizacyjny dla rozkładu a posteriori ma także postać rozkładu normalnego: \[ p(D | M) \equiv p(x| M) = \int p(x| \theta, M) p(\theta |M ) \; d \theta \sim N(\mu, \sigma^2 + \tau^2). \]

Jeśli teraz mielibyśmy, zamiast pojedynczej $n$ – obserwacji, tzn. dane $D = {x_1, \dots , x_n}$, czyli w języku klasycznej statystyki mamy $n$ – krotną realizację zmiennej losowej $X$. Rozkład dla statystyki $\bar{x}$ (średniej arytmetycznej) jest rozkładem normalnym \[ p(\bar{x}|\theta) \sim N(\mu,\, \frac{\sigma^2}{n}),\] wobec tego dla rozkładu a posteriori otrzymamy \[ p(\theta \ \bar{x}) \sim N \left( \frac{\tau^2}{\frac{\sigma^2}{n} + \tau^2} \bar{x} + \frac{\frac{\sigma^2}{n}}{\frac{\sigma^2}{n} + \tau^2} \mu \;,\; \frac{\frac{\sigma^2}{n} \tau^2}{\frac{\sigma^2}{n} + \tau^2} \right). \]

Zauważmy, że rozkład a posteriori ma wartość średnią \[ \frac{\tau^2}{\frac{\sigma^2}{n} + \tau^2} \bar{x} + \frac{\frac{\sigma^2}{n}}{\frac{\sigma^2}{n} + \tau^2} \mu, \] która jest średnią ważoną z estymatora $\bar{x} $ średniej, czyli estymatora MLE największej wiarygodności ( w sensie klasycznym) oraz wartości średniej rozkładu apriorycznego $ \mu $. Z postaci wag wynika, że gdy liczebność próby $n$ staje się duża, wpływ założeń apriorycznych staje się mały. I przeciwnie, gdy $n$ jest małe (próba mało liczebna), a nasze aprioryczne przekonania silne (tzn. wariancja $ \tau^2 $ jest mała) wartość średnia a posteriori bliska jest apriorycznej – badania niewiele wniosły do istniejącej już wiedzy.

Rozkład Beta

Przypuśćmy, że zmienna $X$ ma rozkład dwumianowy – w języku Bayesowskim powiemy, że Bayesowska wiarygodność $X$ uwarunkowana parametrem $ \theta$ ma rozkład $B(n,\theta)$: \[ p(x|\theta) = \left( \begin{array}{c} n \\ x \\ \end{array} \right) \theta^x \; (1- \theta)^{n-x} \]

oraz, że rozkład aprioryczny parametru $\theta$ jest rozkładem beta $ Be(\alpha, \beta)$ przy czym hiperparametry $ \alpha $ i $ \beta $ są zadane - znane (hiperparametrami nazywamy parametry, od których zależą rozkłady aprioryczne): \[ p(\theta) = \frac{1}{Be(\alpha, \beta)} \theta^{\alpha - 1} (1 - \theta)^{\beta -1} , \;\; 0 \leq \theta \leq 1, \;\; Be(\alpha, \beta) = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)} \]

Rozkłady: marginalny (tj. globalna wiarygodność) oraz a posteriori mają następującą postać: \[ p(x) = \frac{\left( \begin{array}{c} n \\ x \\ \end{array} \right) Be(x + \alpha, n - x + \beta)}{Be(\alpha + \beta)}, \;\;\; x = 0,1, \dots , n \]

\[ p(\theta |x) = \frac{1}{Be(x + \alpha, n - x + \beta)} \theta^{\alpha + x -1} (1 - \theta)^{n - x + \beta -1} , \;\;\; 0 \leq \theta \leq 1 \]

Sprawdzenie tych faktów pozostawiam jako ćwiczenie. Widać natomiast, że rozkład a posteriori należąc do klasy rozkładów beta jest tego samego typu co rozkład aprioryczny.

Przykład

Rozważmy ponownie następujący problem: rzucamy monetą 4 razy i za każdym razem wypada orzeł. Jakie jest prawdopodobieństwo $ \theta $ otrzymania reszki w pojedynczym rzucie ? Jest to oczywiście klasyczny przykład sprowadzający się do rozkładu dwumianowego z prawdopodobieństwem sukcesu (reszki) w pojedynczej próbie wynoszącym $ \theta $. Zatem Bayesowska wiarygodność jest opisana rozkładem \[ p(x|\theta) \sim B(n,\theta).\] Klasycznie rozumując, estymator największej wiarygodności (MLE) – równy stosunkowi liczby sukcesów do liczby prób, wynosi $\hat{\theta}=0 $. Z Bayesowskiego punktu widzenia, jeśli przyjąć aprioryczny rozkład dla parametru $\theta $ jako jednorodny $ p(\theta) \sim U(0,1)$ to rozkład a posteriori będzie dany w postaci

\[ p(\theta|x) \sim \theta^0 (1 – \theta)^4. \]

Zatem będzie on typu rozkładu beta $Be(1,5)$

Z własności rozkładu beta wynika \[\langle \theta \rangle_{Beta} = \frac{1}{1+5} = \frac{1}{6}\] co jest znacznie rozsądniejszym wynikiem niż klasyczny estymator MLE.

Przykład

Badanie kliniczne pokazało, że w grupie 980 porodów powikłanych, 437 urodzeń przypadało na dziewczynki. W populacji generalnej częstość urodzeń się dziewczynek wynosi 485 na 1000 porodów. Wyniki badania sugerują, że odsetek dziewczynek w grupie porodów powikłanych 0.446 jest niższy niż w populacji generalnej. Pytanie: Jak znamienna jest ta różnica?

Rozwiązanie

Problem częstości urodzeń jest klasycznym przykładem zastosowania rozkładu Bernoulliego $ B(x;n,\theta) $, gdzie $ \theta $ oznacza prawdopodobieństwo "sukcesu" (w naszym przypadku - urodzenia się dziewczynki) w pojedynczej "próbie'".

Bayesowska wiarygodność modelowana jest zatem rozkładem Bernoulliego \[ p(D|\theta) \sim \theta^{437} (1-\theta)^{543} \sim Be(\theta|437,543). \] Jako rozkład a priori parametru $ \theta $ przyjmijmy rozkład jednorodny w przedziale $[0,1]$ czyniąc zadość zasadzie ekwipartycji racjonalnego zaufania. Szczęśliwą okolicznością jest iż rozkład jednorodny należy do rodziny rozkładów Beta \[ U[0,1] = Be(\theta|1,1). \]

Ponieważ rozkład Beta jest sprzężony z rozkładem Bernoulliego, rozkład a posteriori będzie także rozkładem Bernoulliego \[ p(\theta|D) \sim \theta^{438} (1-\theta)^{544} \sim Be(\theta|438,544) \] Pisząc powyższe wykorzystaliśmy wiedzę z rozdziału 6.5.

Znamienność różnicy pomiędzy odsetkiem dziewczynek rodzących się w powikłanych porodach w porównaniu z populacją generalną najlepiej zobrazuje rysunek przedstawiający gęstość rozkładu a posteriori. Ciągła linia pionowa oznacza wartość częstości urodzeń dziewczynek w populacji generalnej, natomiast linie przerywane pokazują granice 95 % obszaru ufności dla parametru $ \theta $.

Poniżej podany jest program w języku środowiska R, który posłużył do wykonania rysunku.

>theta=seq(0.3,0.6,0.001)
 
>plot(theta,dbeta(theta,437,543),type="l",xlab="theta",ylab="",yaxt="n",bty="n",cex=2) 
 
>lines(theta,dbeta(theta,1,1),lty=2)  
 
>lines(theta,dbeta(theta,438,544), lty=2)  #posterior Be(438,544)
 
>abline(v=0.485)
 
>CI <- quantile(rbeta(1000,438,544),probs=c(0.025,0.975))
 
>abline(v=CI[1],lty=2)
 
>abline(v=CI[2],lty=2)

Widzimy, że częstość urodzeń dziewczynek w populacji generalnej (pionowa linia ciągła) leży na zewnątrz 95% obszaru ufności rozkładu a posteriori (pionowe linie przerywane).

Ćwiczenie

W poprzednim przykładzie przyjęliśmy aprioryczny rozkład parametru $ \theta $ jako płaski argumentując to zasadą ekwipartycji racjonalnego poparcia. Zaproponuj inne postacie rozkładu a priori (w ramach rodziny rozkładów Beta). Używając odpowiednio zmodyfikowanego programu przetestuj wpływ na wynik założeń donośnie rozkładu a priori.

Twierdzenie

Niech zmienna $X $ posiada rozkład należący do rodziny rozkładów wykładniczych, tzn. przy \[ X = (x_1, \dots ,x_n) \]

\[ p(x_i|\theta) \sim A(\theta)\; e^{T^{\star}(x_i) B(\theta)} \; \Psi(x_i) \]

Łączna funkcja wiarygodności jest postaci \[p(D|M) = p(x|\theta) = \prod_{i=1}^{n} p(x_i|\theta) =[A(\theta)]^n\; e^{T\;B(\theta)} \; H(x) \] gdzie \[ T = \sum_i T^{\star}(x_i) \] oraz \[H(x) = \prod_{i=1}^{n} \Psi(x_i).\] Naturalny sprzężony rozkład aprioryczny ma postać \[p(\theta) \sim [A(\theta)]^p \; e^{q\; B(\theta)}\]

i wówczas rozkład a posteriori jest typu \[p(\theta|x) \sim [A(\theta)]^{n+p}\; e^{(T + q)\; B(\theta)}.\]

Rodzina wykładnicza postulowana w powyższym twierdzeniu jest dość ogólna – zawiera w sobie na przykład omawiany wcześniej przypadek rozkładu dwumianowego. Mamy wtedy mianowicie: wiarygodność \[ p(x|\theta) \sim B(n, \theta)\] oraz rozkład aprioryczny \[ p(\theta) \sim Be(\alpha, \beta).\]

Przedstawmy wiarygodność w postaci

\[ p(x | \theta) = \left( \begin{array}{c} n \\ x \\ \end{array} \right) \theta^x\; (1 - \theta)^{n-x} = \left( \begin{array}{c} n \\ x \\ \end{array} \right) \left( \frac{\theta}{1 - \theta} \right)^x\; (1 - \theta)^n = (1 - \theta)^n \; e^{x \ln{\frac{\theta}{1 - \theta}}} \; \left( \begin{array}{c} n \\ x \\ \end{array} \right) \] Spełnione są założenia Twierdzenia, przy oznaczeniu

\[ A(\theta) = (1- \theta),\] \[ T = x, \] \[ B(\theta) = \ln{(\theta / (1-\theta))} \] oraz \[ H(x) = \left( \begin{array}{c} n \\ x \\ \end{array} \right) \]

Rozkład aprioryczny jest typu

\[ p(\theta) \sim \theta^{\alpha - 1} (1 - \theta)^{\beta - 1} = (1 - \theta)^{\alpha + \beta -2} \left( \frac{\theta}{1 - \theta} \right)^{\alpha -1} = (1 - \theta)^{\alpha + \beta -2} \; e^{(\alpha -1) \ln{\frac{\theta}{1 - \theta}}} \]

wobec tego jest on naturalnym rozkładem sprzężonym przy oznaczeniach

\[A(\theta)^p = (1 – \theta)^{\alpha + \beta -2}\] oraz

\[ q = \alpha -1 \]

i na mocy ogólnego twierdzenia rozkład a posteriori ma postać

\[ p(\theta | x) \sim (1- \theta)^{n+\alpha+\beta-2}\;e^{(x + \alpha -1) \ln{\frac{\theta}{1 - \theta}}} = \theta^{x + \alpha - 1} (1 - \theta)^{n + \alpha + \beta - 2 - x - \alpha + 1} = \theta^{x + \alpha -1 } (1 - \theta)^{n + \beta - x - 1} \]

Czyli jest rozkładem beta \[Be(\alpha+x, n-x+\beta).\]

Rozkłady aprioryczne dla parametrów położenia i skalowania

Parametr $\theta$ nazywamy parametrem położenia (ang. location parameter) jeżeli funkcja wiarygodności jest niezmiennicza względem translacji \[ p(x| \theta) \sim f(x -\theta) \]

Od rozkładu apriorycznego oczekujemy także niezmienniczości względem translacji: \[ \forall \theta_0 \;\;\; p(\theta) = p(\theta - \theta_0) \]

Jedyną możliwą tu postacią funkcyjną jest funkcja stała: \[ p(\theta)=c \] czyli naturalnym rozkładem apriorycznym dla parametrów położenia jest rozkład jednorodny. Jest on rozkładem apriorycznym niewłaściwym, gdyż jeśli dziedziną określoności parametru $ \theta $ jest $R$, wówczas takiego rozkładu apriorycznego nie można unormować. Rozkłady aprioryczne niewłaściwe bywają stosowane, gdy wynikają z symetrii problemu (np. jak w naszym przypadku parametru położenia). Od rozkładu niewłaściwego żąda się aby rozkład marginalny (globalna wiarygodność) był skończony: . \[ p(x) = \int p(x| \theta) p(\theta) \; d \theta < \infty \]

Parametr $ \theta $ nazywamy parametrem skalowania (ang. scale parameter) jeżeli funkcja wiarygodności jest typu: \[ p(x | \theta) \sim \frac{1}{\theta} f \left( \frac{x}{\theta}\right) \]

Rozkład aprioryczny powinien więc być niezmienniczy względem przeskalowań: \[ \forall c >0 , \;\;\; p(\theta) = \frac{1}{c} p\left( \frac{\theta}{c} \right) \]

Rozwiązaniem tego równania funkcyjnego jest \[ p(\theta) = \frac{1}{\theta} \] i jest to oczywiście rozkład niewłaściwy.

Rozkład aprioryczny Jeffreys’a

Sir Harold Jeffreys wychodząc z klasycznej teorii funkcji wiarygodności podał pewną propozycję konstruowania rozkładów apriorycznych nazywaną obecnie jego nazwiskiem. Rozważmy, mianowicie Bayesowską funkcję wiarygodności $ p(x|\theta)$ i traktując ją jak klasyczną funkcję wiarygodności policzmy informację Fishera \[ I(\theta) = - E \left( \frac{\partial^2 log(p(x|\theta))}{\partial \theta^2} \right) \]

gdzie wartość oczekiwana $ E[ \cdot ] $ wzięta została w sensie rozkładu $ p(x|\theta)$. Klasycznie informacja Fishera jest miarą krzywizny funkcji wiarygodności w otoczeniu ekstremum tzn. estymatora $ \hat{\theta}$ największej wiarygodności (MLE), czyli innymi słowy mierzy czułość estymatora w otoczeniu $ \hat{\theta}$ . Jeffreys zaproponował, aby rozkład aprioryczny właściwy dla wiarygodności $ p(x|\theta)$ przyjąć jako:

\[ p(\theta) \sim \det{(I(\theta))}^{1/2} \]

oczywiście wyznacznik ma zastosowanie w problemach wielowymiarowych, dla problemu jednowymiarowego informacja Fishera nie ma reprezentacji macierzowej, lecz jest po prostu funkcją parametru. Wybór pierwiastka kwadratowego wynika z zachowania się względem reparametryzacji, mianowicie jeśli $ \psi = h(\theta) $, przy czym funkcja $ h $ jest odwracalna o znaczy istnieje $ \theta = h^{-1}(\psi) $, wówczas \[ p(\psi) = p(\theta) \left| \frac{d \theta}{d \psi} \right| \]

Ale \[ I(\psi) = - E \left[\frac{\partial^2 log(p(x|\psi))}{\partial \psi^2} \right] = - E \left[ \frac{\partial^2 log(p(x|\theta))}{\partial \theta^2} \left| \frac{d \theta}{d \psi} \right|^2 \right] = I(\theta)\; \left| \frac{d \theta}{d \psi} \right|^2 \]

Stąd \[ (I(\theta))^{1/2} = (I(\psi))^{1/2} \left| \frac{d \theta}{d \psi} \right| \]

Przykłady apriorycznych rozkładów Jeffreys’a

Niech $ X_1, \dots ,X_n $ mają Gaussowski rozkład, więc ich łączna wiarygodność jest proporcjonalną do \[ p(x| \theta) \sim e^{ - n( \bar{x} - \mu)^2 / 2 \sigma^2}. \] Jeśli wariancja $ \sigma^2 $ jest znana i ustalona, a parametrem estymowanym jest $ \theta = \mu $ , wówczas \[ I(\theta) = n / \sigma^2 \] czyli rozkład aprioryczny Jeffreys’a jest płaski: \[ p(\theta) \sim (n/ \sigma^2)^{1/2} \sim 1 \]

Wynik jest zgodny z dyskusją poprzedniego paragrafu – parametr $ \theta $ jest parametrem położenia.

Gdyby z kolei znana była wartość $ \mu $ natomiast estymowanym nieznanym parametrem była wariancja $ \theta = \sigma^2 $ wówczas łączna funkcja wiarygodności: \[ p(x| \theta) \sim \theta^{n/2} \; e^{- s /(2 \theta)} \]

gdzie \[ s := \sum_i (x_i - \mu)^2 \] oraz \[ I(\theta) = - E [ n/ (2 \theta^2) - s / \theta^3 ]. \] Ale $ E[s] = n \theta $ więc ostatecznie rozkład a aprioryczny Jeffreys’a ma postać \[ p(\theta) \sim 1 / \theta \] - znowu w zgodzie z faktem, że wariancja jest parametrem skalowania.

Zadanie

Znaleźć aprioryczny rozkład Jeffreys’a dla wiarygodności Gaussowskiej, gdy oba parametry: średnia $ \mu $ i wariancja $ \sigma^2 $ są nieznane (tzn. mają być estymowane).

Odpowiedź: $ p(\theta ,\sigma^2)=1/ \sigma^2 $

Zadanie

Znaleźć aprioryczny rozkład Jeffreys’a dla wiarygodności opisanej rozkładem Poissona $ p(x | \theta) = P_x(\theta) = \frac{\theta^x}{x !} e^{-\theta} $.

Odpowiedź

$ p(\theta) = \frac{1}{\sqrt{\theta}} $

Wiarygodność Bayesowska

Oprócz prawdopodobieństw apriorycznych $ p(D|I)$ kolejnym ważnym elementem jest Bayesowska wiarygodność $ p(D|M,I)$. Podobnie jak w poprzednim rozdziale element uwarunkowania prawdopodobieństw aprioryczną informacją $I$ będziemy opuszczać w notacjach, stale mając jednak ten fakt w pamięci. W przypadku gdy dane $D$ stanowią obserwowane (mierzone) wartości zmiennej $ X$, a model zjawiska $M$ polega na przypuszczeniu iż zmienna $ X$ posiada (konkretny, zadany) rozkład o parametrach $\theta $ będziemy pisać $ p(x|\theta)$.

Tak więc sens wiarygodności Bayesowskiej: prawdopodobieństwo (miara racjonalnego oczekiwania) wygenerowania danych $ x$ z modelowego rozkładu parametryzowanego przez $ \theta$ jest bardzo zbliżony do klasycznej koncepcji rozkładu, z którego wylosowana została próba (ang. sampling distribution). Dla potrzeb obliczeń to jest ten sam rozkład! Zasadnicza różnica jest taka iż w przypadku Bayesowskim, wartości $ x $ (danych $ D$) są konkretne i ustalone – te mianowicie, które zostały faktycznie zaobserwowane. Z tego też powodu (zależności od konkretnych, a nie hipotetycznych – możliwych do zaobserwowania – danych) wiarygodności Bayesowskie nie są rozkładami prawdopodobieństwa sensu stricte. Wiarygodność Bayesowska jest ściśle związana (de facto tożsama) z klasyczną funkcją wiarygodności. Mianowicie w statystyce klasycznej funkcją wiarygodności nazywamy funkcję gęstości rozkładu prawdopodobieństwa zmiennej losowej rozpatrywaną jako funkcję parametrów od której ona zależy. Z rozkładami modelującymi wiarygodność Bayesowską spotkaliśmy się już w rozdziale poświęconym rozkładom apriorycznym. Właściwy wybór rozkładu apriorycznego ściśle związany jest bowiem z postacią wiarygodności. Zakładając znajomość statystyki klasycznej przez czytelnika – przedstawimy poniżej listę najczęściej spotykanych rozkładów mogących służyć do modelowania wiarygodności Bayesowskiej.

Rozkład normalny Gaussa

Stanowi dobre przybliżenie gdy „szum” odpowiedzialny za niepewność, czyli brak idealnej zgodności między wartościami mierzonymi i przewidywaniami teorii, pochodzi od superpozycji wielu niezależnych zmiennych losowych. Centralne Twierdzenie graniczne uprawomocnia to założenie.

Rozkład logarytmiczno-normalny

Stosowany w modelowaniu rozkładów skośnych, gdy dobrą „receptą” na odzyskanie symetrii jest transformacja logarytmiczna oryginalnej zmiennej. Wiele interesujących fizycznie wielkości może przyjmować jedynie wartości dodatnie i z tego powodu (zawężenia dziedziny do liczb rzeczywistych nieujemnych) rozkłady takich zmiennych są skośne.

Rozkład Poissona

Opisuje, w danym przedziale czasu, liczbę dyskretnych zdarzeń występujących z pewną średnią szybkością. Np. liczba jąder ulegających rozpadowi promieniotwórczemu, zliczenia fotonów w detektorze, liczba klientów odwiedzających sklep w danym przedziale czasu, liczba rozmów telefonicznych w danym przedziale czasu itp. Rozkład Poissona jest też granicznym przypadkiem rozkładu Bernouliego przy małych wartościach prawdopodobieństwa sukcesu w pojedynczej próbie i dużych wartościach liczby prób. Z tego powodu rozkładem Poissona modelujemy np. dobową umieralność w populacji ludzkiej.

Rozkład wykładniczy

W technice opisuje niezawodność urządzenia przy zadanej intensywności uszkodzeń (awarii) Inne zastosowania to modelowanie czasu oczekiwania na zgłoszenie się operatora w centrali telefonicznej, odstęp czasowy między wypadkami na kopalni lub trzęsieniami ziemi,

Rozkład dwumianowy (Bernoulliego)

Jest to jeden z najbardziej użytecznych rozkładów dyskretnych. Ma on zastosowanie do wszelkiego rodzaju powtarzalnych zdarzeń o naturze dychotomicznej, tzn. typu: wystąpienie lub niewystąpienie danego zdarzenia. Jako (ciągłe) graniczne przypadki rozkładu dwumianowego wyprowadzić można rozkłady Gaussa i Poissona. Rozkładem dwumianowym modeluje się np. ekspresję genów, rozkład liczby dzieci w rodzinie, prognostyczne znaczenie czynników ryzyka.

Rozkład jednorodny

Modeluje on rozkład błędów zaokrąglania cyfr znaczących w rozwinięciu dziesiętnym. Jest często stosowanym rozkładem apriorycznym modelującym naszą niewiedzę.

Rozkład Beta

Jest to w zasadzie dwu-parametrowa rodzina rozkładów $ Be(x |\alpha, \beta) $, gdzie $ \alpha $ i $ \beta $ są tzw. parametrami kształtu (ang. shape parameters). Rodzina rozkładów Beta jest bardzo użyteczna do modelowania rozkładów apriorycznych sprzężonych z rozkładem Bernoulliego. Zawiera w sobie rozkład jednorodny na przedziale $ [0,1] $ realizowany jako $ Be(x|1,1) $ i jak pokazuje rysunek odpowiednio dobierając parametry kształtu, rozkładem Beta można modelować zarówno symetryczne jak i skośne rozkłady.

Rozkład Cauchy’ego

Rozkład ten jest przykładem rozkładu nie posiadającego wartości oczekiwanej – stosowna całka jest nieokreślona. Nieparametryczne miary tendencji centralnej tj. moda i mediana są tu oczywiście dobrze określone obrazując przewagę podejść nieparametrycznych. Rozkład ten nie jest jedynie czysto akademicki – w fizyce atomowej modeluje on szerokość własną linii widmowych a w fizyce cząstek elementarnych – profil jądrowych poziomów energetycznych wzbudzonych oraz rezonansów cząstkowych. Zwany jest wówczas profilem Lorentza lub rozkładem Breita-Wignera.

Rozkład Weibulla

Stosowany jest w analizie przeżywalności do modelowania czasu dalszego życia po przeżyciu określonego czasu. Analiza przeżywalności (survival analysis ) jest kluczowym zagadnieniem dla firm ubezpieczeniowych. Skuteczność terapeutyczna i prognostyka pewnych zabiegów chirurgicznych, zabiegów typu wszczepienia stymulatora serca, plastyki tętnic wieńcowych jest często obrazowana w terminach analizy przeżywalności.

Metodyka uzyskiwania rozkładów a posteriori

Centralnym wzorem statystyki Bayesowskiej jest

\[ p(\theta | D) \sim p(D | \theta) p(\theta) \]

W jaki sposób przedstawiać prawą stronę tej zależności, omawialiśmy w poprzednich rozdziałach. Jak zatem uzyskiwać rozkład a posteriori? Pewne analityczne metody – oparte na koncepcji sprzężonych rozkładów apriorycznych – zostały już omówione. Miały one bardzo duże znaczenie w statystyce Bayesowskiej przed upowszechnieniem się komputerów. Obecnie w przypadkach ogólnych stosujemy wnioskowanie oparte na metodzie symulacji Monte Carlo – tzw. MCMC (ang. Markov Chain Monte Carlo – metoda Monte Carlo przy użyciu łańcuchów Markowa). Metoda MCMC pozwala wygenerować listę wartości $ {\theta_i}$ próbkujących rozkład a posteriori.

Dalsza idea jest prosta: wartości oczekiwane parametrów, lub funkcji zależnych od parametrów uzyskujemy uśredniając po próbkach wylosowanych w rozkładu a posteriori:

\[ \langle \theta \rangle = \int \theta p(\theta | D, M)\; d \theta \approx \frac{1}{N} \sum_i \theta_i \]

\[ \langle f(\theta) \rangle = \int f(\theta) p(\theta | D, M)\; d \theta \approx \frac{1}{N} \sum_i f(\theta_i) \]

Metody próbkowania

W najprostszych przypadkach, gdy przestrzeń parametrów ma niski wymiar skuteczne okazują się proste metody próbkowania: metoda transformacji odwrotnej oraz metoda odrzucania.

Metoda transformacji odwrotnej – efektywna dla rozkładów jednowymiarowych – polega na spostrzeżeniu, że dystrybuanta rozkładu jest funkcją przyjmującą wartości od 0 do 1. Można zatem utworzyć, przy pomocy generatora liczb losowych, zbiór $ {x_i}$ jednorodnie próbkujący przedział [0,1] a następnie traktując te liczby jako wartości dystrybuanty $ F(\theta)$ próbkowanego rozkładu przekształcić je (odwrócić) w próbkę $ {\theta_i}$ zgodnie ze wzorem: \[ F(\theta_i)=x_i\].

Przykład

Rozkład wykładniczy \[ P(\theta) = 1/2 exp( - \theta/2)\] dla \[ 0 \leq \theta < 1 \] ma dystrybuantę $ F(\theta) = 1 - exp( - \theta/2)$.

Eksplorujemy go tworząc jednorodną próbkę \[ {x_i}\;\; 0 \leq x_i < 1 \] a następnie licząc \[ \theta_i = - 2 ln(1-x_i) \]. Tak uzyskane liczby tworzą próbkę wylosowaną z wyjściowego rozkładu.

Metodę odrzucania dobrze ilustruje problem losowego próbkowania koła jednostkowego. Dysponując standardowym generatorem liczb losowych z przedziału $ [-1,+1] $ łatwo tworzymy punkt $(x, y)$ leżący we wnętrzu kwadratu $ [-1,1] \times [-1,1] $. Sprawdzamy następnie, czy $ x^2 + y^2 < 1 $, jeśli tak zatrzymujemy go na liście, gdy nie odrzucamy i generujemy kolejny punkt, itd.

Porównując powierzchnie koła jednostkowego i opisanego na nim kwadratu widzimy że procedura ta generuje „dobry” punkt z częstością $ \pi/4 $ (wydajność ok. 75%). Problemy pojawiają się w wyżej wymiarowych przestrzeniach parametrów. Gdyby ten sam problem – próbkowania kuli jednostkowej – przenieść do $ n $ wymiarów, łatwo byłoby wygenerować rozkład punktów próbkujących losowo kostkę $n$-wymiarową: $ [-1,1] \times \dots \times [-1,1]$. Objętość $n$ - wymiarowej kostki to $ 2^n $, podczas gdy objętość wpisanej weń $n$ - wymiarowej kuli to $ V_n = 2 \pi^{n/2} / n \Gamma(n/2) $ . Dzieląc te objętości przez siebie, można się przekonać iż stosunek objętości kuli do kostki szybko maleje ze wzrostem wymiaru przestrzeni. I tak gdy wydajność próbkowania sfery metodą odrzucania wynosi $ \pi/4$ dla $ n=2 $, to już dla n=10 jest to ok. 0.025, a dla $ n=100$ jest rzędu $ 10^{-70}$.

Metoda Łańcuchów Markowa.

Metoda ta (MCMC) jest najpowszechniej stosowaną obecnie techniką próbkowania wyżej wymiarowych przestrzeni parametrów. W skrócie: jej idea polega na znalezieniu takiego sposobu eksploracji przestrzeni parametrów (operacyjnie – ciągu punktów) że po dostatecznie długim czasie podróży (operacyjnie dla dostatecznie dużej liczby iteracji $N $ ) prawdopodobieństwo iż znajdujemy się w punkcie $ \theta_i $ jest proporcjonalne do $ p(\theta_i)$ czyli do wartości próbkowanej gęstości prawdopodobieństwa w tym punkcie. Jeżeli reguła przejścia od $ \theta_i$ do $ \theta_{i+1} $ – oznaczmy ją $ T(\theta_i, \theta_{i+1})$ – zależy jedynie od $ \theta_i$ a nie od przeszłej historii naszej podróży, wówczas ciąg generowany taką regułą iteracyjną nazywamy łańcuchem Markowa (ang. Markov Chain). Aby taki ciąg – łańcuch Markowa, miał pożądane własności graniczne prawdopodobieństwo dojścia do punktu $ \theta_i $ musi być równe gęstości próbkowanego prawdopodobieństwa w tym punkcie, czyli innymi słowy: \[ p(\theta_k) = \int p(\theta_i) T(\theta_i , \theta_k) d \theta_i \]

Jednym ze sposobów zapewnienia tej własności jest wybór reguł przejścia spełniających warunek: \[ p(\theta_{i+1}) T(\theta_{i+1}, \theta_i) = p(\theta_i) T(\theta_i, \theta_{i+1} ) \]

Jednym z najpopularniejszych algorytmów konstrukcji łańcuchów Markowa o powyższej własności jest algorytm Metropolisa – Hastingsa. Bardzo dobre pozycje na temat metody MCMC (w języku angielskim) są dostępne „online” – patrz Literatura.

Uwagi

1. Należy pamiętać, że w łańcuchu Markowa, chociaż $ p(\theta_i)$ asymptotycznie dobrze reprezentuje próbkowany rozkład, to kolejne wartości $p(\theta_i)$ oraz $ p(\theta_{i-1}) $ nie są niezależne. Uruchomienie łańcucha i jego zatrzymanie po osiągnięciu równowagi generuje tylko jedną niezależną próbkę.

2. Dla osiągnięcia celu – czyli zbudowania dużej (często o liczebności kilku milionów punktów) próbki musimy za każdym razem uruchamiać kolejny łańcuch.

3. Przy kolejnych uruchomieniach łańcucha należy odczekać aby utracona została informacja o poprzedniej próbce

4. Zawsze trzeba śledzić autokorelacje w próbkach.

5. Nietrywialnym elementem metody MCMC jest reprezentatywność próbki tj. czy łańcuchy Markowa wyczerpująco eksplorują przestrzeń parametrów.

Selekcja modeli

W dotychczasowych rozważaniach koncentrowaliśmy uwagę na problemie estymacji parametrów postrzeganym jako Bayesowskie wnioskowanie w warunkach ciągłej przestrzeni hipotez. Zakładaliśmy słuszność modelu $ M$ parametryzowanego przez $ \theta_i $ . Typowo kilka parametrycznych modeli $ M_i $ (zależnych od różnych zestawów parametrów i różniących się stopniem złożoności) może opisywać badane zjawisko (lepiej rzec kandydować do opisu badanego zjawiska). Podejście Bayesowskie umożliwia hierarchizację poparcia dawanego poszczególnym modelom przez dane eksperymentalne. Jest to przedmiotem tzw. teorii selekcji modeli. W klasycznej statystyce konfrontacja kilku alternatywnych modeli z danymi eksperymentalnymi prowadzi do estymacji ich parametrów (np. drogą minimalizacji funkcji chi-kwadrat lub maksymalizacji funkcji wiarygodności). W ten sposób dostajemy zestaw modeli, z których każdy jest „najlepiej dopasowany” do danych. W ramach klasycznej statystyki nie ma możliwości (tj. rzetelnych podstaw teoretycznych) aby stwierdzić który z owych najlepiej dopasowanych modeli uzyskuje najlepsze poparcie przez dane. Wyjątek stanowi teorio - informacyjne podejście tzw. kryterium Akaike łączące w sobie koncepcje funkcji wiarygodności i tzw. entropii Kullbacka-Leiblera. Pole stosowalności kryterium Akaike jest jednak ograniczone, w porównaniu do metod Bayesowskich.

Podejście Bayesowskie, oparte na rozumieniu prawdopodobieństwa a posteriori modelu jako miary racjonalnego zaufania w słuszność tego modelu w świetle otrzymanych danych, dostarcza spójnej koncepcji dla teorii selekcji modeli. Bayesowska selekcja modeli wymaga specyfikacji modeli $ M_i$ , wraz z ich apriorycznymi prawdopodobieństwami $ p(M_i |I)$ (wracamy w notacji do symbolu $ I $ obrazującego całokształt dotychczasowej wiedzy i kontekstu problemu). W ramach apriorycznej informacji I zakładamy, że rozważane modele wyczerpują one przestrzeń możliwych modeli.

W tym miejscu należy się bardzo istotna uwaga: wynikająca w procesie selekcji modeli hierarchizacja będzie miała zawsze charakter względny. Wskaże ona najlepiej popierany przez dane, model z rozważanego zestawu. Nie oznacza to, że musi to być model „prawdziwy” w sensie absolutnym. Rzeczywistość może być taka, że prawdziwego modelu badanego zjawiska po prostu nie znamy (najczęściej tak właśnie jest) i rozważamy jedynie pewne przybliżenia.

Stosujemy wzór Bayesa \[p(M_i|D, I) = \frac{p(D|M_i I)}{p(D|I)} p(M_i|I) \]

Prawdopodobieństwa p(D|Mi,I) są tu wiarygodnościami modeli, tzn. \[ p(D|M_i) = \int p(D|\theta, M_i, I) p(\theta |M_i, I)\; d \theta \]

Do porównania modeli stosujemy iloraz szans: \[ O_{ij} = \frac{p(M_i |D,I)}{p(M_j |D,I)} = \frac{p(M_i|I)}{p(M_j|I)} \frac{p(D|M_i,I)}{p(D|M_j,I)} = \frac{p(M_i|I)}{p(M_j|I)}\; B_{ij} \] Czynnik normalizacyjny się upraszcza i iloraz szans staje się iloczynem ilorazu prawdopodobieństw a priori oraz tzw. czynnika Bayesa $ B_{ij}$, który jest stosunkiem wiarygodności modeli. Jeżeli prawdopodobieństwa aprioryczne modeli są równe – a tak często się zakłada dokonując ekwipartycji racjonalnego zaufania pomiędzy wszystkie konkurencyjne modele – wówczas iloraz szans jest liczbowo równy czynnikowi Bayesa.

Gdy znamy iloraz szans $ O_{i1} $ modeli $ M_i$ względem pewnego modelu M1, możemy obliczyć prawdopodobieństwa każdego z modeli – czyli innymi słowy stopień racjonalnego poparcia dawanego każdemu z modeli przez dane $ D$.

Mamy mianowicie: \[ \sum_{i=1}^{N_{mod}} p(M_i |D,I) = 1 \]

Dzieląc to wyrażenie przez p(M1|D,I) dostajemy: \[ \sum_{i=1}^{N_{mod}} O_{i1} = \frac{1}{p(M_1 |D,I)} \]

i dalej : \[ \frac{p(M_i|D,I)}{p(M_1 | D,I)} =: O_{i1} = p(M_1 |D,I) \sum_{j=1}^{N_{mod}} O_{j1} \]

Stąd ostatecznie: \[ p(M_i|D,I) = \frac{O_{i1}}{\sum_{j=1}^{N_{mod}} O_{j1}} \]

Jest oczywistym, że $ O_{11} = 1 $. Gdy mamy do porównania dwa modele, to:

\[ p(M_2|D,I) = \frac{O_{21}}{1+O_{21}} \]

W charakterze przykładu rozważmy następującą sytuację: mamy dwa konkurujące modele $M_1$ i $M_2$. Model $M_1$ przewiduje wartość pewnej wielkości jako $d_1 = 100$, natomiast model $M_2$ przewiduje $d_2 = 200$. Załóżmy następnie, że niepewność $u$ pomiaru tej wielkości dana jest rozkładem Gaussa $ p(u) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{u^2}{2 \sigma^2}} $, przy czym odchylenie standardowe wynosi $ \sigma = 40 $.

Załóżmy teraz, że pomiar owej wielkości dał wynik $d = 120 $. Porównajmy teraz oba modele sposobem Bayesowskim, licząc iloraz szans $O_{12} $. Nie mając apriorycznych przesłanek faworyzujących któryś z modeli, przypisujemy im równe prawdopodobieństwa a priori $p(M_1) =p(M_2) = 0.5 $.

Oznacza to, że licząc iloraz szans, prawdopodobieństwa aprioryczne i czynnik normalizacyjny $ p(D)$ ulegną redukcji

\[O_{12} = \frac{p(M_1|D)}{p(M_2|D))} =\frac{p(D|M_1) p(M_1)}{p(D|M_2) p(M_2)} = \frac{p(D|M_1)}{p(D|M_2))} \]

W tym przypadku iloraz szans modeli równa się ilorazowi ich wiarygodności.

Liczymy wiarygodności - dla pierwszego modelu

\[ p(D|M_1) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{(d-d_1)^2}{2 \sigma^2}} = \frac{1}{\sqrt{2 \pi}\, 40} e^{- \frac{(120-100)^2}{2 (40)^2}} = 0.00880 \]

oraz dla drugiego

\[ p(D|M_1) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{(d-d_2)^2}{2 \sigma^2}} = \frac{1}{\sqrt{2 \pi}\, 40} e^{- \frac{(120-200)^2}{2 (40)^2}} = 0.00135 \].

Stąd iloraz szans $O_{12} = 6.52$, co oznacza wyraźne poparcie dla modelu $ M_1$. Opisaną sytuację ilustruje poniższy rysunek.

Literatura

1. W polskojęzycznej literaturze statystykę Bayesowską omawia rozdział 9 książki Romana Nowaka „Statystyka dla fizyków”(PWN, Warszawa 2002). Nie porusza on wielu omawianych w niniejszym skrypcie problemów, jednak zawiera cenną, a tutaj pominiętą, dyskusję wyprowadzenia rachunku prawdopodobieństwa jako sytemu logicznego wnioskowania w warunkach niepewności. Godne uwagi są też propozycje R. Nowaka aby "prawdopodobieństwo aprioryczne" nazywać "zaczątkowym", a "prawdopodobieństwo a posteriori" – "wynikowym".

2. Bogate linki do literatury o tematyce Bayesowskiej zawierają strony:

[1] International Society for Bayesian Analysis

[2] BIPS - Bayesian Inference for Physical Sciences - strona prowadzona przez Toma Loredo - zawiera bogate odnośniki do wykładów i prac oryginalnych stosujących statystykę Bayesowską w astrofizyce i fizyce. Szczególnie godne polecenia są wykłady CERNowskie d'Agostiniego ze statystyki Bayesowskiej.

3. Algorytmy MCMC – literatura “online”:

[3] MacKay, D.J.C. “Information Theory, Inference and Learning Algorithms” Cambridge University Press 2003

[4] Neal R.M. 1993 Probabilistic inference using Markov Chain Monte Carlo methods,

4. Awangardowe zastosowania metod Bayesowskich (w szczególności eksploracji przestrzeni parametrów, marginalizacji i selekcji modeli ) na polu współczesnej kosmologii oraz fizyki cząstek elementarnych ( ciemna materia w supersymetrycznym rozszerzeniu Modelu Standardowego), znaleźć można na stronach: [5] SuperBayeS - Supersymmetry Parameters Extraction Routines for Bayesian Statistics, oraz [6] - Cosmological Monte Carlo

5. W celu ułatwienia uprawiania statystyki Bayesowskiej w oparciu o MCMC, opracowany został publicznie dostępny program WinBUGS. Jest on dostępny darmowo na stronie: [7]

6. Rosnącą popularnością wśród statystyków cieszy się środowisko R. Program instalacjyny, wszelkie informacje związane z R w tym linki do obszernych bibliotek programów zawiera strona [8].