Metropolisin ja Hastingsin algoritmi

Metropolisin ja Hastingsin algoritmi on Markovin ketjuun perustuva Monte Carlo -algoritmi. Se on näytteistysalgoritmi, eli sen avulla voidaan generoida likimäärin tietyn todennäköisyysjakauman mukaisesti jakautunut luku- tai vektorijoukko. Generoitua joukkoa voidaan käyttää esimerkiksi Monte Carlo -integroinnissa.

Yleistä

Metropolisin ja Hastingsin algoritmilla voidaan periaatteessa generoida näytteitä mistä tahansa jakaumasta, kunhan jakauman tiheysfunktion arvot voidaan laskea missä tahansa tila-avaruuden pisteessä. Tiheysfunktio saa myös olla normalisoimaton eli sen integraalin (jatkuvan tila-avaruuden tapauksessa) tai summan (diskreetin tila-avaruuden tapauksessa) koko tila-avaruuden yli ei tarvitse olla yksi. Algoritmi on erityisen hyödyllinen bayesläisessä tilastotieteessä, jossa esiintyy usein analyyttiseltä muodoltaan hyvin monimutkaisia todennäköisyysjakaumia.

Hyvin tunnettu Gibbsin näytteistysalgoritmi on erikoistapaus Metropolisin ja Hastingsin algoritmista. Se on laskennallisesti tehokas mutta vaatii jakaumalta tiettyjä lisäominaisuuksia.

Algoritmin kuvaus

Simuloitavaa jakaumaa kutsutaan kohdejakaumaksi ja merkitään $\pi$. Se saa siis olla normalisoimaton. Algoritmi tuottaa realisoituman sellaiselle Markovin ketjulle, jonka invariantti jakauma kohdejakauma on. Ketjun generoimiseksi on generoitava toinen Markovin ketju, ja sen siirtymätodennäköisyyttä pisteestä $x$ pisteeseen $x'$ merkitään tässä $q(x'|x)$. Funktio $q(x'|x)$ on tiheysfunktio muuttujan $x'$ suhteen, ja sitä kutsutaan ehdokasjakaumaksi. Tämän ketjun määrittelevä ehdokasjakauma voidaan valita mielivaltaisesti, kunhan siitä voidaan (jollakin tähän algoritmiin liittymättömällä tavalla) generoida näytteitä kaikilla nykyisen tilan $x$ arvoilla ja kunhan ehdokasjakauman kantaja sisältää kohdejakauman kantajan. Ehdokasjakauman valinta tosin vaikuttaa huomattavasti algoritmin laskennalliseen tehokkuuteen.

Metropolisin ja Hastingsin algoritmi tarvittavan Markovin ketjun generoimiseen on seuraava: Jos ollaan ketjun alussa, ketjulle on määrättävä alkuarvo sellaiseen pisteeseen, että simuloitava tiheys on siinä nollasta poikkeava. Näin ollen ketjulla on olemassa jokin tila. Olkoon nyt $x$ ketjun nykyinen tila. Generoidaan ehdokasjakaumasta $q(\cdot |x)$ satunnaisluku, jota kutsutaan ehdokasarvoksi ja merkitään $x'$. Kutsutaan nyt Metropolisin ja Hastingsin suhteeksi lukua

$r(x,x')={\frac {\pi (x')q(x|x')}{\pi (x)q(x'|x)}}.$

Seuraavaksi generoidaan satunnaisluku $u$ tasajakaumasta $\operatorname {Tas} (0,1)$. Jos pätee $r(x,x')>u$, niin sanotaan, että ehdokasarvo hyväksytään, ja asetetaan ketjun seuraavaksi tilaksi ehdokasarvo $x'$. Muutoin ehdokasarvo hylätään ja asetetaan ketjun seuraavaksi tilaksi nykyisen tilan arvo $x$. Ehdokasarvon hyväksymisen todennäköisyys on siis

$\alpha (x,x')=\min \left\{{\frac {\pi (x')q(x|x')}{\pi (x)q(x'|x)}},1\right\}.$

Koska ketjun alkupää saattaa riippua voimakkaasti ketjun alkuarvosta, alusta on yleensä syytä poistaa joitakin arvoja. Alkuosan poikkeamista muun ketjun jakaumasta kutsutaan burn in -ilmiöksi. Näin muodostetun ketjun tiloista muodostuu vaaditut ehdot täyttävän Markovin ketjun realisoituma.

Ehdokasjakaumasta

Vaikka ehdokasjakaumalla ei teoriassa juuri ole rajoitteita, algoritmin suorituskyvyn kannalta ehdokasjakauman valinta voi olla hyvin kriittinen. Tavallisia ovat esimerkiksi riippumattomaan satunnaisjonoon sekä satunnaiskulkuun perustuvat Metropolisin ja Hastingsin algoritmit.

Satunnaiskulkuun perustuvassa Metropolisin ja Hastingsin algoritmissa ehdokasarvo generoidaan lausekkeella

$x'=x+w,$

missä $w$ on generoitu jostakin ketjun tilasta riippumattomasta jakaumasta. Tavallisin on normaali satunnaiskulku, jossa luvut $w$ generoidaan nollakeskisestä normaalijakaumasta. Tällöin ehdokasketju on symmetrinen ehdon $q(x|x')=q(x'|x)$ mielessä, joten Metropolisin ja Hastingsin suhde on

$r(x,x')={\frac {\pi (x')}{\pi (x)}}.$

Algoritmin testaaminen

Käytännössä ei koskaan voida varmasti todeta, että jokin Metropolisin ja Hastingsin algoritmin tuottama Markovin ketju simuloi hyvin kohdejakaumaa. Usein luotettavin keino on näytehistorioitten silmämääräinen tarkastelu. Näytehistoria on kuvaaja ketjun tilan arvoista tilan indeksien ("ajanhetkien") funktiona.

Satunnaiskulun tapauksessa keskeinen algoritmin tehokkuutta säätelevä parametri on ehdokasjakauman varianssi $\sigma ^{2}$. Jos varianssi on liian suuri, algoritmi on tehoton, koska hyväksymistodennäköisyys on alhainen. Jos taas varianssi on liian pieni, ehdokasarvot ovat usein lähellä nykyisiä harvoja. Tällöin ketjusta tulee voimakkaasti autokorreloitunut eli ketjuun on generoitava hyvin monta arvoa, jotta se simuloisi hyvin koko jakaumaa.

Taustalla olevaa teoriaa

Yleisesti näytteistysalgoritmin toimimisen takaavat lauseet ovat suurten lukujen laki je keskeinen raja-arvolause.

Metropolisin ja Hastingsin algoritmin tuottamat Markovin ketjun tilat eivät yleisesti ole toisistaan riippumattomia, joten tavallisesti esitetyt versiot suurten lukujen laista ja keskeisestä raja-arvolauseesta eivät riitä.

Artikkelissa todistetaan Metropolisin ja Hastingsin algoritmille oma suurten lukujen lause. Sen mukaan jos ehdokasjakauman kantaja sisältää kohdejakauman kantajan, niin algoritmin tuottaman otoksen keskiarvo lähestyy kohdejakauman odotusarvoa, kun otoskoko lähestyy ääretöntä. Sama koskee kaikkia muitakin jakauman momentteja, siis esimerkiksi varianssia.

Samoin artikkelissa todistetaan keskeinen raja-arvolause Metropolisin ja Hastingsin algoritmille. Siinä otoskeskiarvon poikkeama odotusarvosta riippuu otoskoon lisäksi esimerkiksi ketjun autokorrelaatiofunktiosta, joten se ei tarjoa yhtä käytännöllistä virhearviota kuin riippumattomien jonojen keskeinen raja-arvolause. Se kuitenkin osoittaa, että otoskesarvo on suurilla otosko'oilla likimain normaalijakautunut ja että sen poikkeama jakauman odotusarvosta eli odotusarvoestimaatin virhe on $1/{\sqrt {N}}$-riippuvainen, missä $N$ on ketjun pituus.

Aiheesta muualla

• Metropolis-Hastings algorithm on xβ
• Matlab implementation of Metropolis-Hastings (Arkistoitu – Internet Archive)