Metropolis Hasting Algorithm:

?

MH算法也是一種基于模擬的MCMC技術，一個很重要的應用是從給定的概率分布中抽樣。主要原理是構造了一個精妙的Markov鏈，使得該鏈的穩(wěn)態(tài)是你給定的概率密度。它的好處，不用多說，自然是可以對付數(shù)學形式復雜的概率密度。有人說，單維的MH算法配上Gibbs Sampler幾乎是“無敵”了。

?

今天試驗的過程中發(fā)現(xiàn)，MH算法想用好也還不簡單，里面的轉移參數(shù)設定就不是很好弄。即使用最簡單的高斯漂移項，方差的確定也是個頭疼的問題，需要不同問題不同對待，多試驗幾次。當然你也可以始終選擇“理想”參數(shù)。

?

還是拿上次的混合高斯分布來做模擬，模擬次數(shù)為500000次的時候，概率分布逼近的程度如下圖。雖然幾個明顯的"峰"已經(jīng)出來了，但是數(shù)值上還是有很大差異的。估計是我的漂移方差沒有選好。感覺還是inverse sampling好用，迭代次數(shù)不用很多，就可以達到相當?shù)谋平潭取?/p>

?

?

試了一下MH算法，

?

?

?

R Code：

p=function(x,u1,sig1,u2,sig2){
(1/3)*(1/(sqrt(2*pi)*15)*exp(-0.5*(x-70)^2/15^2)+1/(sqrt(2*pi)*11)*exp(-0.5*(x+80)^2/11^2)+1/(sqrt(2*pi)*sig1)*exp(-0.5*(x-u1)^2/sig1^2)+1/(sqrt(2*pi)*sig2)*exp(-0.5*(x-u2)^2/sig2^2))
}


MH=function(x0,n){
x=NULL
x[1] = x0
for (i in 1:n){
? x_can= x[i]+rnorm(1,0,3.25)
? d= p(x_can,10,30,-10,10)/p(x[i],10,30,-10,10)
? alpha= min(1,d)
? u=runif(1,0,1)
??? if (u<alpha){
??? x[i+1]=x_can}
??? else{
????? x[i+1]=x[i]
???? }
?? if (round(i/100)==i/100) print(i)
}
x
}
z=MH(10,99999)
z=z[-10000]
a=seq(-100,100,0.2)

plot(density(z),col=1,main='Estimated Density',ylim=c(0,0.02),lty=1)
points(a, p(a,10,30,-10,10),pch='.',col=2,lty=2)
legend(60,0.02,c("True","Sim (MH)"),col=c(1,2),lty=c(1,2))

轉載于:https://www.cnblogs.com/juggernaunt/archive/2011/01/16/1936609.html

總結

以上是生活随笔為你收集整理的Metropolis Hasting算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。