Laskennallisen tekniikan laboratorio: Opetus: S-114.600
TKK / Yksiköt / LCE / Opetus / Kurssit / S-114.600
Metropolis- ja Gibbs-otanta
a) Metropolis-otanta.
Toteuta Metropolis-otanta kaksiulotteiselle normaalijakaumalle. Katso Bayesian Data Analysis s. 289.
Ehdotusjakauma q on tässä tapauksessa normaalijakauma, jonka keskipiste on tt_i ja kovarianssimatriisi on diagonaalinen, eli normaalijakauma on isotrooppinen, eli yhtä leveä joka suuntaan. Käytä ehdotusjakaumana normaalijakaumaa jonka skaala on sp. Jos käytät mvnrnd-funktiota, huomaa, että se ottaa argumenttinaan kovarianssimatriisin. Kun sp on ehdotusjakauman hajonta, saat halutun kovarianssimatriisin Matlabissa näin: sp^2*eye(2). Diagonaalisen kovarianssin tapauksessa muuttujat ovat riippumattomia, joten voit arpoa ne myös esim. näin randn(1,2)*sp.
Huomaa, että vaikka sekä kiinnostuksen kohteena oleva jakauma, että ehdotusjakauma tässä ovat molemmat normaalijakaumia näin ei tarvitse olla. Varmista, että ymmärät eron ehdotusjakauman ja kiinnostuksen kohteena olevan jakauman välillä.
Kiinnostuksen kohteena oleva jakauma on normaalijakauma, jonka kovarianssimatriisi on S=[1 r;r 1]. Alussa r on 0.8 ja sitten voit kokeilla miten käy jos r kasavaa, eli muuttujien välinen korrelaatio kasvaa. Leikimme nyt, että emme tästä jakaumasta osaa suoraan poimia näytteitä, mutta osaamme kuitenkin laskea todennäköisyystiheyden halutussa pisteessä (tämä on se yleinen tilanne hankalammilla jakaumilla kuin normaalijakauma). mnorm_lpdf(ttnew,[y1 y2],S) laskee moniulotteisen normaalijakauman todennäköisyystiheyden logaritmin.
Todennäköisyystiheydet (pdf) on hyvä laskea logaritmeina lukutarkkuuden vuoksi ja myös laskea laskut riittävän pitkälle logaritmeina. Matlabin omat pdf-funktiotkin laskevat sisäisesti logaritmeilla, mutta lopuksi palauttavat tuloksen exp:in läpi. Tämän vuoksi olen tarjonnut kaikki tarvittavat log-pdf-funktiot käyttöönne, jotta laskut voitaisiin laskea riittävän pitkälle logaritmeina ilman, että välillä otetaan exp.
Metropolis-algoritmin oleellisena osana on termi alpha=min(1,a/b). Tämä voidaan laskea suoraan näin, mutta parempaan lukutarkkuuteen päästään jos tämä lasketaan alpha=min(1,exp(log(a)-log(b))). Ja nyt lpdf-funktioiden avulla saadaan siis suoraan laskettua log(a) ja log(b). Tyypillisesti a ja b voivat olla hyvin isoja tai hyvin pieniä lukuja ja silloin niiden välinen jakolasku ei ole tarkka (yli/alivuodot mahdollisia). Jos on suoraan laskettu log(a) ja log(b) niin tarkempi tulee kun lasketaan logaritimien vähennyslasku. Tämän jälkeen voidaan kuitenkin siirtyä pois logaritmeista, koska alpha on luku 0 ja 1 välissä, ja algoritmin kannalta ei ole väliä onko alpha kovin tarkka lähellä 0 tai 1 suuremilla arvoilla.
Kokeile eri r:n arvoilla, eri ketjun alkupisteillä, eri pituisilla ketjuilla ja eri kokoisilla ehdotusjakaumilla (esim. 0.1<sp<4). Työn palautukseen liitä mukaan algoritmin toteutus, esimerkkikuva ja lyhyet vastaukset seuraaviin kysymyksiin:
Miten autokorrelaation (acorr,acorrtime) käy r:n kasvaessa?
Miten alkuarvaus vaikuttaa burn-in aikaan?
Miten ehdotusjakauman skaala sp vaikuttaa hylkäystaajuteen, burn-in aikaan ja autokorrelaatioon?
sampling1.m sisältää alkutiedot ja koodin esimerkkikuvaa varten. Lisää vain itse algoritmin toteuttava koodi (alle 15 riviä).
b) Gibbs-otanta.
Toteuta Gibbs-otanta kaksiulotteiselle normaalijakaumalle. Katso Bayesian Data Analysis s. 287.
Gibbs-otantaa varten tarvitset ehdolliset jakaumat (s. 288) Huomaa, että 1-r^2 on normaalijakauman varianssi, ja että esim. normrnd ottaa toisena argumenttinaan hajonnan sqrt(1-r^2).
Huomaa, että Metropolis-algoritmissa riitti, että osasimme laskea tiheysfunktion arvon halutussa pisteessä. Gibbs-otannassa on osattava pomia näytteitä jakauman ehdollisista jakaumista.
Kokeile eri r:n arvoilla, eri ketjun alkupisteillä ja eri pituisilla ketjuilla. Työn palautukseen liitä mukaan algoritmin toteutus, esimerkkikuva ja lyhyet vastaukset seuraaviin kysmyksiin:
sampling2.m sisältää alkutiedot ja koodin esimerkkikuvaa varten. Lisää vain itse algoritmin toteuttava koodi (alle 10 riviä).
Tarvitset toteutuksessa ja esimerkkikuvien teossa seuraavia funktioita, jotka eivät tule Matlabin mukana:
ellipse.m
acorr.m
acorrtime.m
mnorm_lpdf.m
Sivua on viimeksi päivitetty 26.8.2003
URL: http://www.lce.hut.fi/teaching/S-114.600/ex/sampling.html