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.
Pochopitelně v praxi bychom raději použili nějakou přímou metodu, v R např. funkci mvrnorm z knihovny MASS.
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)
.
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.
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)
.
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.