3 Markovin ketju Monte Carlo

Bayes-tilastotieteessä Markovin ketju Monte Carloa (Markov chain Monte Carlo, MCMC) käytetään havaintojen simulointiin posteriorijakaumasta. Tämä tapahtuu rakentamalla Markovin ketju, jonka tasapainojakauma on kiinnostuksen kohteena oleva posteriorijakauma.

3.1 Markovin ketju

Määritelmä 3.1 (Markovin ketju) Satunnaismuuttujien jonoa \(\theta^{(0)},\theta^{(1)},\theta^{(2)},\ldots\,\) kutsutaan Markovin ketjuksi mikäli kaikilla \(t > 0\) \[ p(\theta^{(t+1)} | \theta^{(0)},\ldots,\theta^{(t)}) = p(\theta^{(t+1)} | \theta^{(t)}). \]

Satunnaismuuttujien \(\theta^{(0)},\theta^{(1)},\theta^{(2)},\ldots\) mahdollisia arvoja kutsutaan usein tiloiksi ja näiden muodostamaa joukkoa tila-avaruudeksi. Markovin ketjussa tilasiirtymien todennäköisyydet (tai tiheydet, jos tila-avaruus on jatkuva) riippuvat ainoastaan nykyisestä tilasta. Sillä, kuinka nykyiseen tilaan on päädytty, ei ole merkitystä tuleviin tilasiirtymiin. Tiettyjen ehtojen vallitessa voidaan osoittaa, että kun \(t \rightarrow \infty\), satunnaismuuttujan \(\theta^{(t+1)}\) jakauma suppenee kohti tasapainojakaumaa, joka ei riipu ketjun alkuarvosta \(\theta^{(0)}\).

Tehtävänä on siis luoda Markovin ketju, jonka tasapainojakauma on haluttu posteriorijakauma. Tämän toteuttamiseen on olemassa useita algoritmeja, joista tärkeimpiä käsitellään jatkossa. Lisäksi tulee valita Markovin ketjun alkuarvo ja selvittää milloin ketju on supennut tasapainojakaumaansa.

3.2 Globaali ja lokaali tasapainoehto

Tarkastellaan Markovin ketjua, jonka tila-avaruus \(S\) on diskreetti. Tällöin siirtymää tilasta \(i\) tilaan \(j\) voidaan merkitä todennäköisyydellä \(p_{ij}\).

Määritelmä 3.2 (Tasapainojakauma) Markovin ketjun tasapainojakauma on jakauma \(\{\pi_i, i \in S\}\), joka toteuttaa globaalin tasapainoehdon \[ \sum_{i \in S} \pi_i p_{ij} = \pi_j \textrm{ kaikille } j\in S. \]

Globaalissa tasapainoehdossa tarkastellaan siis erilaisia tapoja saapua tilaan \(j\). Kun kunkin lähtötilan \(i\) todennäköisyys määrätään tasapainojakaumasta, päädytään tilaan \(j\) tasapainojakauman mukaisella todennäköisyydellä \(\pi_j\).

Tasapainojakauman olemassaolo voidaan osoittaa, mikäli tietyt ketjua koskevat ehdot täyttyvät. Näistä ehdoista tärkein on pelkistymättömyys (irreducibility), joka takaa sen, että siirtyminen mistä tahansa tilasta \(i\) mihin tahansa tilaan \(j\) on mahdollista äärellisellä määrällä siirtymiä.

Globaalin tasapainoehto ei ole kovinkaan käyttökelpoinen, kun tavoitteena on konstruoida Markovin ketju, jolla on annettu tasapainojakauma. Simulointialgoritmit perustuvatkin vahvempaan lokaaliin tasapainoehtoon \[ \pi_i p_{ij} = \pi_j p_{ji} \textrm{ kaikilla } i,j \in S. \] Lokaalissa tasapainoehdossa tarkastellaan siis siirtymää kahden tilan välillä. Ehdon mukaan tasapainotilanteessa on yhtä todennäköistä havaita siirtymä tilasta \(i\) tilaan \(j\) kuin tilasta \(j\) tilaan \(i\). Lokaalista tasapainosta seuraa globaali tasapaino seuraavalla tavalla \[ \begin{aligned} \sum_{i \in S} \pi_i p_{ij} = \sum_{i \in S} \pi_j p_{ji} = \pi_j \sum_{i \in S} p_{ji} = \pi_j. \end{aligned} \]

3.3 Gibbsin algoritmi

Gibbsin algoritmissa Markovin ketju rakennetaan siten, että mallin parametreja päivitetään yksi kerrallaan. Uusi arvo generoidaan ehdollisesta jakaumasta, jossa ehtoina ovat muut parametrit. Aluksi mallin parametrit jaetaan \(k\) komponenttiin \(\theta=(\theta_1,\theta_2,\ldots,\theta_k)\). Yleensä komponentit ovat skalaariparametreja, mutta ne voivat olla myös vektoreita. Gibbsin algoritmin vaiheet ovat:

Aseta alkuarvot \(\theta_1^{(0)},\theta_2^{(0)},\ldots,\theta_k^{(0)}\) ja aseta \(t=0\)
Generoi ensin \(\theta_1^{(t+1)}\) jakaumasta \(p(\theta_1 | \theta_2^{(t)},\ldots,\theta_k^{(t)})\) ja sitten \(\theta_2^{(t+1)}\) jakaumasta \(p(\theta_2 | \theta_1^{(t+1)},\theta_3^{(t)},\ldots,\theta_k^{(t)})\) (jossa \(\theta_1\) on jo siis päivitetty). Jatka näin kunnes kaikki komponentit on päivitetty ja kasvata sitten iteraatiota \(t\) yhdellä.
Toista vaihetta 2 riittävän monta kertaa (yleensä tarvitaan tuhansia iteraatioita).

Esimerkki 3.1 (Kaksiulotteisen normaalijakauman simulointi) Oletetaan, että \[ \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} \sim \operatorname{N}_2(\mu, \Sigma),\, \mu = \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix},\, \Sigma = \begin{pmatrix} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_2^2 \end{pmatrix}. \] On laskettava ehdolliset jakaumat. Tässä tapauksessa saadaan: \[ \begin{aligned} X_1|X_2 &= x_2 \sim \operatorname{N}\left(\mu_1 + \frac{\sigma_1}{\sigma_2}\rho(x_2-\mu_2),(1-\rho^2)\sigma_1^2 \right) \\ X_2|X_1 &= x_1 \sim \operatorname{N}\left(\mu_2 + \frac{\sigma_2}{\sigma_1}\rho(x_1-\mu_1),(1-\rho^2)\sigma_2^2 \right), \end{aligned} \] missä \(\rho = \dfrac{\sigma_{12}}{\sigma_1 \sigma_2}\) on korrelaatiokerroin.

Simulointi etenee seuraavasti: Olkoon alkuarvo \((x_1^{(0)},x_2^{(0)})^T\). \[ \begin{pmatrix} x_1^{(0)} \\ x_2^{(0)} \end{pmatrix} \rightarrow \begin{pmatrix} x_1^{(1)} \\ x_2^{(0)} \end{pmatrix} \rightarrow \begin{pmatrix} x_1^{(1)} \\ x_2^{(1)} \end{pmatrix} \rightarrow \begin{pmatrix} x_1^{(2)} \\ x_2^{(1)} \end{pmatrix} \rightarrow \cdots \] Simuloinnin voi toteuttaa esimerkiksi R-ohjelmistolla seuraavasti

norm_sim <- function(mu, S, x_0, n_sim, n_burn) {
  x <- matrix(ncol = 2, nrow = n_burn + n_sim)
  x[1, ] <- x_0          # alkuarvo
  s1 <- sqrt(S[1, 1])    # x1:n keskihajonta
  s2 <- sqrt(S[2, 2])    # x2:n keskihajonta
  s12 <- S[1, 2]         # kovarianssi
  rho <- s12 / (s1 * s2) # korrelaatiokerroin
  for (i in seq(2, n_burn + n_sim)) {
    x[i, 1] <- rnorm(
      n = 1, 
      mean = mu[1] + s1 / s2 * rho * (x[i - 1, 2] - mu[2]), 
      sd = sqrt(1 - rho^2) * s1
    )
    x[i, 2] <- rnorm(
      n = 1, 
      mean = mu[2] + s2 / s1 * rho * (x[i, 1] - mu[1]), 
      sd = sqrt(1 - rho^2) * s2
    )
  }
  x[seq(n_burn + 1, n_burn + n_sim), ]
}

mu <- c(4, 1)
S <- matrix(c(25, 10.5, 10.5, 9), ncol = 2)
out <- norm_sim(mu, S, c(0, 0), 5e3, 1e3)

Jakaumaa voidaan kuvailla esimerkiksi sen tunnuksilla, marginaalijakaumien histogrammeilla sekä otoksen hajontakuviolla.

	Mean	SD	2.5 %	97.5 %
\(x_1\)	3.7759299	4.954002	-6.150616	13.694780
\(x_2\)	0.8515133	2.959983	-4.922519	6.578646

3.4 Metropolis-algoritmi

Metropolis-algoritmissa käytetään ehdotusjakaumaa, josta generoidaan ehdotuksia Markovin ketjun seuraavaksi tilaksi. Ehdotus voidaan hyväksyä tai hylätä. Jälkimmäisessä tapauksessa seuraava tila on sama kuin nykyinen tila. Algoritmi perustuu siihen, että posteriorijakaumien suhde voidaan laskea normeeraamattomien posteriorijakaumien avulla, vaikka Bayesin kaavan nimittäjää \(p(y)\) ei tunnettaisikaan. Olkoon \(\theta_a\) ja \(\theta_b\) kaksi parametrin \(\theta\) arvoa. Posteriorijakaumien suhteelle pätee \[ \frac{ p(\theta_b | y) }{ p(\theta_a | y)} = \frac{ p(\theta_b)p(y | \theta_b)/p(y) }{ p(\theta_a)p(y | \theta_a)/p(y)} = \frac{ p(\theta_b)p(y | \theta_b)}{p(\theta_a)p(y | \theta_a)}. \] Algoritmi voidaan esittää seuraavasti:

Aseta alkuarvo \(\theta^{(0)}\) ja aseta \(t=0\).
Toista vaiheita 3–5, iteraatioilla \(t = 1, 2, \ldots, T\).
Generoi ehdotus \(\theta^*\) ehdotusjakaumasta \(q(\theta^* | \theta^{(t)})\). Ehdotusjakauman tulee täyttää symmetrisyysehto \(q(\theta_a | \theta_b) = q(\theta_b | \theta_a)\) kaikilla \(\theta_a\), \(\theta_b\).
Laske hyväksymistodennäköisyys normeeraamattomien posteriorijakaumien avulla \[ \alpha = \min \left\{1, \frac{p(\theta^*)p(y | \theta^*)}{p(\theta^{(t)})p(y | \theta^{(t)})} \right\} \]
Aseta \[ \theta^{(t+1)} = \begin{cases} \theta^* & \text{todennäköisyydellä } \alpha, \\ \theta^{(t)} & \text{muutoin.} \end{cases} \]

Metropolis-algoritmi toteuttaa lokaalin tasapainoehdon. Tämän osoittamiseksi tarkastellaan Markovin ketjun siirtymätodennäköisyyksiä
\[ P(\theta^* | \theta^{(t)}) = \begin{cases} q(\theta^* | \theta^{(t)}) \min \left\{1, \frac{p(\theta^*)p(y | \theta^*)}{p(\theta^{(t)})p(y | \theta^{(t)})} \right\}, & \theta^* \neq \theta^{(t)}, \\ q(\theta^{(t)} | \theta^{(t)}) + \sum_{\theta' \neq \theta^{(t)}} q(\theta' | \theta^{(t)}) \left(1- \min \left\{1, \frac{p(\theta')p(y | \theta')}{p(\theta^{(t)})p(y | \theta^{(t)})} \right\}\right), & \theta^* = \theta^{(t)}. \end{cases} \] Lokaali tasapaino edellyttää ehdon \(p(\theta^{(t)} | y)P(\theta^* | \theta^{(t)})=p(\theta^* | y)P(\theta^{(t)} | \theta^*)\) voimassaoloa. Jos \(\theta^* = \theta^{(t)}\), ehto toteutuu suoraan. Tilanteessa \(\theta^* \neq \theta^{(t)}\) saadaan \[ \begin{aligned} p(\theta^{(t)} | y)P(\theta^* | \theta^{(t)}) & = p(\theta^{(t)} | y) q(\theta^* | \theta^{(t)}) \min \left\{1, \frac{p(\theta^* | y)}{p(\theta^{(t)} | y)} \right\} \\ & = q(\theta^{(t)} | \theta^*) \min \left\{p(\theta^* | y),p(\theta^{(t)} | y) \right\} \\ & = p(\theta^* | y) q(\theta^{(t)} | \theta^*) \min \left\{1, \frac{p(\theta^{(t)} | y)}{p(\theta^* | y)} \right\} \\ & = p(\theta^* | y)P(\theta^{(t)} | \theta^*), \end{aligned} \] missä toisella rivillä on käytetty ehdotusjakauman symmetrisyysominaisuutta \(q(\theta^* | \theta^{(t)})=q(\theta^{(t)} | \theta^*)\).

3.5 Metropolis-Hastings-algoritmi

Metropolis-Hastings-algoritmissa Metropolis-algoritmia muokataan siten, että ehdotusjakauman ei tarvitse täyttää symmetrisyysehtoa. Algoritmi voidaan esittää seuraavasti:

Aseta alkuarvo \(\theta^{(0)}\) ja aseta \(t=0\).
Toista vaiheita 3–5, iteraatioilla \(t=1,2,\ldots,T\).
Generoi ehdotus \(\theta^*\) ehdotusjakaumasta \(q(\theta^* | \theta^{(t)})\). Ehdotusjakauman ei tarvitse täyttää symmetrisyysehtoa.
Laske hyväksymistodennäköisyys ehdotusjakauman ja normeeraamattomien posteriorijakaumien avulla \[ \alpha = \min \left\{1, \frac{p(\theta^*)p(y | \theta^*)q(\theta^{(t)} | \theta^*)}{p(\theta^{(t)})p(y | \theta^{(t)})q(\theta^* | \theta^{(t)})} \right\} \]
Aseta \[ \theta^{(t+1)}= \begin{cases} \theta^* & \text{todennäköisyydellä } \alpha, \\ \theta^{(t)} & \text{muutoin.} \end{cases} \]

Metropolis-Hastings-algoritmia voidaan käyttää Gibbsin algoritmin osana. Parametri, jonka ehdollista jakaumaa ei kyetä laskemaan, voidaan päivittää Metropolis-Hastings-algoritmilla. Toisaalta Gibbsin algoritmi voidaan nähdä Metropolis-Hastings-algoritmin varianttina, jossa ehdotukset tuotetaan komponenteittain ehdollisista jakaumista ja hyväksytään todennäköisyydellä 1.

3.6 Ketjun suppenemisen tarkastelu

MCMC-menetelmiä käytettäessä tulee varmistaa, että ketju on supennut tasapainojakaumaan. Ketjun alkupään havainnot eivät ole peräisin tasapainojakaumasta eikä niitä tästä syystä tule käyttää posteriorijakaumaa estimoitaessa. Ketjun alkua kutsutaan nimellä lämmitys tai sisäänajo (englanniksi burn-in tai warm-up). Monet ohjelmistot poistavat oletusarvoisesti ketjun ensimmäisen puoliskon, ellei käyttäjä toisin määrää.

Ensimmäinen tapa arvioida ketjun suppenemista on silmämääräinen tarkastelu. Ketjun tulisi näyttää satunnaiselta vaihtelulta vakiona pysyvän keskiarvon ympärillä. Silmämääräisen tarkastelun avulla on usein mahdollista todeta, että ketju ei ole vielä supennut, mutta suppenemisen varmistaminen pelkästään tällä tavalla ei ole täysin luotettavaa. On yleistä, että jokin mallin parametri on supennut mutta jokin toinen parametri ei ole. Kaikkien parametrien systemaattinen tarkastelu kuvaajien avulla voi viedä paljon aikaa, koska Bayes-malleissa usein satoja tai tuhansia parametreja.

Edistyneemmät suppenemistarkastelut vaativat usean eri ketjun käyttämistä. Ketjujen alkuarvojen tulisi selvästi poiketa toisistaan, jotta alkuarvojen vaikutusta suppenemiseen voidaan tarkastella. Samaa alkuarvoa käytettäessä osa mahdollisista ongelmista voi jäädä huomaamatta. Ketjujen käyttäytymistä voi edelleen tarkastella silmämääräisesti esimerkiksi piirtämällä (saman parametrin) eri ketjut samaan kuvaan. Suppenemisen tarkastelemiseen on myös ehdotettu erilaisia tunnuslukuja, joista seuraavaksi esitellään usean ohjelmiston käyttämä tunnusluku mahdollisen skaalan pienenemisen kerroin \(\hat{R}\) (potential scale reduction factor).

Tunnusluku perustuu ketjujen sisäisen ja ketjujen välisen vaihtelun vertaamiseen ja siitä on esitetty useita eri versioita kirjallisuudessa. Seuraava määritelmä on niin sanottu jaettujen ketjujen \(\hat{R}\) (split \(\hat{R}\)). Idea tässä tunnusluvun variantissa on, että sen tulisi huomioida samanaikaisesti ketjujen mahdollinen epästationaarisuus ja huono sekoittuminen. Olkoon käytössä \(m_0\) ketjua, joista jokaisesta on simuloitu \(n_0\) havaintoa. Jätetään nyt ensimmäinen puolikas havainnoista pois lämmittelyjaksona (eli \(n_0 / 2\) havaintoa). Jaetaan nyt jokainen ketju kahteen osaan siten, että ensimmäiseen osaan valitaan ensimmäiset \(n_0/4\) jäljelle jääneistä havainnoista, ja toiseen osaan viimeiset \(n_0/4\) havaintoa. Ajatellaan nyt jaettuja ketjuja erillisinä siten, että käsiteltävä ketjujen määrä on \(m = 2m_0\) joista jokaisesta on \(n = n_0/4\) havaintoa.

Skalaariparametrista \(\theta\) on käytössä simuloidut arvot \(\theta^{(ij)}\), \(i=1,\ldots,n\), \(j=1,\ldots,m\). Ketjujen välistä vaihtelua kuvaa neliösumma \[ B = \frac{n}{m-1} \sum_{j=1}^m (\overline{\theta}^{(\cdot j)}-\overline{\theta}^{(\cdot \cdot)})^2, \textrm{ missä } \overline \theta^{(\cdot j)}=\frac{1}{n} \sum_{i=1}^n \theta^{(ij)} \textrm{ ja } \overline \theta^{(\cdot \cdot)}=\frac{1}{m} \sum_{j=1}^m \overline \theta^{(\cdot j)}. \tag{3.1}\] Ketjujen sisäistä vaihtelua kuvaa neliösumma \[ W=\frac{1}{m} \sum_{j=1}^m s_j^2, \textrm{ missä } s_j^2=\frac{1}{n-1} \sum_{i=1}^n (\theta^{(ij)}-\overline{\theta}^{(\cdot j)})^2. \] Sisäinen varianssi \(W\) aliestimoi posteriorivarianssia \(\operatorname{Var}(\theta | y)\), koska äärellinen ketju ei välttämättä ole käynyt läpi koko parametriavaruutta. Sen sijaan estimaattori \[ \widehat{\operatorname{Var}}^+(\theta|y) = \frac{n-1}{n}W + \frac{1}{n}B \] yliestimoi posteriorivarianssia. Sekä \(W\) että \(V\) ovat posteriorivarianssin tarkentuvia estimaattoreita. Tunnusluku \(\hat{R}\) perustuu näiden kahden estimaattorin suhteeseen \[ \hat{R} = \sqrt{\frac{\widehat{\operatorname{Var}}^+(\theta|y)}{W}}, \] joka suppenee kohti arvoa 1, kun \(n \rightarrow \infty\). Käytännössä ketjun katsotaan supenneen, jos \(\hat{R}<1.05\) kaikille parametreille.

3.7 Tehokas otoskoko

Koska Markovin ketjut tuottavat riippuvia havaintoja, otos sisältää tyypillisesti vähemmän informaatiota kuin samankokoinen otos riippumattomia havaintoja. Tehokkaan otoskoon tarkoituksena on kertoa, kuinka suurta riippumatonta otosta käytettävissä oleva ketju vastaa. Esimerkiksi Stan estimoi tehokkaan otoskoon automaattisesti. Tehokkaan otoskoon voi estimoida esimerkiksi myös coda-paketin funktiolla effectiveSize.

Mikäli ketju on tasapainotilassa, riippuvuudella ei ole vaikutusta odotusarvoon. Ketjujen välinen vaihtelu \(B\) (kaava 3.1) sen sijaan on suurempaa kuin riippumattomalla otoksella. Mitä suurempi ketjujen autokorrelaatio on, sitä enemmän \(B\) poikkeaa riippumattoman otoksen vaihtelusta.

Tehokas otoskoko määritellään ketjun autokorrelaation avulla eri viiveitä (lags) tarkastelemalla. Viiveen \(k \geq 0\) autokorrelaatio \(\rho_k\) ketjulle, jonka yhteisjakauman on \(p(\theta)\) odotusarvolla \(\mu\) ja varianssilla \(\sigma^2\), määritellään seuraavasti \[ \rho_k = \frac1{\sigma^2} \int_S (\theta^{(t)} - \mu)(\theta^{(t + k)} - \mu)p(\theta)\, d\theta. \tag{3.2}\] Autokorrelaatio mittaa ketjujen \(\theta^{(t)}\) ja \(\theta^{(t+k)}\) välistä korrelaatiota, eli siis alkuperäisen ketjun \(\theta^{(t)}\) ja siitä \(k\) indeksiä siirtämällä saadun ketjun välistä korrelaatiota.

Otoskokoa \(n\) vastaava tehokas otoskoko määritellään autokorrelaatioiden avulla seuraavasti: \[ n_\textrm{eff} = \frac{n}{\sum_{k = -\infty}^\infty \rho_k} = \frac{n}{1 + 2\sum_{k = 1}^\infty \rho_k}. \] Otoksesta laskettuja autokorrelaatioita käyttämällä tehokkaalle otoskoolle voidaan johtaa estimaattori \[ \hat{n}_\textrm{eff}=\frac{mn}{1+2 \sum_{k=1}^K \hat{\rho}_k}, \] missä \(\hat{\rho}_k\) on estimaatti viiveen \(k\) autokorrelaatiolle joka on laskettu \(m\) ketjusta. Estimoinnissa käytetty maksimiviive \(K\) on pienin pariton luku, jolle \(\hat{\rho}_{K+1}+\hat{\rho}_{K+2}<0\).