MCMC-simuloinnin lisäksi posteriorijakauman estimointiin on olemassa muitakin menetelmiä. Tiettyihin tilanteisiin on olemassa menetelmiä, jotka saattavat olla huomattavasti MCMC-simulointia nopeampia. Toisaalta approksimointi on ainoa vaihtoehto silloin, kun malli niin monimutkainen ettei sitä laskennallisten rajoitteiden takia ole käytännössä mahdollista sovittaa MCMC-menetelmin.
Kun mallin parametrien lukumäärä on pieni, käytännössä yksi tai kaksi, integraalin \(\int p(\theta)p(y | \theta)\, d\theta\) voi laskea numeerisella integroinnilla. Tähän soveltuu esimerkiksi R-ohjelman integrate-funktio. Parametrien lukumäärän kasvaessa numeerisen integroinnin laskennallinen vaativuus kasvaa nopeasti. Moniulotteisessa tapauksessa voidaan käyttää esimerkiksi R pakettia cubature, joka sisältää useita adaptiivisia menetelmiä numeeriseen integrointiin.
Esimerkki 7.1 (Genetiikkaesimerkki) Koe-eläimet \((n = 197)\) on jaettu neljään luokkaan geneettisin perustein:
| Kategoria | AA | AB | BA | BB |
|---|---|---|---|---|
| Frekvenssi | 125 | 18 | 20 | 34 |
Olkoon \(y_i\) luokan \(i\) frekvenssi. Oletetaan, että luokittelu tapahtuu riippumattomasti, jolloin yleisin malli havaintovektorille \(y = (y_1,y_2,y_3,y_4)\) on multinomiaalimalli: \[ y_1,y_2,y_3,y_4|\pi_1,\pi_2,\pi_3,\pi_4 \sim \operatorname{Mult}(n; \pi_1,\pi_2,\pi_3,\pi_4), \] missä \(\sum_{i=1}^4 y_i\) = n, \(0 < \pi_i < 1,\, i = 1,2,3,4\) ja \(\sum_{i=1}^4 \pi_i = 1\).
Geneettisen teorian perusteella asetetaan malli todennäköisyyksille \(\pi_i\) \[ \pi = \left(\frac12 + \frac{\theta}4, \frac14(1-\theta),\frac14(1-\theta),\frac{\theta}4 \right). \] Mallissa on siten yksiulotteinen parametri \(\theta\). Huomataan, että summaehto täyttyy aina ja lisäksi voidaan todeta, että \(0 < \theta < 1\).
Valitaan prioriksi tasajakauma \(p(\theta) \propto 1\). Posterioriksi saadaan (ilman normeerausta) \[ p(\theta|y) \propto \left(\frac12 + \frac{\theta}4 \right)^{y_1}\left(\frac14(1-\theta)\right)^{y_2} \left(\frac14(1-\theta)\right)^{y_3} \left(\frac{\theta}4\right)^{y_4}. \] Lasketaan \[ p(y) = \int_0^1 \left(\frac12 + \frac{\theta}4 \right)^{y_1}\left(\frac14(1-\theta)\right)^{y_2} \left(\frac14(1-\theta)\right)^{y_3} \left(\frac{\theta}4\right)^{y_4}\, d\theta \] numeerisella integroinnilla.
y <- c(125, 18, 20, 34)
# Normeeraamaton posteriori
p_un <- function(theta, y) {
(0.5 + 0.25 * theta)^y[1] * (0.25 * (1 - theta))^y[2] *
(0.25 * (1 - theta))^y[3] * (0.25 * theta)^y[4]
}
I <- integrate(p_un, lower = 0, upper = 1, y = y)$value
I[1] 5.84345e-91Parametrin \(\theta\) posteriorijakauman kuvaaja on
Numeerisella integroinnilla saadaan posteriorijakaumalle posteriorikeskiarvo \(\approx 0.623\) ja posteriorikeskihajonta \(\approx 0.051\).
Tiettyjen säännöllisyysehtojen vallitessa voidaan osoittaa, että posteriorijakauma lähestyy asymptoottisesti normaalijakaumaa, kun havaintojen määrä kasvaa. Tätä ominaisuutta hyödynnetään Laplace-approksimaatiossa. Jos havaintojen määrä on pieni, approksimaatio ei välttämättä toimi luotettavasti.
Laplace-approksimaatio perustuu posteriorimaksimin \(\hat{\theta}\) ympärille kehitettyyn Taylorin sarjaan. Merkitään \(\theta = (\theta_1,\ldots,\theta_m)^\top\). Taylorin kehitelmässä käytetään normeeraamattoman posteriorijakauman logaritmia \[ h(\theta) = \log [p(y | \theta) p(\theta)]. \] Lasketaan ensimmäisen ja toisen kertaluvun osittaisderivaatat \[ \begin{aligned} h'(\theta) =& \frac{\partial h(\theta)}{\partial \theta} = \left(\frac{\partial h(\theta)}{\partial \theta_1}, \ldots, \frac{\partial h(\theta)}{\partial \theta_m} \right)^T \\ h''(\theta) =& \frac{\partial^2 h(\theta)}{\partial \theta^2} = \left( \frac{\partial^2 h(\theta)}{\partial \theta_r \partial \theta_s} \right), \quad r,s=1,\ldots,m \end{aligned} \] missä \(h''(\theta)\) on siis \(m \times m\)-matriisi. Posteriorimaksimille \(\hat{\theta}\) pätee \(h'(\hat{\theta})=0\). Taylorin kehitelmästä saadaan \[ h(\theta) \approx h(\hat{\theta})+(\theta-\hat{\theta})^\top h'(\hat\theta)+ \frac{1}{2} (\theta-\hat{\theta})^\top h''(\hat\theta) (\theta-\hat{\theta}). \] Kun \(\theta-\hat{\theta}\) on pieni, korkeamman kertaluvun termit voidaan jättää ottamatta huomioon. Tällöin posteriorijakaumalle pätee \[ p(\theta | y) \propto \exp(h(\theta)) \propto \exp \left(\frac{1}{2}(\theta-\hat{\theta})^\top h''(\hat\theta) (\theta-\hat{\theta}) \right). \] Kirjoittamalla \(V=(-h''(\hat\theta))^{-1}\) saadaan muoto \[ p(\theta | y) \propto \exp \left( -\frac{1}{2}(\theta-\hat{\theta})^\top V^{-1} (\theta-\hat{\theta}) \right), \] joka vastaa normeeraamatonta normaalijakaumaa. Laplace-approksimaatioksi saadaan tämän perusteella \[ \theta | y \sim N_m(\hat{\theta},V). \] On siis laskettava funktion \(h(\theta)\) maksimi (joko derivoimalla tai numeerisella optimoinnilla) sekä sen toiset derivaatat maksimikohdassa. On huomattava, että tämä on approksimaatio. Jos normeeraustekijää tarvitaan, niin se voidaan myös laskea \[ \begin{aligned} p(y) &\approx \int \exp(h(\hat\theta)) \exp \left( -\frac{1}{2}(\theta-\hat{\theta})^T V^{-1} (\theta-\hat{\theta}) \right)\, d\theta\\ &= p(\hat\theta) p(y|\hat\theta) (2\pi)^{\frac{m}2} \det(V)^{\frac12} \times \int (2\pi)^{-\frac{m}2}\det(V)^{-\frac12} \exp \left( -\frac{1}{2}(\theta-\hat{\theta})^T V^{-1} (\theta-\hat{\theta}) \right)\, d\theta \\ &= p(\hat\theta) p(y|\hat\theta) (2\pi)^{\frac{m}2} \det(-h''(\hat\theta))^{-\frac12}. \end{aligned} \] Approksimaatio voidaan kytkeä Newtonin algoritmiin:
Jos algoritmi konvergoi kun \(k \to \infty\), niin silloin
Esimerkki 7.2 (Genetiikkaesimerkki (jatk.)) Palataan genetiikkaesimerkkiin, jossa \[ \begin{aligned} p(\theta|y) &\propto \left(\frac12 + \frac{\theta}4 \right)^{y_1}\left(\frac14(1-\theta)\right)^{y_2} \left(\frac14(1-\theta)\right)^{y_3} \left(\frac{\theta}4\right)^{y_4} \\ &\propto (2 + \theta)^{y_1}(1-\theta)^{y_2 + y_3}\theta^{y_4}. \end{aligned} \] Nyt \[ \begin{aligned} h(\theta) &= y_1 \log(2 + \theta) + (y_2 + y_3)\log(1 - \theta) + y_4\log(\theta) \\ h'(\theta) &= \frac{y_1}{2+\theta} - \frac{y_2 + y_3}{1 - \theta} + \frac{y_4}{\theta} \\ h''(\theta) &= -\frac{y_1}{2 + \theta} - \frac{y_2 + y_3}{(1 - \theta)^2} - \frac{y_4}{\theta^2}. \end{aligned} \] Newtonin algoritmilla saadaan \(\hat\theta = 0.6268215\) ja \(h''(\hat\theta) = -377.5169.\) Posteriori on approksimatiivisesti normaalijakauma, jonka odotusarvo on \(0.627\) ja varianssi \(0.0026\).
Posteriorista voidaan laskea esimerkiksi approksimatiivinen 90 %:n Bayes-väli parametrille \(\theta\): \(0.627 \pm 1.644 \times 0.0516\) ja saadaan tulokseksi \((0.590, 0.711)\). Tulos vastaa varsin hyvin numeerisella integroinnilla saatua (jota voidaan pitää tarkkana). Tässä tapauksessa Laplace-approksimaatio toimii siis hyvin.
Normaalijakauman sijaan posteriorijakaumaa voi approksimoida myös jollain muulla jakaumalla. Bayes-tilastotieteessä voidaan käyttää variaatiolaskentaa (variational inference) parhaan approksimaation löytämiseen. Approksimaation ominaisuuksia ei vielä tunneta perinpohjaisesti.
Menetelmässä posteriorijakaumaa \(p(\theta | y)\) approksimoidaan parametrisella jakaumalla \(q_\eta(\theta)\). Useimmiten käytetään tulomuotoisia jakaumia \[ q_\eta(\theta) = \prod_{i=1}^m q_\eta^{(i)}(\theta_i) \] jollekin parametrivektorin ositukselle \(\{\theta_1,\ldots,\theta_m\}\). Approksimaatioon liittyvät parametrit \(\eta\) määritetään siten, että approksimaation \(q_\eta(\theta)\) etäisyys posteriorista \(p(\theta | y)\) minimoituu. Jakaumien välistä etäisyyttä mitataan Kullback–Leibler-divergenssilla \[ \int q_\eta(\theta) \log \left( \frac{q_\eta(\theta)}{p(\theta | y)}\right)\, d\theta. \] Tulomuotoisen approksimaation tapauksessa variaatiolaskennan avulla voidaan osoittaa, että etäisyys minimoituu, kun
\[ \log q_\eta^{(i)}(\theta_i) \propto \int q_\eta^{(-i)}(\theta_{-i}) \log(p(y,\theta))\, d\theta_{-i}, \tag{7.1}\] missä alaindeksi \(-i\) tarkoittaa ositteen \(i\) poisjättämistä. Ratkaisu 7.1 tulee normalisoida jakaumaksi. Käytännössä laskenta vaatii usein iterointia, koska ratkaisu ositteelle \(\theta_i\) riippuu muiden ositteiden ratkaisuista. Yleisemmin Kullback–Leibler-divergenssi voidaan Bayesin kaavan mukaan kirjoittaa muodossa \[ \int q_\eta(\theta) \log \left( \frac{q_\eta(\theta)}{p(\theta | y)}\right)\, d\theta = -\int q_\eta(\theta) \log p(y | \theta) \,d\theta + \int q_\eta(\theta) \log \left( \frac{q_\eta(\theta)}{p(\theta)}\right)\, d\theta + \log p(y), \] mistä nähdään, että optimaalinen variationaalinen approksimaatio \(q_\eta(\theta)\) tasapainottaa uskottavuusfunktion maksimoimista ja pientä etäisyyttä priorijakaumasta. Variaatioapproksimaation parametrit \(\eta\) voidaan optimoida esim. stokastista gradienttimenetelmää käyttäen.
Integroitu sisäkkäinen Laplace-approksimaatio (Integrated Nested Laplace Approximation, INLA) on Laplace-approksimaation yleistys. Menetelmää voi käyttää latenteille gaussisille malleille, joita ovat muun muassa yleistetyt lineaariset mallit sekä monet spatiaaliset mallit. INLA on käyttökelpoinen menetelmä silloin, kun hyperparametrien määrä \(\theta\) on pieni, tyypillisesti alle 6. Latenttien muuttujien määrä sen sijaan saa olla suuri.
Latentti gaussinen malli koostuu kahdesta osasta: latentit muuttujat \(x\) ovat normaalijakautuneita mallin \(p(x | \theta)\) mukaisesti ja havainnot \(y\) noudattavat ehdollista jakaumaa \(p(y | x,\theta)\). Keskeinen idea on kirjoittaa hyperparametrien posteriori muodossa \[ p(\theta | y) \propto \frac{p(\theta)p(x | \theta)p(y | x, \theta)}{p(x | \theta,y)} \tag{7.2}\] ja käyttää Laplace-approksimaatiota nimittäjälle \(p(x | \theta, y)\). Approksimaatio on varsin tarkka, koska jakauma \(p(x | \theta)\) on normaalijakauma. Kaava 7.2 tuottaa approksimaation normeeraamattomalle posteriorille \(p(\theta)p(y | \theta)\). Lisää approksimaatioita tarvitaan, kun halutaan määrittää marginaaliposteriorijakaumat parametrivektorin \(\theta\) komponenteille. Nämä approksimaatiot ovat teknisesti hankalia ja toimivat ainoastaan silloin, kun parametrivektorin \(\theta\) dimensio on pieni.
Lisätietoja menetelmästä on saatavilla INLAn kotisivulta (https://www.r-inla.org/home).
Hylkäysotannassa (rejection sampling) pyritään tuottamaan havaintoja kiinnostuksen kohteena olevasta posteriorijakaumasta \(p(\theta|y)\) käyttämällä hyväksi ehdotusjakaumaa \(g(\theta)\). Hylkäysotanta edellyttää, että seuraavat ominaisuudet ovat voimassa ehdotusjakaumalle:
Havaintoja voidaan tuottaa jakaumasta \(g(\theta)\).
Tärkeyssuhteelle on voimassa \(p(\theta|y)/g(\theta) \leq M\) kaikilla \(\theta\), missä \(M\) on jokin positiivinen vakio.
Hylkäysotanta etenee seuraavasti (yhden realisaation tuottamiseksi):
Generoi ehdotus \(\theta^*\) jakaumasta \(g(\theta)\).
Hyväksy \(\theta^*\) realisaatioksi jakaumasta \(p(\theta|y)\) todennäköisyydellä \(p(\theta|y)/(Mg(\theta))\). Jos ehdotus hylätään, palaa kohtaan 1.
Tärkeyssuhteelle voimassaoleva ehto \(p(\theta|y)/g(\theta) \leq M\) takaa sen, että ehdotuksen hyväksymistodennäköisyyden lauseke on aina \(\leq 1\) kaikilla \(\theta\):n arvoilla. Menetelmän toimivuutta voi arvioida hyväksyttyjen realisaatioiden osuuden avulla. Hylkäysotanta toimii sitä tehokkaammin, mitä enemmän \(g(\theta)\) muistuttaa posterioria \(p(\theta|y)\).
Tärkeysotannassa (importance sampling) jostakin sopivasta ehdotusjakaumasta generoituja havaintoja painotetaan siten, että ne edustavat kiinnostuksen kohteena olevaa posteriorijakaumaa. Ehdotusjakauman \(g(\theta)\) tulisi muistuttaa posteriorijakaumaa ja havaintojen generoimisen tulisi olla helppoa.
Olkoon kiinnostuksen kohteena odotusarvo \(\operatorname{E}(\xi(\theta) | y)\), missä \(\xi\) on jokin tunnettu funktio. Tärkeysotanta perustuu tulokseen \[ \begin{aligned} \operatorname{E}(\xi(\theta) | y) &=\int \xi(\theta) p(\theta | y)\, d\theta = \frac{\int \xi(\theta)p(\theta)p(y | \theta)\, d\theta}{\int p(\theta)p(y | \theta)\, d\theta} \\ &= \frac{\int \xi(\theta)\frac{p(\theta)p(y | \theta)}{g(\theta)}g(\theta)\, d\theta}{\int \frac{p(\theta)p(y | \theta)}{g(\theta)} g(\theta)\, d\theta} = \frac{\operatorname{E}_g(w(\theta)\xi(\theta))}{\operatorname{E}_g(w(\theta))}, \end{aligned} \] missä painot määräytyvät kaavasta \(w(\theta) = p(\theta)p(y | \theta)/g(\theta)\).
Generoidaan nyt riippumattomia havaintoja \(\theta^{(1)},\dots, \theta^{(n)}\) jakaumasta \(g(\theta)\). Keskiarvon \(\overline \xi=\operatorname{E}(\xi(\theta)|y)\) tärkeysotantaan perustuva estimaattori on painotettu keskiarvo \[ \overline \xi_{IS}=\frac{ \sum_{j=1}^n w(\theta^{(j)})\xi(\theta^{(j)}) }{ \sum_{j=1}^n w(\theta^{(j)}) }, \] ja estimoitu keskivirhe saadaan painotettujen havaintojen varianssin avulla \[ \operatorname{SE}(\overline \xi_{IS})=\frac {\sqrt{\sum_{j=1}^n [w(\theta^{(j)})]^2\, (\xi(\theta^{(j)})-\overline \xi_{IS})^2}} {\sum_{j=1}^n w(\theta^{(j)})}\,. \] Jos ehdotusjakauma on valittu huonosti, suurin osa painoista voi olla lähellä nollaa samalla, kun muutamat painot ovat todella suuria. Tällöin tärkeysotanta toimii huonosti. Ehdotusjakauman voi valita esimerkiksi Laplace-approksimaation pohjalta. Joissain tapauksissa varianssin suurentaminen saattaa olla tarpeen liian suurien painojen välttämiseksi. Ehdotusjakaumana voi käyttää myös informatiivista prioria. Tällöin painot yksinkertaistuvat uskottavuudeksi \[ w(\theta)=\frac{p(\theta)p(y | \theta)}{g(\theta)}=\frac{p(\theta)p(y | \theta)}{p(\theta)}=p(y | \theta). \]
Esimerkki 7.3 (Genetiikkaesimerkki (jatk.)) Lasketaan posterioriodotusarvo tärkeysotannalla käyttäen simulointijakaumana Laplacen approksimaatiota \(\operatorname{N}(0.6268215, 0.0514673^2)\):
importance_sampling <- function(y, mu, sigma, n_sim) {
theta <- rnorm(n_sim, mu, sigma)
dens <- dnorm(theta, mu, sigma)
p_un <- (0.5 + 0.25 * theta)^y[1] * (0.25 * (1 - theta))^y[2] *
(0.25 * (1 - theta))^y[3] * (0.25 * theta)^y[4]
w <- p_un / dens
theta_is <- sum(w * theta) / sum(w)
theta_is
}
mu <- 0.6268215
sigma <- 0.0514673
y <- c(125, 18, 20, 34)
importance_sampling(y, mu, sigma, 1e4)[1] 0.622691Mikä on hyvin lähellä numeerisella integroinnilla saatua arvoa.