Markov Chain Monte Carlo (MCMC) je sve popularnija metoda za dobijanje informacija o raspodelama, posebno za procenu aposteriornih raspodela u Bajesovskom pristupu. Sada ćemo opisati šta je zapravo MCMC i za šta se može koristiti, sa jednostavnim ilustrativnim primerima. Istaći ćemo neke od prednosti i ograničenja za MCMC metodu uzorkovanja, kao i različite pristupe zaobilaženju nedostataka koji bi najverovatnije doveli do problema.

Tokom dvadeset prvog veka, upotreba Markov Chain Monte Carlo uzorkovanja, ili samo MCMC, dramatično je porasla. Ali, šta je zapravo MCMC? I zašto njegova popularnost tako brzo raste? Mi ćemo prikazati osnovu, a ko bude zainteresovan, može se posvetiti detaljnijem proučavanju, jako je interesantno!

Markov Chain Monte Carlo (MCMC) metode obuhvataju klasu algoritama za uzorkovanje iz raspodele verovatnoća, kada je direktno uzorkovanje teško. Dobijeni uzorak se može koristiti za aproksimiranje raspodele (npr. za generisanje histograma) ili za izračunavanje integrala (npr. očekivana vrednost). Konstruisanjem lanca Markova koji ima našu željenu raspodelu kao stacionarnu raspodelu, može se dobiti uzorak iz željene raspodele posmatranjem lanca nakon određenog broja koraka (što je više koraka, raspodela je bliža stvarnoj raspodeli koja nas interesuje).

Podsećanje: Lanac Markova je uređen, indeksiran skup slučajnih veličina (stohastički proces) gde vrednost svake veličine verovatnosno zavisi samo od prethodne vrednosti. Konkretno, ako je \(\{X_0, X_1, X_2, ...\}\) lanac Markova, onda ima Markovsko svojstvo: \[P\{X_t \in \mathcal{A}|X_0,X_1,...,X_{t-1}\} = P\{X_t \in \mathcal{A}|X_{t-1}\}\]

Dakle, \(X_t\) je je uslovno nezavisno od svih ranijih vrednosti osim prethodne. Tako da vrednosti u Markovskom lancu nisu nezavisne, ali jesu “skoro nezavisne”.

MCMC algoritmi se obično koriste za uzorkovanje iz višedimenzionalnih raspodela, posebno kada je broj dimenzija visok. Za jednodimenzione raspodele, obično postoje druge, jednostavnije metode koje direktno vraćaju nezavisne uzorke iz raspodele i one su slobodne od problema korelisanih uzoraka koji su svojstveni MCMC metodama.

MCMC metode se prvenstveno koriste za izračunavanje numeričkih aproksimacija višedimenzionalnih integrala. Ove metode imaju široku primenu. Nama će od najvećeg značaja biti njihova primena u Bajesovskoj statistici, ali one imaju primene i u biologiji, fizici, itd. U Bajesovskoj statistici, razvoj MCMC metoda predstavlja ključni korak u omogućavanju izračunavanja velikih hijerarhijskih modela koji zahtevaju integracije preko velikog broja nepoznatih parametara.

Znamo da u mnogim slučajevima aposteriorna raspodela nema jednostavnu prepoznatljivu formu, tako da iz nje ne možemo uzorkovati korišćenjem ugrađenih R funkcija poput rgamma. U ovom slučaju, kao što smo rekli, koriste se Markov Chain Monte Carlo (MCMC) metode.

Ukratko, cilj Bajesovog zaključivanja je da računamo aposteriornu raspodelu nekog parametra (ili parametara), ako nam je dat nekakav uzorak. Međutim, računanje i izražavanje ove raspodele često uključuje računanje komplikovanih integrala, koji su, za većinu netrivijalnih modela, nerešivi. Algoritmi uzorkovanja bazirani na MCMC (Markov Chain Monte Carlo) metodama su jedan od mogućih načina za izvođenje zaključaka o takvim modelima.

Na primer, možemo proceniti bilo koje očekivanje aposteriorne raspodele dok god imamo N simuliranih uzoraka iz te raspodele kao: \[E[f(s)]_{P}\approx\frac{1}{N}\sum_{i=1}^{N}f(s^{(i)})\] gde je \(P\) aposteriorna raspodela, \(f(s)\) željena veličina čije očekivanje nas interesuje, a \(f(s^{(i)})\), je njena vrednost na \(i\)-tom simuliranom uzorku iz raspodele \(P\). Na primer, možemo da računamo srednju vrednost kao \[E[x]_{P} = \frac{1}{N}\sum_{i=1}^{N}x^{(i)}\]

U praksi, lanac se postepeno kreira tako što se počne od skupa tačaka proizvoljno odabranih i dovoljno udaljenih. Zatim se kreira lanac, tako što se predlažu nova stanja i odlučuje se da li će se ona prihvatiti ili ne. Ovi lanci su stohastički procesi koji se kreću po uzoračkom prostoru na slučajan način, ali tako da je veća verovatnoća da će tačke u regionima velike verovatnoće biti prihvaćene.

Ovi algoritmi stvaraju Markovske lance, takve da imaju ravnotežnu raspodelu koja je proporcionalna datoj funkciji.

Obično nije teško konstruisati Markovski lanac sa željenim svojstvima. Veći problem je odrediti koliko je koraka potrebno za konvergenciju stacionarnoj raspodeli unutar prihvatljive greške. Dobar lanac će imati “brzo mešanje”: stacionarna raspodela se brzo dostiže počevši od proizvoljne pozicije.

Tipično, MCMC može samo aproksimirati ciljnu raspodelu, jer uvek postoji neki rezidualni efekat početne pozicije. Sofisticiraniji MCMC algoritmi (poput coupling from the past) mogu da proizvode i tačne uzorke, ali po cenu dodatnog računanja i neograničenog (iako ograničenog u očekivanju) vremena rada.

Osnovni primeri MCMC metoda, koje ćemo mi koristiti su:

Pored ovih koje ćemo mi pominjati, postoje i mnoge druge metode koje se koriste poput Slice sampling, Reversible-jump, itd.

Metropolis-Hastings algoritam

Metropolis-Hastings algoritam je jedna od osnovnih MCMC metoda za dobijanje slučajnih uzoraka iz raspodela verovatnoća iz kojih je teško direktno uzorkovanje.

Ideja

Recimo da želimo da vadimo uzorke iz raspodele verovatnoća \(f(x)\), međutim, ne možemo to da uradimo jer je direktno uzorkovanje teško. Ono što prvo pada na pamet je da pokušamo nekako da aproksimiramo.

Metropolis-Hastings kaže da ako znamo kako da računamo vrednosti neke funkcije \(\tilde{f}(x)\), koja je proporcionalna raspodeli \(f\), moći ćemo i da izvadimo uzorke iz naše ciljne raspodele. Činjenica da je potrebno da \(\tilde{f}(x)\) bude samo proporcionalna željenoj gustini, ali ne i jednaka u potpunosti nam omogućava ovo, jer je računanje normalizacione konstante obično najveći problem (posebno u Bajesovskom pristupu, kada nam je ciljna raspodela aposteriorna raspodela parametara od interesa).

Zašto ovo funkcioniše?

Pogledajmo specijalan slučaj Metropolis-Hastings algoritma - Metropolis algoritam, kod koga je raspodela predloga simretrična:

  1. Prvo se vrši inicijalizacija - biramo početnu tačku \(x_0\) kao prvu vrednost uzorka.

  2. Biramo raspodelu predloga, odaberemo proizvoljnu raspodelu - nazovimo je \(q(y|x)\) koja predlaže kandidata \(y\) za sledeću vrednost uzorka (za datu prethodnu vrednost uzorka \(x\)). Kod Metropolis algoritma, mora da važi da je \(q(y|x) =q(x|y)\). Čest izbor za funkciju \(q\) je normalna raspodela centrirana oko \(x\), tako da je verovatnije da će tačke blizu \(x\) biti posećene sledeće.

  3. Sada pravimo lanac. Ponavljamo sledeći postupak dok ne popunimo naš uzorak unapred određenog obima (velikog!). Dakle, u \(i\)-toj iteraciji:

  • Generišemo kandidata \(y\) za sledeću vrednost u uzorku, uzorkujući ga iz raspodele \(q(y|x_i)\).

  • Računamo koeficijent prihvatanja \(r = \frac{\tilde{f}(y)}{\tilde{f}(x_i)}\), koji koristimo da odlučimo da li ćemo prihvatiti ili odbaciti kandidata. Dovoljno je da ovako računamo, jer znamo da je \(\tilde{f}\) proporcionalno \(f\), pa je \(r = \frac{\tilde{f}(y)}{\tilde{f}(x_i)}=\frac{f(y)}{f(x_i)}\).

  • Generišemo broj \(u\) iz \(\mathcal{U}[0,1]\) raspodele. Ako je \(u \leq r\) prihvatamo kandidata, i dodajemo predlog kao sledeći element uzorka (\(x_{i+1}=y\)), a ako je \(u \geq r\) odbacujemo kandidata, i kao sledeći element uzorka prenosimo prethodni (\(x_{i+1}=x_{i}\)).

Dakle, algoritam radi tako što u svakoj iteraciji bira kandidata za sledeću vrednost uzorka na osnovu trenutne vrednosti uzorka, pri čemu raspodela sledeće vrednosti uzroka zavisi samo od trenutne vrednosti (čime se niz uzoraka pretvara u Markovski lanac). Zatim se, sa određenom verovatnoćom, kandidat prihvata (u tom slučaju se vrednost kandidata koristi u sledećoj iteraciji), ili odbacuje (u tom slučaju se vrednost kandidata odbacuje, a trenutna vrednost se ponovo upotrebljava u sledećoj iteraciji). Što je više iteracija izvršeno, odnosno, što je više vrednosti u uzorku kreirano, uzoračka raspodela će bolje da aproksimira ciljnu raspodelu.

MCMC pokušava da aproksimira plavu raspodelu narandžastom raspodelom

Prema tome, algoritam se odvija slučajnim kretanjem po uzoračkom prostoru, tako što se predlozi ponekad prihvataju, i skače se na novu vrednost, a ponekad se ostaje na istom mestu.

Koeficijent \(r\) pokazuje koliko je verovatno da je novi predloženi uzorak iz ciljne raspodele, u odnosu na trenutnu vrednost u uzorku. Ako pokušamo da pređemo na tačku koja je verovatnija od postojeće tačke, prihvatićemo predlog. Međutim, ako pokušamo da pređemo na manje verovatnu tačku, ponekad ćemo prihvatiti, a ponekad odbaciti predlog (što je \(r\) manje, veća je verovatnoća da će kandidat biti odbačen). Na taj način nastojimo da ostanemo u regionima velike verovatnoće, dok samo povremeno posećujemo regione male verovatnoće. Intuitivno, ovo je razlog zašto ovaj algoritam funkcioniše, i daje uzorke iz željene raspodele \(f(x)\).

Pojašnjenje

Kao što smo već rekli, ideja koja stoji iza Metropolis-Hastings algoritma je da simulira lanac Markova u prostoru stanja tako da stacionarna raspodela lanca bude ciljna raspodela. Kod tradicionalnih analiza lanaca Markova obično su date verovatnoće prelaza, a nepoznata je stacionarna raspodela, dok je kod MCMC lanaca poznata ravnotežna raspodela, a verovatnoće prelaza se opisuju tako da se ta ravnoteža postigne.

U MCMC algoritmima, uobičajeno je da se koriste reverzibilni lanci Markova, jer važenje detailed balance uslova za željenu raspodelu implicira da je lanac Markova konstruisan tako da mu je ta raspodela stacionarna raspodela. (Lanac Markova je reverzibilan ako zadovoljava detailed balance jednačinu, tj. ako \[\pi_i \cdot p_{ij}=\pi_j \cdot p_{ji}\] gde je \(p_{ij} = P(X_n=j|X_{n-1}=i)\) verovatnoća prelaza iz stanja \(i\) na stanje \(j\), a \(\pi_i\) i \(\pi_j\) su verovatnoće da lanac u ravnotežnom stanju bude u stanju \(i\), odnosno \(j\).)

Prema tome, Metropolis-Hastings algoritam uključuje kreiranje Markovskog procesa (konstruisanjem verovatnoća prelaza) koji zadovoljava gornji uslov, tako da je njegova stacionarna raspodela \(\pi\) baš naša željena raspodela \(f\).

Počinjemo prvim uslovom - detailed balance. Dakle, treba da važi da je: \[\pi(x) \cdot p_{xy} = \pi(y) \cdot p_{yx}\] odnsono \[f(x) \cdot p_{xy} = f(y) \cdot p_{yx}\] Dalje, dobijamo \[\frac{ f(y)}{f(x)} = \frac{p_{xy}}{p_{yx}}\]

Sada verovatnoće prelaza sa stanja \(x\) na stanje \(y\) i obrnuto, možemo razdvojiti na dva dela: predlog novog stanja i verovatnoća prihvatanja ili odbijanja, tj. \[p_{xy} = q(y|x) \cdot r(x \rightarrow y)\]

gde je \(q(y|x)\) gore pomenuta raspodela predloga, odnosno verovatnoća predlaganja \(y\) za dato trenutno stanje \(x\), a \(r(x \rightarrow y)\) je verovatnoća prihvatanja predloženog stanja \(y\).

Prema tome, dobijamo: \[\frac{r(x \rightarrow y)}{r(y \rightarrow x)} = \frac{p_{xy}}{q(y|x)} \cdot \frac{q(x|y)}{p_{yx}} = \frac{f(y)q(x|y)}{f(x)q(y|x)}\] Sada još treba da odredimo verovatnoću prihvatanja \(r\) takvu da zadovoljava navedene uslove. Čest izbor je:

\[r(x \rightarrow y) = min \big(1, \frac{f(y)q(x|y)}{f(x)q(y|x)} \big)\] Za ovako odabrano \(r\), vidimo da mora biti \(r(x \rightarrow y)=1\) ili \(r(y \rightarrow x)=1\).

Ako koristimo Metropolis algoritam (\(q(y|x)=q(x|y)\)), \[r(x \rightarrow y) = min \big(1, \frac{f(y)}{f(x)} \big)\]

Sada sprovodimo algoritam na već opisani način:

ALGORITAM:

  1. Izaberemo početno stanje \(x_0\).

  2. U svakoj iteraciji predlažemo nove vrednosti i prihvatamo ih ili odbacujemo u zavisnosti od \(r\). To jest, u \(i\)-toj iteraciji:

  • Generišemo kandidata \(y\) na osnovu raspodele predloga \(q(y|x_i)\);

  • Računamo verovatnoću prihvatanja \(r(x \rightarrow y) = min \big(1, \frac{f(y)q(x|y)}{f(x)q(y|x)} \big)\);

  • Generišemo slučajan broj \(u\) iz \(\mathcal{U}[0,1]\) raspodele - ako je \(u \leq r(x_i \rightarrow y)\) prihvatamo novo stanje i postavljamo \(x_{i+1}=y\), a ako je \(u \geq r(x_i \rightarrow y)\) odbacujemo novo stanje i postavljamo \(x_{i+1}=x_{i}\);

Ponavljamo prethodne korake \(N\) puta, gde je \(N\) veliki broj. Raspodela uzoračkih vrednosti \(x_0,...,x_N\) će aproksimirati našu željenu raspodelu. Koliko treba da bude \(N\), odnosno koliko je iteracija potrebno da se efektivno procenila željena raspodlela razlikuje se od slučaja do slučaja i zavisi od nekoliko faktora, uključujući vezu između željene raspodele i raspodele predloga i željene tačnosti procene.

Nedostaci

Postoji niz problema koji se mogu pojaviti u radu sa MCMC algoritmima.

Pre svega, jedan od najvažnijih problema je što su vrednosti u uzorku dobijene na ovaj način korelisane. Dakle, iako posle velikog broja iteracija dobijamo približno uzorak iz ciljne raspodele, bliske vrednosti su korelisane, i neće ispravno odražavati raspodelu.

Dodatno, iako se lanac na kraju približava ciljnoj raspodeli, početne vrednosti obično ne potiču iz nje, pa je potrebno odbaciti određeni broj početnih vrednosti. Dakle, Markovski lanac se pokreće iz proizvoljne tačke i algoritam se ponavlja za mnogo iteracija sve dok ovo početno stanje ne bude “zaboravljeno”. Ovi uzorci, koji se odbacuju, poznati su kao burn-in. Preostali skup prihvaćenih vrednosti predstavlja uzorak iz željene raspodele.

Kao što smo već pomenuli, u opštem slučaju ne znamo koji je broj iteracija potreban za pravilnu procenu. Takođe, ne znamo ni koju raspodelu predloga odbrati u opštem slučaju - i jedno i drugo se prilagođavaju konkretnom problemu. Pokazalo se da algoritam radi najbolje ako raspodela predloga odgovara obliku željene raspodele tj. ako je \(q(y|x_t) \approx f(y)\).

Dodatno, ako se koristi normalna raspodela kao raspodela predloga, parametar disperzije mora biti podešen tokom burn-in perioda. Ako je odabrana disperzija mala, lanac će se sporo mešati (stopa prihvatanja će biti visoka, ali će se uzastopne vrednosti kretati sporo po uzoračkom prostoru i lanac će sporo konvergirati ka željenoj raspodeli). S druge strane, ako je disperzija prevelika, stopa prihvatanja će biti veoma niska, jer je verovatnije da predlozi leže u regionima manje verovatnoće, pa će opet lanac konvergirati veoma sporo. Obično se podešava diperzija raspodele predloga tako da algoritam prihvata približno \(30\%\) svih uzoraka, ali opet, zavisi od konkretnog problema.(Obično se širina skoka određuje izračunavanjem stope prihvatanja, koja predstavlja udeo predloženih vrednosti koje su prihvaćene. Željena stopa prihvatanja zavisi od ciljne raspodele, međutim, teorijski je pokazano da je idealna stopa prihvatanja za jednodimenzionu normalnu raspodelu približno \(50\%\), smanjujući se do približno \(23\%\) normalnu ciljnu raspodelu većih dimenzija.)

Vratimo se na prvi problem - bliske vrednosti su korelisane, a mi obično želimo skup nezavisnih vrednosti. Ovo se može ublažiti tako što ćemo odbaciti većinu vrednosti i uzeti svaku \(k\)-tu vrednost u uzorak, za neko \(k\), koje se obično određuje ispitivanjem autokorelacije između susednih vrednosti. Autokorelacija se može smanjiti povećanjem širine skoka (koja je povezana sa disperzijom raspodele skoka, odnosno raspodele predloga), međutim, to će povećati i verovatnoću odbacivanja predloga. Kao što smo već rekli, prevelika ili premala širina skoka će dovesti do sporog mešanja, tako da će biti potreban veliki broj iteracija da bi se dobila razumna procena željene raspodele.

Međutim, i uz sve probleme i nedostatke, MCMC metode se jako često koriste i jako su moćne. Jedna od glavnih prednosti je to što, dok većina jednostavnih metoda uzorkovanja, koje daju nezavisne uzorke, prati od problema dimenzionalnosti, Metropolis-Hastings i slične metode, nemaju taj problem u tolikoj meri, pa su stoga često jedino rešenje kada želimo da uzorkujemo iz raspodela velikih dimenzija, na primer iz hijerarhijskih Bajesovih modela i drugih visokodimenzionalnih statističkih modela koji se danas koriste u mnogim oblastima.

Kada uzorkujemo iz višedimenzionih raspodela, klasični Metropolis-Hastings algoritam u svakom koraku uključuje izbor nove višedimenzionalne tačke uzorka. Međutim, kada je broj dimenzija veliki, pronalaženje odgovarajuće raspodele predloga može da bude teško, jer širina skoka mora biti dobro odabrana za sve komponente odjednom da bi se izbeglo previše sporo mešanje. To je često teško uraditi, pa se često prelazi na Gibbs-ovo uzorkovanje, koje podrazumeva uzorkovanje iz svake komponente posebno, umesto za sve odjednom. Dakle, promenljive se uzorkuju jedna po jedna, pri čemu je svaka promenljiva uslovljena najskorijim vrednostima svih ostalih. za izbor ovih pojedinačnih uzoraka , mogu se koristiti razni algoritmi, u zavisnosti od tačnog oblika višedimenzione raspodele (npr. slice sampling, adaptive rejection Metropolis sampling, itd.).

Ovo je posebno primenljivo kada imamo skup slučajnih promenljivih gde je svaka promenljiva uslovljena samo malim brojem drugih promenljivih, kao što je slučaj u većini hijerarhijskih modela.

Primer 1

Ilustrovaćemo Metropolis algoritam na sledećem primeru. Vadimo uzorak iz normalne raspodele sa srednjom vrednošću 0 i jediničnom disperzijom koristeći Metropolis algoritam sa uniformnom raspodelom predloga.

Lanac je inicijalizovan na 0, i u svakoj iteraciji se predlaže nova vrednost uz pomoć \(\mathcal{U}(-\alpha,\alpha)\) raspodele.

norm = function (n, alpha) {
        vec = rep(0, n)
        x = 0
        vec[1] = x
        for (i in 2:n) {
                innov = runif(1, -alpha, alpha)
                predlog = x + innov
                r = min(1, dnorm(predlog)/dnorm(x))
                u = runif(1)
                if (u < r) 
                        x = predlog
                vec[i] = x
        }
        return(vec)
}

Dakle, predlog je kandidat za novi element uzorka, a r je verovatnoća prihvatanja. Odluka o prihvatanju ili odbacivanju se onda donosi na osnovu toga da li je slučajno odbran broj između 0 i 1 manji od verovatnoće prihvatanja. Možemo testirati na sledeći način:

vec = norm(10000, 1)
par(mfrow = c(2,1))
plot(ts(vec))
hist(vec, 30)

Ovo nam prikazuje brzo mešajući lanac i približno normalnu raspodelu vrednosti. Međutim, ovo je zasnovano na izboru \(\alpha\) koje je jednako 1. Druge vrednosti \(\alpha\) neće uticati na stacionarnu raspodelu, ali će uticati na brzinu mešanja lanca.

Primer 2

Pokušajmo sada da simuliramo iz gama raspodele sa proizvoljnim parametrima, koristeći Metropolis-Hastings algoritam sa normalnom raspodelom predloga, sa istom srednjom vrednošću i disperzijom kao i željena gama raspodela.

Lanac inicijalizujemo na \(\frac{a}{b}\), i u svakom koraku predlažemo kandidata iz \(\mathcal{N}(\frac{a}{b}, \frac{a}{b^2})\) raspodele.

gamm = function (n, a, b) {
        mu = a/b
        sig = sqrt(a/(b^2))
        vec = rep(0, n)
        x = a/b
        vec[1] = x
        for (i in 2:n) {
                predlog =  rnorm(1, mu, sig)
                r = min(1, (dgamma(predlog, a, b)*dnorm(x, mu, sig))/(dgamma(x, a, b)*dnorm(predlog, mu, sig)))
                u = runif(1)
                if (u < r) 
                        x = predlog
                vec[i] = x
        }
        return(vec)
}

Dakle, predlog je kandidat za novi element uzorka, a r je verovatnoća prihvatanja. Odluka o prihvatanju ili odbacivanju se onda donosi na osnovu toga da li je slučajno odbran broj između 0 i 1 manji od verovatnoće prihvatanja. Možemo testirati na sledeći način:

vec = gamm(10000, 2.3, 2.7)
par(mfrow = c(2,1))
plot(ts(vec))
hist(vec, 30)

Ovo nam prikazuje brzo mešajući lanac i približno gama raspodelu vrednosti.

Primer 3

Jedna od najkorisnijih primena MCMC metoda je u Bajesovskom pristupu za dobijanje uzoraka iz aposteriorne raspodele, jer nju obično nije lako odrediti analitički.

MCMC generiše uzorke iz aposteriorne raspodele konstruisanjem reverzibilnog lanca Markova koji za ravnotežnu raspodelu ima našu ciljnu aposteriornu raspodelu.

Prisetimo se Bajesove formule: \[P(\theta|X) = \frac{P(X|\theta)P(\theta)}{P(X)}\] Imamo \(P(\theta|X)\), verovatnoću naših parametara modela \(\theta\) za date podatke \(X\), što je naša raspodela od interesa. Da bismo je izračunali, pomnožimo apriornu \(P(\theta)\) (naša uverenja o parametru \(\theta\) pre nego što smo videli bilo koji podatak) i verovatnoću \(P(X|\theta)\) tj. kako su naši podaci raspodeljeni. Ovaj brojilac je obično prilično lako izračunati.

Pogledajmo malo bolje imenilac: \(P(X)\). Prisetimo se, ova vrednost se dobija integracijom po svim mogućim vrednostima parametra \(\theta\): \[P(X)=\int_{\Theta} P(X, \theta)P(\theta)d\theta\]

Ovo je glavna poteškoća u Bajesovoj formuli - formula izgleda naivno, ali čak i za malo komplikovanije modele od onih trivijalnih, ne možemo izračunati aposteriornu raspodelu ekspilicitno.

Sada bismo mogli reći - U redu, ne možemo da je izračunamo, ali možemo li da je aproksimiramo? Na primer, ako bismo nekako mogli izvući uzorke iz te aposteriorne raspodele, možemo je aproksimirati Monte Carlo metodom. Nažalost, da bismo direktno izvukli uzorak iz te raspodele, moramo da rešimo Bajesovu formulu.

Onda možemo reći - “Dobro, hajde onda da konstruišemo ergodični, reverzibilni lanac Markova koji ima stacionarnu raspodelu koja odgovara našoj aposteriornoj raspodeli”. Ako ne možemo izračunati raspodelu, ne možemo iz nje izvući uzorak, onda konstrisanje Markovskog lanca sa svim tim svojstvima deluje još teže. Međutim, uz pomoć MCMC-a, ovo je veoma lako postići.

Postavka problema

Hajde prvo da generišemo neke podatke: 20 tačaka iz normalne raspodele centrirane oko 0. Naš cilj je da izračunamo aposterironu raspodelu parametra \(\mu\) (pretpostavićemo, radi jednostavnosti, da je standardna devijacija jednaka 1).

data = rnorm(20)
hist(data)

Zatim moramo definisati naš model. U ovom jednostavnom primeru, pretpostavićemo da su naši podaci normalno raspodeljeni. Kao što znamo, normalna raspodela ima dva parametra - očekivanje \(\mu\) i disperziju \(\sigma^2\). Radi jednostavnosti, pretpostavićemo da je \(\sigma^2=1\) i da želimo da odredimo aposteriornu raspodelu za \(\mu\). Za svaki parametar o kome želimo da zaključujemo, moramo izabrati apriornu raspodelu. Radi jednostavnosti, pretpostavimo normalnu raspodelu i za \(\mu\). Dakle, naš statistički model je: \[\mu \sim \mathcal{N}(0,1)\] \[X|\mu \sim \mathcal{N}(\mu,1)\]

Ono što je zgodno za ovaj model je to što možemo da izračunamo analitički aposteriornu raspodelu. To je zato što za normalnu raspodelu sa poznatom disperzijom, normalna apriorna za \(\mu\) je konjugovana (aposteriorna će imati isti oblik kao apriorna), tako da znamo da je naša aposteriorna raspodela za \(\mu\) takođe normalna. Parametre aposteriorne raspodele možemo lako izračunati (Podsećanje).

aposteriorna = function(x, data) {
  n = length(data)
  mu = n*mean(data)/(n+1)
  sigma = sqrt(1/(n+1))
  return(dnorm(x, mu, sigma))
}

mu = seq(-1, 1, length.out = 500)
plot(mu, aposteriorna(mu, data), type = "l", ylab = "Analiticka aposteriorna")

Ovo nam prikazuje našu raspodelu od interesa, vrednosti \(\mu\) nakon što smo videli podatke, uzimajući u obzir naša prethodna uverenja. Pretpostavimo, međutim, da naše raspodele nisu bile konjugovane i da nismo mogli da rešimo ovo ovako lako, što je obično i slučaj.

Pređimo na logiku uzorkovanja. Za početak, biramo početnu poziciju parametra, možemo je odabrati proizvoljno, pa hajde da je fiksiramo na:

mu_current = 1

Zatim, predlažemo da se pomerimo (tj. skočimo) sa te pozicije na neku drugu (to je Markovski deo). Kao predlog uzećemo uzorak iz normalne raspodele (nema veze sa normalnom raspodelom pretpostavljenom za model), centrirane oko naše trenutne vrednosti za \(\mu\), sa određenom disperzijom (proposal_width) koja će odrediti koliko daleko predlažemo skokove.

proposal = rnorm(1, mu_current, proposal_width)

Nakon toga, procenjujemo da li je dobro napraviti skok ili ne. Ako rezultujuća normalna raspodela sa predloženim \(\mu\) bolje opisuje podatke od starog \(\mu\), definitivno ćemo želeti da se prebacimo na tu vrednost. Međutim, kako da odredimo šta bolje opisuje podatke? Kvantifikujemo uklapanje izračunavajući verovatnoću podataka (likelihood), sa predloženim vrednostima parametara (predloženo \(\mu\) i \(\sigma=1\)). Ovo se lako može izračunati izračunavanjem verovatnoće za svaku tačku podataka, a onda množenjem pojedinačnih verovatnoća, tj. računanjem funkcije verodostojnosti (obično koristimo logaritme verovatnoća, ali ovde nećemo tako raditi):

likelihood_current = prod(dnorm(data, mu_current, 1))
likelihood_proposal = prod(dnorm(data, mu_proposal, 1))

# Racunamo apriorne gustine predlozenog i trenutnog mu
prior_current = dnorm(mu_current, mu_prior_mu, mu_prior_sd)
prior_proposal = dnorm(mu_proposal, mu_prior_mu, mu_prior_sd)

# Brojilac u Bajesovoj formuli
p_current = likelihood_current*prior_current
p_proposal = likelihood_proposal*prior_proposal

Računamo koeficijent prihvatanja:

p_accept = p_proposal/p_current

Vidimo da ako je verovatnoća predloga veća, ova verovatnoća će biti veća od 1, pa ćemo sigurno prihvatiti predlog. Međutim, ako je verovatnoća predloga manja, onda ćemo možda skočiti na novu poziciju, a možda ostati na staroj, što zavisi od koeficijenta prihvatanja:

accept = runif(1) < p_accept

if(accept) 
  # Pomeramo se na novu poziciju
  cur_pos = proposal

Ova jednostavna procedura nam daje uzorke iz aposteriorne raspodele.

Zašto ovo ima smisla?

Ako se vratimo korak unazad, vidimo da je gore navedeni koeficijent prihvatanja razlog zašto ovo radi i zašto možemo da zaobiđemo integraciju. To možemo pokazati izračunavanjem koeficijenta prihvatanja za normalizovanu aposteriornu raspodelu. Vidimo kako je on ekvivalentan sa koeficijentom prihvatanja za nenormalizovanu aposteriornu raspodelu (recimo da je \(\mu_0\) naša trenutna pozicija, a \(\mu\) naš predlog): \[\frac{\frac{P(x|\mu)P(\mu)}{P(x)}}{\frac{P(x|\mu_0)P(\mu_0)}{P(x)}} = \frac{P(x|\mu)P(\mu)}{P(x|\mu_0)P(\mu_0)}\] Deljenjem aposteriorne raspodele predloženog parametra aposteriornom raspodelom trenutnog, \(P(X)\) - gadni integral koji ne možemo da izračunamo, biva izbačen. Prema tome, možemo reći da smo zapravo podelili pravu aposteriornu raspodelu u jednoj poziciji sa pravom aposteriornom u drugoj (dakle, nema nikakve magije ovde). Na taj način, relativno češće posećujemo regije visoke aposteriorne verovatnoće od onih sa manjom aposteriornom verovatnoćom.

Dakle, kada sve spojimo u jedan kod:

sampler = function(data, samples = 4, mu_init=5, proposal_width=5, mu_prior_mu=0, mu_prior_sd=1) {

mu_current = mu_init

posterior = c()
posterior[1] = mu_current

for (i in 1:samples) {
# predlazemo novu poziciju
mu_proposal = rnorm(1, mu_current, proposal_width)

# racunamo likelihood mnozenjem verovatnoca za svaku tacku uzorka
likelihood_current = prod(dnorm(data, mu_current, 1))
likelihood_proposal = prod(dnorm(data, mu_proposal, 1))

# racunamo apriorne raspodele trenutnog i predlozenog mu
prior_current = dnorm(mu_current, mu_prior_mu, mu_prior_sd)
prior_proposal = dnorm(mu_proposal, mu_prior_mu, mu_prior_sd)

p_current = likelihood_current * prior_current
p_proposal = likelihood_proposal * prior_proposal

# verovatnoca prihvatanja
p_accept = p_proposal / p_current

accept = runif(1) < p_accept
if(accept) {
  mu_current = mu_proposal
}

posterior = c(posterior, mu_current)
i = i+1
}
return(posterior)
}

Da bismo videli kako se vrši uzorkovanje, prikazaćemo neke uobičajene grafike. Svaki red ispod predstavlja jednu iteraciju u našem Metropolis sampler-u.

Prva kolona predstavlja apriornu raspodelu - to je naše uverenje o \(\mu\) pre nego što vidimo podatke. Možemo da vidimo da je raspodela statična i da samo uključujemo \(\mu\) predloge. Vertikalne linije predstavljaju trenutno \(\mu\) u plavoj, a predloženo \(\mu\) u zelenoj ili crvenoj (prihvaćeno ili odbijeno).

Druga kolona je funkcija verodostojnosti i ono što koristimo da procenimo koliko dobro naš model objašnjava podatke. Možemo videti da se likelihood menja u zavisnosti od predloženog \(\mu\). Plavi histogram predstavlja podatke. Puna linija u zelenoj ili crvenoj predstavlja likelihood sa trenutno predloženim \(\mu\). Intuitivno, što je više preklapanja između likelihood-a i podataka, bolje model objašnjava podatke i veća će biti rezultujuća verovatnoća. Isprekidana linija iste boje je predloženo \(\mu\) , a isprekidana linija plave boje je trenutno \(\mu\).

Treća kolona predstavlja aposteriornu raspodelu. Ovde prikazujemo normalizovanu aposteriornu, ali kako smo videli gore, možemo samo da pomnožimo vrednosti apriorne raspodele sa vrednostima likelihooda za predloženo i trenutno \(\mu\), da bismo dobili nenormalizovane vrednosti aposteriorne raspodele (koje zapravo koristimo u računanju), i podelimo jednu drugom da bismo dobili našu verovatnoću prihvatanja.

Četvrta kolona predstavlja trag (tj. aposteriorne uzorke \(\mu\) koje generišemo) gde prikazujemo da li je predlog prihvaćen ili odbačen (u tom slučaju linija samo ostaje konstantna).

Treba imati na umu da mi uvek prelazimo na relativno verovatnije \(\mu\) vrednosti (u smislu njihove aposteriorne verovatnoće), a samo ponekad na relativno manje verovatne vrednosti \(\mu\).

Sada je magija MCMC-a u tome što samo moramo da radimo ovo dugo vremena, i uzorci koje smo generisali na ovaj način biće približno iz aposteriorne raspodele našeg modela.

Da bismo stekli osećaj šta je rezultat ovoga, izvucimo puno uzoraka i nacrtajmo ih:

posterior = sampler(data, samples = 15000, mu_init = 1)
plot(ts(posterior), xlab = "sample", ylab = "mu")

Da bismo dobili aproksimaciju aposteriorne raspodele (što je i razlog zbog kojeg radimo sve ovo), jednostavno nacrtamo histogram traga.

Važno je napomenuti da iako ovo izgleda slično podacima koje smo uzorkovali gore da bismo napravili model, ta dva su potpuno različita. Donji grafik predstavalja naše verovanje o \(\mu\). U ovom slučaju, zadesilo se da je opet normalna raspodela, međutim, za neki drugi model, može imati potpuno drugačiji oblik od likelihood-a ili apriorne raspodele.

hist(posterior[500:length(posterior)], freq = FALSE)
x = seq(-5, 5, length.out = 500)
lines(x, aposteriorna(x,data), type = "l", xlim = c(-0.5,0.5), col = "red")

Kao što možemo videti, sledeći gornju proceduru, dobijamo uzorke iz iste raspodele, kao iz one koju smo analitički izveli.

Širina predloga

Gore smo postavili širinu predloga na 0.5. Pokazalo se da je to prilično dobra vrednost. U opštem slučaju, ne želimo da širina bude suviše mala, jer će uzorkovanje biti neefikasno zato što je potrebno puno vremena za istraživanje celog parametarskog prostora i uzorkovanje pokazuje tipično ponašanje slučajnog lutanja:

posterior_small = sampler(data, samples = 10000, mu_init = 1, proposal_width = 0.01)
plot(ts(posterior_small), xlab = "sample", ylab = "mu")

Međutim, ne želimo da bude ni prevelika pa da nikada ne prihvatamo predlog:

posterior_large = sampler(data, samples = 10000, mu_init = 1, proposal_width = 3)
plot(ts(posterior_large), xlab = "sample", ylab = "mu")

Primetimo da i dalje uzorkujemo iz naše ciljne raspodele garantovane matematičkim dokazom, samo manje efikasno:

hist(posterior_small[1001:length(posterior_small)], xlim = c(-0.5, 0.5), freq = FALSE)
lines(x, aposteriorna(x,data), type = "l", xlim = c(-0.5,0.5), col = "red")

hist(posterior_large[1001:length(posterior_large)], freq = FALSE)
lines(x, aposteriorna(x,data), type = "l", xlim = c(-0.5,0.5), col = "red")

Sa više uzoraka, ovo će na kraju izgledati kao prava aposteriorna raspodela.

Proširenje na složenije modele

Sada možemo lako zamisliti da možemo dodati i \(\sigma\) parametar za standardno odstupanje i slediti istu proceduru za ovaj drugi parametar. U tom slučaju, mi bismo generisali predloge za \(\mu\) i \(\sigma\), ali logika algoritma bi bila skoro identična. Ili, mogli bismo imati podatke iz veoma različite raspodele kao što je Binomna i opet bismo koristili isti algoritam i dobili ispravnu aposteriornu. To je prilično zgodno i velika prednost probabilističkog programiranja: Samo definišemo model koji želimo i pustimo MCMC da se pobrine za zaključak.