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

最近我們被客戶要求撰寫關(guān)于貝葉斯簡單線性回歸的研究報告，包括一些圖形和統(tǒng)計輸出。

貝葉斯分析的許多介紹都使用了相對簡單的教學(xué)實例（例如，根據(jù)伯努利數(shù)據(jù)給出成功概率的推理）。雖然這很好地介紹了貝葉斯原理，但是這些原則的擴展并不是直截了當(dāng)?shù)?/h1>

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

貝葉斯模型

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

對于

我們的模型是

有興趣的是作出推論

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

假設(shè)超參數(shù)

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

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

伴隨的R代碼的第0部分為該指定的“真實”參數(shù)從該模型生成數(shù)據(jù)。我們稍后將用這個數(shù)據(jù)估計一個貝葉斯回歸模型來檢查我們是否可以恢復(fù)這些真實的參數(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ù)的后驗分布產(chǎn)生樣本。它通過按照以下方式從每個參數(shù)的條件后面依次繪制：

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

點擊標(biāo)題查閱往期內(nèi)容

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

左右滑動查看更多

01

02

03

04

條件后驗分布

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

它有助于從完全非標(biāo)準化的后驗開始：

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

條件后驗：

同樣的，

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

條件后驗

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

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

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

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

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

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

??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)。

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

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

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

因為它來自已知的分布

請注意，這種網(wǎng)格方法有一些缺點。

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

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

仿真結(jié)果

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

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

結(jié)果很好。下圖顯示了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)

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

點擊文末?“閱讀原文”

獲取全文完整資料。

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

點擊標(biāo)題查閱往期內(nèi)容

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