原文鏈接：http://tecdat.cn/?p=23019?

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

蒙特卡洛方法利用隨機(jī)數(shù)從概率分布P(x)中生成樣本，并從該分布中評(píng)估期望值，該期望值通常很復(fù)雜，不能用精確方法評(píng)估。在貝葉斯推理中，P（x）通常是定義在一組隨機(jī)變量上的聯(lián)合后驗(yàn)分布。然而，從這個(gè)分布中獲得獨(dú)立樣本并不容易，這取決于取樣空間的維度。因此，我們需要借助更復(fù)雜的蒙特卡洛方法來(lái)幫助簡(jiǎn)化這個(gè)問(wèn)題；例如，重要性抽樣、拒絕抽樣、吉布斯抽樣和Metropolis Hastings抽樣。這些方法通常涉及從建議密度Q(x)中取樣，以代替P(x)。

在重要性抽樣中，我們從Q(x)中產(chǎn)生樣本，并引入權(quán)重以考慮從不正確的分布中抽樣。然后，我們對(duì)我們需要評(píng)估的估計(jì)器中的每個(gè)點(diǎn)的重要性進(jìn)行調(diào)整。在拒絕抽樣中，我們從提議分布Q(x)中抽取一個(gè)點(diǎn)，并計(jì)算出P(x)/Q(x)的比率。然后我們從U(0,1)分布中抽取一個(gè)隨機(jī)數(shù)u；如果

，我們就接受這個(gè)點(diǎn)x，否則就拒絕并回到Q(x)中抽取另一個(gè)點(diǎn)。吉布斯抽樣是一種從至少兩個(gè)維度的分布中抽樣的方法。這里，提議分布Q(x)是以聯(lián)合分布P(x)的條件分布來(lái)定義的。我們通過(guò)從后驗(yàn)條件中迭代抽樣來(lái)模擬P(x)的后驗(yàn)樣本，同時(shí)將其他變量設(shè)置在其當(dāng)前值。

雖然，重要性抽樣和拒絕抽樣需要Q(x)與P(x)相似（在高維問(wèn)題中很難創(chuàng)建這樣的密度），但當(dāng)條件后驗(yàn)沒(méi)有已知形式時(shí)，吉布斯抽樣很難應(yīng)用。這一假設(shè)在更普遍的Metropolis-Hastings算法中可以放寬，在該算法中，候選樣本被概率性地接受或拒絕。這種算法可以容納對(duì)稱(chēng)和不對(duì)稱(chēng)的提議分布。該算法可以描述如下?

初始化

抽取

計(jì)算

從

中抽取

如果

?設(shè)

否則，設(shè)置

結(jié)束?

吉布斯抽樣是Metropolis Hastings的一個(gè)特例。它涉及一個(gè)總是被接受的提議（總是有一個(gè)Metropolis-Hastings比率為1）。

我們應(yīng)用Metropolis Hastings算法來(lái)估計(jì)標(biāo)準(zhǔn)G-BLUP模型中回歸系數(shù)的方差成分。

對(duì)于G-BLUP模型。

其中

，

代表表型的向量和基因型的矩陣。?

是標(biāo)記效應(yīng)的向量，

是模型殘差的向量，殘差為正態(tài)分布，均值為0，方差為

和

。

考慮到其余參數(shù)，

的條件后驗(yàn)密度為：

這是一個(gè)逆卡方分布。

假設(shè)我們需要使

的先驗(yàn)盡可能地不具信息性。一種選擇是設(shè)置

，

，并使用拒絕抽樣來(lái)估計(jì)

；但是，設(shè)置S0=0可能會(huì)導(dǎo)致算法卡在0處。 因此，我們需要一個(gè)可以替代逆卡方分布的先驗(yàn)，并且可以非常靈活。為此，我們建議使用β分布。由于所得到的后驗(yàn)不是一個(gè)合適的分布，Metropolis Hastings算法將是獲得

后驗(yàn)樣本的一個(gè)好選擇。

這里我們把

作為我們的提議分布Q。因此。

我們的目標(biāo)分布是

的正態(tài)似然與

的β先驗(yàn)的乘積。由于β分布的域在0和1之間，我們用變量

來(lái)代替β先驗(yàn)，其中MAX是一個(gè)確保大于

的數(shù)字，這樣

。

其中α1和α2是β分布的形狀參數(shù)，其平均值由

給出。

我們按照上面的算法步驟，計(jì)算出我們的接受率，如下所示。

然后我們從均勻分布中抽取一個(gè)隨機(jī)數(shù)u，如果

，則接受樣本點(diǎn)

，否則我們拒絕該點(diǎn)并保留當(dāng)前值，再次迭代直至收斂。

Metropolis Hastings 算法

1. MetropolisHastings=function(p, ...)

2. 


3. chain[1]=x

4. for (i in 1:nIter) {

5. y[i] <-(SS+S0)/rchisq(df=DF,n=1,...)

6. 


7. logp.old[i]=-(p/2)*log(chai) - (SS/(2*chain) + (shape1-1)*(log(chain[i]/(MAX)))+(shape2-1)*(log(1-(chain[i]/(MAX))

8. 


9. logp.new[i]=-(p/2)*log(y[i]) - (SS/(2*y[i])) + (shape1-1)*(log(y[i]/(MAX)))+(shape2-1)*(log(1-(y[i]/(MAX))

10. chain[i+1] = ifelse (runif(1)<AP[i] , y[i], chain[i],...)

吉布斯采樣器?

1. gibbs=function(p,...)

2. b = rnorm(p,0,sqrt(varb),...)

3. for (i in 1:Iter) {

4. chain[i] <-(S+S0)/rchisq(df=DF,n=1,...)

繪制圖

1. plot = function(out1,out2)

2. plot(density(chain1),xlim=xlim)

3. lines(density(chain2),xlim=xlim)

4. abline(v=varb,col="red",lwd=3)

5. 


設(shè)置參數(shù)?

運(yùn)行吉布斯采樣器

1. ##################

2. out1=gibbs(p=sample.small,...)

3. out2=gibbs(p=sample.large,...)

在不同的情況下運(yùn)行METROPOLIS HASTINGS

小樣本量，先驗(yàn)

out.mh=mh(p=sample.small,nIter=nIter,varb=varb,shape1=shape.flat,shape2=shape.flat, MAX=MAX)

?

樣本量小，β值的形狀1參數(shù)大

1. p=sample.small

2. nIter

3. varb

4. shape.skew[1]

5. shape.skew[2]

6. MAX

plot(out.mh, out.gs_1)

?

樣本量小，β值的形狀1參數(shù)大

MetropolisHastings(p)

makeplot(out.mh, out.gs_1)

1. ## Summary of chain for MH:

2. ## ? ?Min. 1st Qu. ?Median ? ?Mean 3rd Qu. ? ?Max.

3. ## ?0.2097 ?0.2436 ?0.2524 ?0.2698 ?0.2978 ?0.4658

樣本量小，β的形狀參數(shù)相同（大）

plot(out.mh, out1)

大的樣本量，先驗(yàn)

plot(out.mh, out2)

大樣本量，形狀1參數(shù)的β

plot(out.mh, out2)

大樣本量，β值的大形狀2參數(shù)

plot(out.mh, out_2)

大樣本量，β的形狀參數(shù)相同（大）

plot(out.mh, out2)

參考文獻(xiàn)

1. Gelman, Andrew, et al. Bayesian data analysis. Vol. 2. London: Chapman & Hall/CRC, 2014.

最受歡迎的見(jiàn)解

1.使用R語(yǔ)言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行

2.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型

3.R語(yǔ)言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

4.R語(yǔ)言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

5.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

6.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析

7.R語(yǔ)言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)

8.R語(yǔ)言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸

9.R語(yǔ)言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長(zhǎng)