使用R語言進(jìn)行Metroplis-in-Gibbs采樣和MCMC運(yùn)行分析
原文鏈接：http://tecdat.cn/?p=12200
對(duì)于許多模型，例如邏輯模型，沒有共軛先驗(yàn)分布。因此，吉布斯采樣不適用。
這篇文章展示了我們?nèi)绾问褂肕etropolis-Hastings（MH）從每次Gibbs迭代中的非共軛條件后驗(yàn)對(duì)象中進(jìn)行采樣–比網(wǎng)格方法更好的替代方法。
我將說明該算法，給出一些R代碼結(jié)果，然后分析R代碼以識(shí)別MH算法中的瓶頸。
模型
此示例的模擬數(shù)據(jù)是包含
患者的橫截面數(shù)據(jù)集。有一個(gè)二元因變量Y，一個(gè)二元處理變量A，一個(gè)因子變量age。年齡是具有3個(gè)等級(jí)的分類變量。我用貝葉斯邏輯回歸建模：
對(duì)于Metroplis-in-Gibbs吉布斯采樣來說，這是一個(gè)相當(dāng)不錯(cuò)的示例：
我們有一個(gè)二進(jìn)制結(jié)果，為此我們采用了非線性鏈接函數(shù)。
我們有一個(gè)需要調(diào)整的因素。
我們正在估計(jì)我們關(guān)心的更多參數(shù)，但肯定會(huì)給采樣器帶來壓力。
非規(guī)范條件后驗(yàn)
讓我們看一下該模型的（非標(biāo)準(zhǔn)化）條件后驗(yàn)。
此條件分布不是已知分布，因此我們不能簡單地使用Gibbs從中進(jìn)行采樣。相反，在每個(gè)gibbs迭代中，我們需要另一個(gè)采樣步驟來從該條件后驗(yàn)中提取。第二個(gè)采樣器將是MH采樣器。
?
Metroplis-in-Gibbs采樣
目標(biāo)是從中取樣
。
MH采樣器的工作方式如下：
開始采樣。
讓我們假設(shè)將提議分布的方差設(shè)置為某個(gè)常數(shù)。
我們計(jì)算在上一次繪制時(shí)評(píng)估的非標(biāo)準(zhǔn)化密度與當(dāng)前提議的比率：? ?
如果該比率大于1，則當(dāng)前提議分布的密度高于先前值的密度。因此，我們“接受”了提議并確定了
。然后，我們使用以提議為中心的提議分布重復(fù)步驟2-4??，然后生成新提議。如果該比率小于1，則當(dāng)前建議值的密度低于先前建議。
因此，總是接受產(chǎn)生更高條件的后驗(yàn)評(píng)估的提議。但是，有時(shí)僅接受具有較低密度評(píng)估的提議-提議的相對(duì)密度評(píng)估越低，其接受的可能性就越低。
經(jīng)過多次迭代，從后驗(yàn)的高密度區(qū)域開始的抽樣被接受，并且被接受的序列“爬升”到高密度區(qū)域。一旦序列到達(dá)此高密度區(qū)域，它將趨于保持在那里。因此，這也類似于模擬退火。
這種表示法很容易擴(kuò)展到我們的4維示例：提議分布現(xiàn)在是4維多元高斯模型。代替標(biāo)量方差參數(shù)，我們有一個(gè)協(xié)方差矩陣。因此，我們的建議是系數(shù)的向量。從這個(gè)意義上講，我們運(yùn)行的是Gibbs –使用MH每次迭代繪制整個(gè)系數(shù)。
?
跳躍分布的方差是重要的參數(shù)。如果方差太小，則當(dāng)前提議可能會(huì)非常接近最后一個(gè)值，因此
也很可能接近1。因此，我們會(huì)非常頻繁地接受，但由于接受的值彼此之間非常接近，因此我們會(huì)攀升至較高在許多次迭代中慢慢降低密度區(qū)域。如果方差太大，則序列到達(dá)高密度區(qū)域后可能無法保留在該區(qū)域。
許多“自適應(yīng)” MH方法是此處描述的基本算法的變體，但包括調(diào)整周期以找到產(chǎn)生最佳接受率的跳躍分布方差。
MH中計(jì)算量最大的部分是密度評(píng)估。對(duì)于每個(gè)Gibbs迭代，我們必須兩次評(píng)估4維密度。
盡管很容易擴(kuò)展到高維度，但性能本身在高維度上會(huì)變差。
結(jié)果
這是我們感興趣的4個(gè)參數(shù)的MCMC鏈。紅線表示真實(shí)值。
for(i in 2:gibbs_iter){# 來自 phi 后驗(yàn)分布的樣本gibbs_res[i,p+1] <- rcond_post_phi(gibbs_res[i-1,1:p],alpha, gamma, lambda, p)# ?來自beta后驗(yàn)分布的樣本 （ 使用 Mh ）mh_draw <- rcond_post_beta_mh(gibbs_res[i-1,1:p], gibbs_res[i,p+1],lambda, X, Y, mh_trials=5, jump_v=.01)}par(mfrow=c(2,2))plot(gibbs_res[,1],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),main='Intercept')abline(h=-1,col='red')plot(gibbs_res[,2],type='l',xlab='MCMC Iterations',ylab=c('Coefficient Draw'),main='Age1')# 計(jì)算后驗(yàn)分布和置信區(qū)間post_burn_trim<-gibbs_res[seq(1000,gibbs_iter,100),]colMeans(post_burn_trim)apply(post_burn_trim, 2, quantile, p=c(.025,.975))
?
?
有一些改進(jìn)的空間：
接受率只有18％，我本可以調(diào)整跳躍分布協(xié)方差矩陣來獲得更好的接受率。
我認(rèn)為更多的迭代肯定會(huì)在這里有所幫助。這些鏈看起來不錯(cuò)，但仍然是自相關(guān)的。
關(guān)于貝葉斯范式的好處是，所有推斷都是使用后驗(yàn)分布完成的?，F(xiàn)在，系數(shù)估計(jì)值是對(duì)數(shù)化，但是如果我們需要比值，則只需對(duì)后驗(yàn)取冪。如果我們想要對(duì)比值進(jìn)行區(qū)間估計(jì)，那么我們就可以獲取指數(shù)后驗(yàn)的2.5％和97.5％。
下面是使用R分析，顯示了這一點(diǎn)。for循環(huán)運(yùn)行Gibbs迭代。在每個(gè)Gibbs迭代中，我都調(diào)用函數(shù)使用MH從參數(shù)向量的條件后驗(yàn)中得出圖形。
?
我們看到子例程log_cond（）是MH運(yùn)行中的瓶頸。此函數(shù)是beta的對(duì)數(shù)條件后驗(yàn)密度。
最受歡迎的見解
1.使用R語言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行
2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
3.R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
7.R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
8.R語言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸
9.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長
?
