全文下載鏈接：http://tecdat.cn/?p=4612

最近我們被客戶要求撰寫關于Gibbs抽樣的研究報告，包括一些圖形和統(tǒng)計輸出。 貝葉斯分析的許多介紹都使用了相對簡單的教學實例（例如，根據(jù)伯努利數(shù)據(jù)給出成功概率的推理）。雖然這很好地介紹了貝葉斯原理，但是這些原則的擴展并不是直截了當?shù)?/p>

這篇文章將概述這些原理如何擴展到簡單的線性回歸。我將導出感興趣參數(shù)的后驗條件分布，給出用于實現(xiàn)Gibbs采樣器的R代碼，并提出所謂的網格點方法。

貝葉斯模型

假設我們觀察數(shù)據(jù)

對于

我們的模型是

有興趣的是作出推論

如果我們在方差項之前放置正態(tài)前向系數(shù)和反伽馬，那么這個數(shù)據(jù)的完整貝葉斯模型可以寫成：

假設超參數(shù)

是已知的，后面可以寫成一個常數(shù)的比例，

括號中的術語是數(shù)據(jù)或可能性的聯(lián)合分布。其他條款包括參數(shù)的聯(lián)合先驗分布（因為我們隱含地假設獨立前，聯(lián)合先驗因素）。

伴隨的R代碼的第0部分為該指定的“真實”參數(shù)從該模型生成數(shù)據(jù)。我們稍后將用這個數(shù)據(jù)估計一個貝葉斯回歸模型來檢查我們是否可以恢復這些真實的參數(shù)。

tphi<-rinvgamma(1,?shape=a,?rate=g) tb0<-rnorm(1,?m0,?sqrt(t0)?) tb1<-rnorm(1,?m1,?sqrt(t1)?) tphi;?tb0;?tb1; y<-rnorm(n,?tb0?+?tb1*x,?sqrt(tphi))

吉布斯采樣器

為了從這個后驗分布中得出，我們可以使用Gibbs抽樣算法。吉布斯采樣是一種迭代算法，從每個感興趣的參數(shù)的后驗分布產生樣本。它通過按照以下方式從每個參數(shù)的條件后面依次繪制：

可以看出，剩下的1,000個抽簽是從后驗分布中抽取的。這些樣本不是獨立的。繪制順序是隨機游走在后空間，空間中的每一步取決于前一個位置。通常還會使用間隔期（這里不做）。這個想法是，每一個平局可能依賴于以前的平局，但不能作為依賴于10日以前的平局。

點擊標題查閱往期內容

使用R語言進行Metroplis-in-Gibbs采樣和MCMC運行分析

左右滑動查看更多

01

02

03

04

條件后驗分布

要使用Gibbs，我們需要確定每個參數(shù)的條件后驗。

它有助于從完全非標準化的后驗開始：

為了找到參數(shù)的條件后驗，我們簡單地刪除不包含該參數(shù)的關節(jié)后驗的所有項。例如，常數(shù)項

條件后驗：

同樣的，

條件后驗可以被認為是另一個逆伽馬分布，有一些代數(shù)操作。

條件后驗

不那么容易識別。但是如果我們愿意使用網格方法，我們并不需要經過任何代數(shù)。

考慮網格方法。網格方法是非常暴力的方式（在我看來）從其條件后驗分布進行抽樣。這個條件分布只是一個函數(shù)

。所以我們可以評估一定的密度

值。在R表示法中，這可以是grid = seq（-10，10，by = .001）。這個序列是點的“網格”。

那么在每個網格點評估的條件后驗分布告訴我們這個抽取的相對可能性。

然后，我們可以使用R中的sample（）函數(shù)從這些網格點中抽取，抽樣概率與網格點處的密度評估成比例。

??for(i?in?1:length(p)?){ ????p[i]<-?(-(1/(2*phi))*sum(?(y?-?(grid[i]+b1*x))^2?))??+?(?-(1/(2*t0))*(grid[i]?-?m0)^2) ??} ?? ??draw<-sample(grid,?size?=?1,?prob?=?exp(1-p/max(p)))

這在R代碼的第一部分的函數(shù)rb0cond（）和rb1cond（）中實現(xiàn)。

使用網格方法時遇到數(shù)值問題是很常見的。由于我們正在評估網格中未標準化的后驗，因此結果可能會變得相當大或很小。這可能會在R中產生Inf和-Inf值。

例如，在函數(shù)rb0cond（）和rb1cond（）中，我實際上評估了派生的條件后驗分布的對數(shù)。然后，我通過從所有評估的最大值減去每個評估之前歸一化，然后從對數(shù)刻度取回。

我們不需要使用網格方法來從條件的后面繪制。

因為它來自已知的分布

請注意，這種網格方法有一些缺點。

首先，這在計算上是復雜的。通過代數(shù)，希望得到一個已知的后驗分布，從而在計算上更有效率。

其次，網格方法需要指定網格點的區(qū)域。如果條件后驗在我們指定的[-10,10]的網格間隔之外具有顯著的密度？在這種情況下，我們不會從條件后驗得到準確的樣本。記住這一點非常重要，并且需要廣泛的網格間隔進行實驗。所以，我們需要聰明地處理數(shù)字問題，例如在R中接近Inf和-Inf值的數(shù)字。

仿真結果

現(xiàn)在我們可以從每個參數(shù)的條件后驗進行采樣，我們可以實現(xiàn)Gibbs采樣器。這是在附帶的R代碼的第2部分中完成的。它編碼上面在R中概述的相同的算法。

iter<-1000burnin<-101phi<-b0<-b1<-numeric(iter) phi[1]<-b0[1]<-b1[1]<-6

結果很好。下圖顯示了1000個吉布斯（Gibbs）樣品的序列。紅線表示我們模擬數(shù)據(jù)的真實參數(shù)值。第四幅圖顯示了截距和斜率項的后面聯(lián)合，紅線表示輪廓。

z?<-?kde2d(b0,?b1,?n=50) plot(b0,b1,?pch=19,?cex=.4) contour(z,?drawlabels=FALSE,?nlevels=10,?col='red',?add=TRUE)

總結一下，我們首先推導了一個表達式，用于參數(shù)的聯(lián)合分布。然后我們概述了從后面抽取樣本的Gibbs算法。在這個過程中，我們認識到Gibbs方法依賴于每個參數(shù)的條件后驗分布的順序繪制。這是一個容易識別的已知的分布。對于斜率和截距項，我們決定用網格方法來規(guī)避代數(shù)。

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析》。

點擊標題查閱往期內容

python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉換隨機波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言中貝葉斯網絡（BN）、動態(tài)貝葉斯網絡、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估