Uvažujte situaci z úlohy 5 ze čtvrtého listu ke cvičením, tj. cílem je sestavit MH-algoritmus s nezávislými návrhy pro simulaci z N(0,1).
Napište funkci, která bude generovat markovský řetězec s cílovým
rozdělením s hustotou N(0,1). Návrhovou hustotu volte:
a) normální rozdělení s nulovou střední hodnotou a rozptylem
σ2,
b) Laplaceovo rozdělení s parametrem a.
Pozorujte, jaký je vliv rozptylu návrhového rozdělení na průběh vygenerovaného řetězce a pravděpodobnost přijetí.
Zkuste také použít algoritmus symetrické náhodné procházky, opět s různými volbami návrhového rozdělení (rovnoměrné, normální) a různými startovacími hodnotami.
Chování řetězce - tedy mixování a rychlost konvergence řetězce k cílovému rozdělení souvisí s ergodicitou použitého MCMC algoritmu. Prakticky se parametry algoritmu vybírají po prozkoumání pilotního běhu algoritmu - zkoumá se např. průměrná přijímací pravděpodobnost, autokorelace v řetězci atp.
Nicméně, i když generované řetězce vypadají OK, ještě to neznamená, že máme dobře volený MCMC algoritmus, který rozumně rychle konverguje k cílovému rozdělení a pro ergodické průměry platí CLV. Vždy je dobré znát teoretické vlastnosti použitého MCMC řetězce abychom se vyvarovali nejhoršího. Např. geometrická ergodicita nám zajišťuje CLV pro ergodické průměry. Když neplatí, ani CLV pro ergodické průměry nemusí platit (čímž se algoritmus stává zcela nepoužitelným !)
Nyní chceme sestavit algoritmus pro simulování z exponenciálního rozdělení s parametrem 1 a opět použít MH algoritmus s nezávislými návrhy. Jako návrhové rozdělení volíme opět exponenciální s parametrem k.
S využitím předchozího lehce řešíme např. takto
Zkuste volit k=0.01 a k=5. Které návrhové rozdělení vypadá jako vhodnější (podle průběhu řetězce)?
Přestože by se mohlo zdát, že první varianta algoritmu mixuje mnohem pomaleji a má malé pravděpodobnosti přijetí (a výsledný řetězec velké hodnoty autokorelační funkce i pro velké rozsahy), tak problémový je algoritmus pro k=5.
Že je něco špatně je vidět na obrázku s jádrovým odhadem rozdělení ergodického průměru přes 1000000 iterací, získaný pomocí 55 nezávislých běhů algoritmu pro k=0.01 a k=5, startovaných z 1.
Pro k=5 to nevypadá, že by ergodické průměry konvergovaly k normálnímu rozdělení. A opravdu nekonvergují - CLV tu neplatí (Roberts(96)) (neplatí pro žádné k větší nebo rovno 2). Pro k=0.01 platí (Roberts a Rosenthal(97)).
Také MH algoritmus je geometricky ergodický pro k=0.01, ale pro k>2 není geometricky ergodický (Smith a Tierney(96)) (algoritmus je geometricky ergodický pro k menší nebo rovné 1). Nicméně pro k=0.01 je řád konvergence velmi velký, cca. 0.99. Proto máme v řetězci velké autokorelace a pomalé mixování.
Pro porovnání dva různé běhy algoritmu pro k=5 (počet iterací = 1000000) a odhadnuté empirické autokorelační funkce.
Tady je vidět, v čem spočívá problém algoritmu pro k=2. Většina běhů algoritmu bude vypadat cca jako druhý průběh, bude se držet v oblasti menších hodnot a výsledný ergodický průměr bude značně podhodnocovat teoretickou střední hodnotu 1. Nicméně občas nějaký průběh bude vypadat jako průběh a - algoritmus se dostane (i když je to málo pravděpodobné) do velké hodnoty (míněno cca větší než 3) a díky divergenci vah pro velké hodnoty výchozího stavu se v této hodnotě "zasekne" na velký počet iterací, než se mu podaří nezamítnout návrh a tuto pozici opět opustit. Výsledný ergodický průměr přes běh potom bude znatelně větší než teoretická hodnota 1.
Několik různých variací na Metropolis-Hastingsův algoritmus pro generování z různých jednorozměrných hustot je naimplementováno zde. A to včetně online průběhu algoritmu, histogramů a odhadu autokorelační funkce a možnosti měnit parametry používaných hustot.