R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
原文鏈接：http://tecdat.cn/?p=19889
?
如果您可以寫出模型的似然函數，則?Metropolis-Hastings算法可以負責其余部分（即MCMC ）。我寫了r代碼來簡化對任意模型的后驗分布的估計。具體如下：
1）定義模型（即概率先驗）。在此示例中，讓我們構建一個簡單的線性回歸模型（對數）。
a<-pars[1]????? #截距b<-pars[2]????? #斜率sd_e<-pars[3]?? #殘差if(sd_e<=0){return(NaN)}log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
先驗：
epsilon<-pars[3]??? #殘差prior_a<-dnorm(a,0,100,log=TRUE)???? ##所有的非信息性先驗prior_b<-dnorm(b,0,100,log=TRUE)???? ## 參數.prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)
現在讓我們模擬一些數據以進行運行測試：
x<-runif(30,5,15)y<-x+rnorm(30,0,5) ##斜率=1, 截距=0, epsilon=5
2）Metro Hastings 完成所有工作。
MH(li_func=li_reg,pars=c(0,1,1),
3）您可以使用plotMH（）查看所有模型參數的后驗
plot(mcmc)
繪制所有參數之間的相關性。
4）輸出后驗置信區(qū)間。
BCI#? ??????????? 0.025??? 0.975# a?????? -5.3345970 6.841016# b??????? 0.4216079 1.690075# epsilon? 3.8863393 6.660037
接下來，我想提供一種直觀的方法來可視化此算法運行的情況。
主要思想是從分布中抽取樣本。積分很重要，貝葉斯定理本身：
P（θ| D）= P（D |θ）P（θ）/ P（D）
其中P（D）是觀察數據的無條件概率。由于這不依賴于推斷的模型（θ）參數，因此P（D）是歸一化常數。
因此，我們有一個非歸一化的概率密度函數，我們希望通過隨機抽樣來估計。對于復雜的模型而言，隨機抽樣本身的過程通常很困難，因此，我們使用馬爾可夫鏈來探索分布。我們需要一個鏈，如果運行時間足夠長，它將作為目標分布的隨機樣本整體。我們構建的馬爾可夫鏈的這種特性稱為?遍歷性。Metropolis-Hastings算法是構建這種鏈的一種方法。
步驟：
在參數空間k_X中選擇一些起點
選擇一個候選點k_Y?N（k_X，σ）。這通常稱為提議分布。
移至候選點的概率為：min（π（k_Y）/π（K_X），1）
重復。
以下代碼通過簡單的正態(tài)目標分布演示了此過程。
###???? Metropolis-Hastings 可視化??? ??? ??? ??? #######k_X = seed; ##將k_X設置為種子位置for(i in 1:iter){track<-c(track,k_X)??? ## 鏈k_Y = rnorm(1,k_X,prop_sd) ##候選點## -- 繪制鏈的核密度估計lines(density(track,adjust=1.5),col='red',lwd=2)## -- 繪制鏈plot(track,1:i,xlim=plot_range,main='',type='l',ylab='Trace')## -- 繪制目標分布和提議分布curve(dnorm(x,k_X,prop_sd),col='black',add=TRUE)abline(v=k_X,lwd=2)## 接受概率為a_X_Yif (log(runif(1))<=a_X_Y)points(k_Y,0,pch=19,col='green',cex=2)## 調整提議if(i>100)prop_sd=sd(track[floor(i/2):i])
?
該算法實現中的一個普遍問題是σ的選擇。當σ接近目標分布的標準偏差時，將發(fā)生有效混合（鏈收斂到目標分布）。當我們不知道這個值時。我們可以允許σ根據到目前為止的鏈歷史記錄進行調整。在上面的示例中，將σ更新為鏈中某些先驗點的標準偏差值。
輸出為多頁pdf，可以滾動瀏覽。
頂部顯示了目標分布（藍色虛線）和通過MCMC樣本對目標進行的核平滑估計。第二面板顯示了鏈的軌跡，底部顯示了算法本身的步驟。
注意：請注意，前100次左右的迭代是目標分布的較差表示。在實踐中，我們將“預燒”該鏈的前n個迭代-通常是前100-1000個。
?
?
?
最受歡迎的見解
1.使用R語言進行METROPLIS-IN-GIBBS采樣和MCMC運行
2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
3.R語言實現MCMC中的Metropolis–Hastings算法與吉布斯采樣
4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
7.R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數
8.R語言使用Metropolis- Hasting抽樣算法進行邏輯回歸
9.R語言中基于混合數據抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長
?
