Gibbs-ovo uzorkovanje je jedna od Markov Chain Monte Carlo (MCMC) metoda za dobijanje uzoraka koji su aproksimativno iz višedimenzione raspodele verovatnoće, kada je direktno uzorkovanje iz nje teško. To je MCMC algoritam koji aproksimira zajedničku raspodelu više slučajnih veličina uzimanjem uzoraka iz uslovne raspodele svake promenljive.

Gibbs-ovo uzorkovanje je primenljivo kada zajednička raspodela nije eksplicitno poznata ili je teško uzorkovati iz nje direktno, ali je uslovna raspodela svake promenljive poznata i lako je (ili barem lakše) iz nje uzorkovati. Glavna ideja je da se razbije problem uzorkovanja iz višedimenzione raspodele na serije uzoraka iz uslovnih raspodela manjih dimenzija. Gibbs-ov algoritam za uzorkovanje generiše vrednost iz raspodele svake promenljive pod uslovom trenutnih vrednosti drugih promenljivih. Kao i kod drugih MCMC algoritama, Gibbs-ovim uzorkovanjem se generiše Markovski lanac, td. je svaka uzorkovana vrednost povezana sa bliskim vrednostima (zato se mora paziti ako su potrebni nezavisni uzorci!). Dakle, može se pokazati da niz uzoraka čini Markovski lanac, a stacionarna raspodela tog Markovskog lanca je tražena zajednička raspodela.

Niz uzoraka dobijenih ovom metodom se može koristiti za aproksimaciju zajedničke raspodele (npr. za generisanje histograma raspodele), za aproksimiranje marginalne raspodele neke promenljive, ili skupa promenljivih, ili računanje integrala (npr. očekivanje neke promenljive). Obično neke od promenljivih odgovaraju opservacijama čije su vrednosti poznate, pa ne moraju biti uzorkovane.

Uobičajeno je da uzorkovane vrednosti sa početka lanca (burn-in period) odbacuju, jer verovatno neće tačno predstavljati željenu raspodelu. Ipak, pokazano je da upotreba dužeg lanca dovodi do dobrih procena prave željene raspodele.

Mi ćemo se baviti najosnovnijim tipom Gibbs-ovog uzorkovanja, ali postoje mnoga njegova poboljšanja i proširenja, i ovaj algoritam se, u proširenim verzijama može se smatrati opštim načinom za uzorkovanje iz velikog skupa promenljivih uzorkovanjem svake promenljive (ili u nekim slučajevima, grupa promenljivih), a takođe, može se i ugraditi u Metropolis-Hastings algoritam, ili neke druge metode poput slice sampling, itd.

Gibbs-ovo uzorkovanje se obično koristi kao pomoć kod statističkog zaključivanja, a njegova možda najvažnija primena odnosi se na procenjivanje aposteriorne raspodele kod Bajesovog zaključivanja.

Gibbs-ovo uzorkovanje je predloženo ranih 90-ih i fundamentalno je izmenilo Bajesovsko izračunavanje.

Algoritam:

Ideja kod Gibbs-ovog uzorkovanja je generisanje uzoraka, tako što se prolazi kroz svaku promenljivu i uzorkovanje vrši iz njene uslovne raspodele sa preostalim promenljivama fiksiranim na neke vrednosti.

Na primer, razmotrimo slučajne veličine \(X_1, X_2\) i \(X_3\).

Počinjemo tako što ćemo ih sve postaviti na neke početne vrednosti \(x^{(0)}_1, x^{(0)}_2\) i \(x^{(0)}_3\) (često su to vrednosti dobijene uzorkovanjem iz apriorne raspodele).

U \(i\)-toj iteraciji, uzorkujemo \(x^{(i)}_1 \sim p(x_1|X_2=x^{(i-1)}_2, X_3=x^{(i-1)}_3)\), pa onda \(x^{(i)}_2 \sim p(x_2|X_1=x^{(i)}_1, X_3=x^{(i-1)}_3)\) i na kraju \(x^{(i)}_3 \sim p(x_3|X_1=x^{(i)}_1, X_2=x^{(i)}_2)\).

Ovaj proces se nastavlja do konvergencije (dok vrednosti iz uzorka ne budu imale istu raspodelu kao da su uzorkovane iz prave zajedničke raspodele). Sledeći algoritam prikazuje korake u Gibbs-ovom uzorkovanju (takozvani Gibbs-ov algoritam):

Gibbs-ov algoritam:

Pretpostavimo da želimo da dobijemo \(N\) uzoraka iz zajedničke raspodele slučajnog vektora \(\mathbf{\theta} = (\theta_1,\theta_2,...,\theta_k)\). Algoritam se odvija na sledeći način:

  1. Na početku je potrebno inicijalizovati lanac. Dakle, počinjemo sa nekom inicijalnom vrednošću \(\mathbf{\theta^{(1)}}=(\theta_1^{(1)},\theta_2^{(1)},...,\theta_k^{(1)})\). Zatim, u svakoj iteraciji primenjujemo sledeće:

  2. Pretpostavimo da imamo uzorak \(\mathbf{\theta^{(i)}}\), i želimo sledeći uzorak, \(\mathbf{\theta^{(i+1)}}=(\theta_1^{(i+1)},\theta_2^{(i+1)},...,\theta_n^{(i+1)})\). Pošto je \(\mathbf{\theta^{(i+1)}}\) vektor, potrebno je da uzorkujemo svaku komponentu vektora. To radimo na sledeći način:

  • Prođemo kroz svaku uslovnu raspodelu, uzorkovanjem, \(i=1,2,..\) \[\theta^{(i+1)}_1 \sim p(\theta_1|\theta^{(i)}_2, \theta^{(i)}_3,...,\theta^{(i)}_k)\] \[\theta^{(i+1)}_2 \sim p(\theta_2|\theta^{(i+1)}_1, \theta^{(i)}_3,...,\theta^{(i)}_k)\] \[\cdot\] \[\cdot\] \[\cdot\] \[\theta^{(i+1)}_k \sim p(\theta_k|\theta^{(i+1)}_1, \theta^{(i+1)}_2,...,\theta^{(i+1)}_{k-1})\]

Primetimo da se u svakom koraku koristi najskorija vrednost \(\theta_j\). Dakle, koristimo do sada uzorkovane komponente vektora \(\mathbf{\theta^{(i+1)}}\) , a one koje još uvek nisu uzorkovane prenosimo iz vektora \(\mathbf{\theta^{(i)}}\), odnosno od prve do \((j-1)\)-e komponente koristimo vrednosti koje će biti u trenutnom uzorku, a od \((j+1)\)-ve komponente do poslednje koristimo vrednosti koje su te komponente imale u \(i\)-tom uzorku, a ne u \((i+1)\)-om (njih tek treba da dobijemo). Da bismo ovo postigli, uzorkujemo komponente po redu, počevši od prve komponente.

  1. Ponavljamo korak (2) \(N\) puta.

Dakle, moramo biti u stanju da uzorkujemo iz svake uslovne raspodele da bismo mogli da koristimo Gibbs Sampler. Olakšavajuća okolnost je ta što je uslovna raspodela jedne promenljive (u odnosu na sve ostale) proporcionalna zajedničkoj raspodeli, jer imenilac ne zavisi od nje, pa je uvek isti (bez obzira na to koju vrednost uzme ta slučajna veličina), odnosno on čini deo normalizacione konstante za raspodelu. Npr., za \(\theta_j\) važi da je: \[p(\theta_j|\theta_1,...,\theta_{j-1},\theta_{j+1},...,\theta_k)=\frac{p(\theta_1,...,\theta_k)}{p(\theta_1,...,\theta_{j-1},\theta_{j+1},\theta_k)} \propto p(\theta_1,...,\theta_k)\]

Dakle, da bi se odredio oblik uslovne raspodele za pojedinačnu komponentu, najjednostavnije je zajedničku raspodelu razložiti na pojedinačne uslovne raspodele i ignorisati sve faktore koji ne zavise od te komponente, a zatim ponovo dodati konstantu normalizacije na kraju, po potrebi (ako je raspodela diskretna, izračunaju se pojedinačne verovatnoće svih mogućih vrednosti koje može uzeti ta slučajna veličina, a zatim se sabiraju kako bi se utvrdila konstanta normalizacije, ako je raspodela neprekidna i poznatog oblika, konstanta normalizacije se takođe lako izračunava, a u drugim slučajevima, konstanta normalizacije se obično može ignorisati, jer većina metoda uzorkovanja to ne zahteva).

Nakon puno ciklusa (tj. za veliko \(N\)), uzorkovane vrednosti će aproksimirati slučajne uzorke iz zajedničke raspodele \((\theta_1,...,\theta_k)\).

Dakle, u ovom algoritmu ne uzimamo uzorke iz same aposteriorne raspodele. Umesto toga, simuliramo uzorke tako što prolazimo kroz sve aposteriorne uslovne raspodele, za svaku od slučajnih veličina, kada su preostale fiksirane. Početne vrednosti promenljivih mogu se odrediti nasumično ili uz pomoć nekog algoritma koji daje bržu konvergenciju, a nije potrebno odrediti početnu vrednost za prvu promenljivu. Budući da na početku inicijalizujemo slučajnim početnim vrednostima, uzorci simulirani pri početnim iteracijama ne moraju biti reprezentativni za pravu aposteriornu raspodelu. Međutim, teorija MCMC metoda garantuje da će stacionarna raspodela uzoraka generisanih ovim algoritmom biti ciljna aposteriorna raspodela.

Iz tog razloga, MCMC algoritmi se obično pokreću za veliki broj iteracija (u nadi da će se konvergencija ka ciljnoj aposteriornoj raspodeli postići). Pošto uzorci iz početnih iteracija nisu bliski onima iz aposteriorne raspodele, praksa je, kao i kod drugih MCMC metoda, da se ignoriše određeni broj uzoraka na početku (burn-in). Takođe, često je da se uzme u obzir samo svaki \(n\)-ti kada se, recimo, usrednjavaju vrednosti da bi se izračunalo očekivanje. Na primer, prvih 1000 uzoraka može biti ignorisano, a onda je svaki 100. uzorak uzet u prosek, dok se ostalo odbacuje. Razlog za to je to što: (1) stacionarna raspodela Markovskog lanca je željena zajednička raspodela promenljivih, ali može potrajati dok se ta stacionarna raspodela ne dostigne; (2) uzastopni uzorci nisu nezavisni jedan od drugog, već formiraju lanac Markova sa nekom količinom korelacije. Ponekad se mogu koristiti neki algoritmi za određivanje količine autokorelacije između uzoraka i vrednosti \(n\) (period između uzoraka koji se stvarno koriste), ali u praksi je obično to produkt “crne magije” (kod koji se zasniva na tehnikama koje deluju, ali koje nemaju teorijsko objašnjenje).

Primer 1

Recimo da želimo da izvadimo uzorak iz dvodimenzionalne normalne raspodele sa očekivanjem nula, i kovarijacionom matricom sa jedinicama na glavnoj dijagonali i \(\rho\) na sporednoj. Naravno, da bismo simulirali uzorak iz ove raspodele, ne moramo da koristimo Gibbs-ovo uzorkovanje, možemo npr. simulirati iz marginalne raspodele za \(X\), a onda iz uslovne raspodele za \(Y|X\). To možemo uraditi na sledeći način:

binorm = function (n, rho) {
        x = rnorm(n, 0, 1)
        y = rnorm(n, rho * x, sqrt(1 - rho^2))
        cbind(x, y)
}

Ovo kreira vektor vrednosti za \(X\), a onda njih koristi za kreiranje vrednosti vektora \(Y\), pod uslovom dobijenih vrednosti za \(X\). Smeštamo rezultat u \(n \times2\) matricu:

uzorak = binorm(10000,0.98)

par(mfrow=c(2,1))
plot(uzorak,col=1:10000)
plot(uzorak,type="l")

plot(ts(uzorak[,1]))
plot(ts(uzorak[,2]))

# histogrami pojedinacnih komponenti
hist(uzorak[,1],40)
hist(uzorak[,2],40)

Hajde sada da vidimo kakve bismo rezultate dobili ako bismo koristili Gibbs-ovo uzorkovanje. Pravimo Gibbs-ov uzorak uzastopnim uzorkovanjem iz uslovnih raspodela.

gibbs = function(n, rho) {
  m = matrix(nrow = n, ncol = 2)
  x = 0
  y = 0
  m[1, ] = c(x, y)
  for(i in 2:n) {
    x = rnorm(1, rho*y, sqrt(1-rho^2))
    y = rnorm(1, rho*x, sqrt(1-rho^2))
    m[i, ] = c(x, y)
  }
  m
}

Dakle, prvo smo kreirali matricu u kojoj čuvamo rezultate, a zatim inicijalizovali lanac na \((0,0)\). Petlja uzostopno uzorkuje iz uslovnih raspodela, čuvajući rezulatate u matrici. Prikažimo rezultate:

uzorak = gibbs(10000, 0.98)

par(mfrow=c(2,1))
plot(uzorak,col=1:10000)
plot(uzorak,type="l")

plot(ts(uzorak[,1]))
plot(ts(uzorak[,2]))

hist(uzorak[,1],40)
hist(uzorak[,2],40)

Uz malo sreće, ovo daje rezultate koji izgledaju veoma slično onima koji su dobijeni ranije, osim grafika vremenskih serija marginalnih raspodela, koji pokazuju različitu autokorelaciju između uzastopnih vrednosti.

Primer 2

Pretpostavimo da: \[Y_1,...,Y_n \overset{iid}{\sim} \mathcal{N}(\mu, \tau^{-1})\] \[\mu \sim \mathcal{N}(0,1)\] \[\tau \sim \mathcal{\gamma}(2,1)\] Pretpostavimo još i da su \(\mu\) i \(\tau\) nezavisne (tj. da je apriorna gustina za \((\mu,\tau)\) jednaka proizvodu pojedinačnih gustina).

Hajde sada da pronađemo uslovne raspodele za za \(\mu\) i \(\tau\).

Za početak, izračunajmo funkciju verodostojnosti (likelihood) (što je ukupna gustina podataka (za date parametre)), kao funkciju od parametara: \[L(\mu,\tau) = \prod_{i=1}^n f(Y_i|\mu,\tau)=\prod_{i=1}^n \frac{\sqrt{\tau}}{\sqrt{2\pi}}e^{-\frac{\tau(Y_i-\mu)^2}{2}} = (2\pi)^{-\frac{n}{2}}\tau^{\frac{n}{2}}e^{-\tau \frac{ \sum_{i=1}^n(Y_i-\mu)^2}{2}}\] Da bismo pronašli zajedničku gustinu \((\textbf{Y},\mu,\tau)\), pomnožimo \(f(\textbf{Y}|\mu,\tau)\) (likelihood) aproirnom gustinom: \[p(\textbf{Y},\mu,\tau) = (constant) \times \tau^{\frac{n}{2}}e^{-\tau \frac{ \sum_{i=1}^n(Y_i-\mu)^2}{2}} e^{-\frac{\mu^2}{2}}\tau e^{-\tau}\] Za dobijene podatke \(\textbf{Y}\), u Bajesovskom pristupu, naša pažnja se fokusira na uslovnu raspodelu od \((\mu,\tau)\) za dato \(\textbf{Y}\). Nju možemo računati na sledeći način: \[p(\mu,\tau|\textbf{Y}) = \frac{p(\textbf{Y},\mu,\tau)}{p(\textbf{Y})} = \frac{p(\textbf{Y},\mu,\tau)}{\int \int p(\textbf{Y},\mu^*,\tau^*)d\mu^*d\tau^*}\]

Međutim, imenilac nam nije od nekog značaja trenutno (jer ne uključuje \(\mu\) ili \(\tau\)), tako da možemo jednostavno pisati: \[p(\mu,\tau|\textbf{Y}) \propto p(\textbf{Y},\mu,\tau)\] Sada konačno možemo govoriti o uslovnim raspodelama. Za \(\mu\) dobijamo: \[p(\mu|\tau,\textbf{Y}) = \frac{p(\mu,\tau|\textbf{Y})}{p(\tau|\textbf{Y})}\] Međutim, imenilac ponovo ne uključuje \(\mu\) pa možemo da zapišemo jednostavnije, kao funkciju od \(\mu\): \[p(\mu|\tau,\textbf{Y}) \propto p(\textbf{Y},\mu,\tau) \propto e^{-\tau \frac{ \sum_{i=1}^n(Y_i-\mu)^2}{2}} e^{-\frac{\mu^2}{2}}\]

(Zanemarili smo sve faktore u \(p(\textbf{Y},\mu,\tau)\) koji ne uključuju \(\mu\).)

Slično, kao funkciju od \(\tau\),uslovnu gustinu za \(\tau\) možemo pisati na sledeći način: \[p(\tau|\mu,\textbf{Y}) \propto p(\textbf{Y},\mu,\tau) \propto \tau^{\frac{n}{2}} e^{-\tau \frac{ \sum_{i=1}^n(Y_i-\mu)^2}{2}} \tau e^{-\tau}\]

Ako ovo malo sredimo, dobijamo sledeći izraz za uslovnu gustinu za \(\mu\): \[p(\mu|\tau,\textbf{Y}) = e^{-\frac{1+n\tau}{2} (\mu - \frac{\tau\sum_{i=1}^nY_i}{1+n\tau})^2} \times \{\text{nešto što ne zavisi od }\mu\}\]

Ovo prepoznajemo! Ovo je gustina normalne raspodele u tački \(\mu\), sa očekivanjem \(\tau\frac{\sum_{i=1}^nY_i}{1+n\tau}\) i disperzijom \(\frac{1}{1+n\tau}\). Dakle: \[\mu|\tau,\textbf{Y} \sim \mathcal{N}(\tau\frac{\sum_{i=1}^nY_i}{1+n\tau}, \frac{1}{1+n\tau})\]

Slično, nakon malo izračunavanja dobijamo izraz za uslovnu gustinu za \(\tau\):

\[p(\tau|\mu,\textbf{Y}) = \tau^{1+\frac{n}{2}} e^{-\tau(1+\frac{1}{2}\sum_{i=1}^n(Y_i-\mu)^2) e^{-\frac{1}{\tau}}} \times \{\text{nešto što ne zavisi od }\tau\}\] I ovo nam je poznato - u pitanju je gama raspodela sa odgovarajućim parametrima. Dakle: \[\tau|\mu,\textbf{Y} \sim \mathcal{\gamma}(2+\frac{n}{2},1+\frac{1}{2}\sum_{i=1}^n(Y_i-\mu)^2)\] Iskoristimo sada ove činjenice da napravimo Markovski lanac. Naš cilj je da napravimo lanac takav da je njegova stacionarna raspodela aposteriorna raspodela \(p(\mu,\tau|\textbf{Y})\) koja nas interesuje.

# podaci Y
y = scan("http://personal.psu.edu/drh20/515/hw/MCMCexampleData.txt")

Prvo, treba da postavimo neke početne vrednosti \(\mu^{(1)}\) i \(\tau^{(1)}\). U idealnom slučaju, uzeli bismo uzorak iz prave stacionarne (aposteriorne) raspodele, ali ako bismo to mogli da uradimo, ne bi nam ni bio potreban MCMC!

Uzmimo jednostavno \((\mu^{(1)}, \tau^{(1)}) = (1,2)\), što su sredine apriornih raspodela. (Mogli smo, takođe, da pogledamo podatke i kao početne vrednosti uzmemo uzoračku sredinu i recipročnu vrednost uzoračke disperzije, na primer.)

i = 1 # iteracije
mu = 1
tau = 2
sumOfY = sum(y) # Ova vrednost ce nam trebati vise puta
n = length(y)

Sada ćemo pokrenuti 1000 iteracija da dobijemo lanac Markova, tako što uzorkujemo \(\mu^{(i+1)}\) iz uslovne raspodele za \(\mu\) u odnosu na \(\tau^{(i)}\) i \(\textbf{Y}\), a onda uzorkujemo \(\tau^{(i+1)}\) iz uslovne raspodele za \(\tau\), u odnosu na \(\mu^{(i+1)}\) i \(\textbf{Y}\):

while (i <= 1000) {
  tmp = tau[i]
  mu = c(mu, rnorm(1, mean = sumOfY*tmp/(1 + n*tmp), sd = sqrt(1/(1 + n*tmp))))
  tau = c(tau, rgamma(1, shape = 2 + n/2, rate = 1+sum((y - mu[i+1])^2)/2))
  i = i + 1
}

Proverimo sada neka od svojstava našeg MCMC uzorka. Da bismo videli kako naš uzorak izgleda, nacrtaćemo jednostavni grafik:

plot(cbind(mu, tau))

Prva stvar koju vidimo je da je početna tačka daleko od ostalih tačaka. To je delom razlog zbog koga se često koristi burn-in period. U našem slučaju, izostavljanje prvih nekoliko tačaka čini se dovoljnim.

Pogledajmo sada još neke grafike da bismo procenili da li se Markovski lanac “dobro meša”. Jednostavni grafici iteracija u odnosu na \(\mu\) ili u odnosu na \(\tau\) (bez prvih 5 tačaka) daju prilično dobar uvid u ovo.

par(mfrow=c(2,1))
plot(mu[-(1:5)], type="l")
plot(tau[-(1:5)], type="l")

Grafici ne prikazuju nikakav opšti trend, koji bi mogao ukazivati na to da lanac ipak nije imao priliku da adekvatno istraži aposteriorni parametarski prostor. (Napomena: Grafik ne daje nikakav garanciju! Neki vrste neobičnog ponašanja se možda neće pojaviti na grafiku.)

Korisno je da pogledamo i korelogram. U ovom slučaju, autokorelacija je bliska nuli, čak i za korak 1, što znači da su uzastopne uzorkovane vrednosti ovog lanca skoro nekorelisane. Ovo je vrlo dobra stvar!

par(mfrow=c(1,2))
acf(mu)
acf(tau)