Cvičení z MCMC 14.11.2018 - Gibbsův výběrový plán


Jednoduchý ilustrační příklad

Nechť náhodný vektor (X,Y) má dvourozměrné normální rozdělení s nulovou střední hodnotou. Potom platí, že podmíněná rozdělení X|Y a Y|X jsou normální. Můžeme tedy jednoduše sestavit Gibbsův výběrový plán pro simulaci ze sdruženého rozdělení. Vyzkoušejte různé volby počátečních stavů a varianční matice.

norm-gibbs.R

Pochopitelně v praxi bychom raději použili nějakou přímou metodu, v R např. funkci mvrnorm z knihovny MASS.


Domácí úloha - počet důlních katastrof

Odhadnutí change pointu v Poissonově procesu pro data .
Data (Jarret, 1979) představují časovou řadu s počty důlních katastrof za rok ve Velké Británii v průběhu let 1851-1962. Pomocí Gibbsova výběrového plánu odhadneme aposteriorní rozdělení parametrů modelu.

Možné řešení

Zjistili jsme, že s velkou pravděpodobností je bod zvratu kolem roku 1891 a že intenzity Poissonova rozdělení jsou dost odlišné před a za tímto bodem. Střední hodnota, modus a medián spočtené ze skutečného aposteriorního rozdělení jsou shrnuty v tabulce. Rovněž jsou uvedeny 95% oboustranné aposteriorní intervaly spolehlivosti (v diskrétním případě přibližně 95%).

Aposteriorní střední hodnota MAP (Maximum aposterior) Aposteriorní medián 95% interval spolehlivosti
Change point 1890 1891 1889 (1890) (1886, 1895)
k 39.95 41 39 (40) (36,45)
θ 3.12 3.09 3.11 (2.60, 3.69)
λ 0.92 0.91 0.92 (0.71,1.16)

MCMC aproximace (pro 5000 iterací) vycházejí přibližně kolem těchto hodnot.

Mohla by nás zajímat otázka predikce počtu katastrof v dalších letech. Prediktivní rozdělení dalšího pozorování z při daných minulých pozorování x je dáno pravděpodobnostmi
f(z|x) = ∫ f(z|θ,λ,k,x) π(θ,λ,k|x) ν(d(θ,λ,k)).
Protože máme k dispozici realizaci řetězce, jehož stacionární rozdělení je π(θ,λ,k|x), a f(z|θ,λ,k,x) = exp(-λ) λz/z!, můžeme f(z|x) odhadnout jako
(1/n) ∑ f(z|λi).

gibbs-du-ext.R

V tomto případě malá variabilita v aposteriorním rozdělení λ znamená, že odhad prediktivního rozdělení a odhad založený pouze na aposteriorní střední hodnotě jsou skoro stejné. Obecně však prediktivní rozdělení má větší variabilitu, protože bere do úvahy nejistotu o parametru λ.

Vhodný software pro provádění MCMC v bayesovských modelech je BUGS (Bayesian Inference using Gibbs Sampling) nebo JAGS (Just Another Gibbs Sampler). Oba je možné spouštět i z R.


Domácí úloha - Gibbsův algoritmus pro "rozdělení" s exponenciálními úplnými podmíněnými hustotami

V příkladě P2 z třetího listu úloh jsme diskutovali Gibbsův výběrový plán pro "dvourozměrné rozdělení" s podmíněnými rozděleními složek X|Y ~ Exp(Y) a Y|X ~ Exp(X).

gibbs-exp.R

ukazuje chování Gibsova algoritmu pro toto zadání i pro rozdělení náhodného vektoru se stejnými podmíněnými rozděleními, ale omezenými jen na interval [0,B], B >0.