R語言Gibbs抽樣的貝葉斯簡單線性迴歸仿真分析

原文連接：http://tecdat.cn/?p=4612

貝葉斯分析的許多介紹都使用了相對簡單的教學實例（例如根據伯努利數據給出成功機率的推理）。雖然這能夠很好地介紹貝葉斯原理，可是將這些原理擴展到迴歸並非直接的。

這篇文章將概述這些原理如何擴展到簡單的線性迴歸。我將導出感興趣參數的後驗條件分佈，給出用於實現Gibbs採樣器的R代碼，並提出所謂的網格方法。

貝葉斯模型

假設咱們觀察數據​

對於咱們的模型​是

​

有趣的是推斷​和​。若是咱們將正態先驗放在係數上，將反伽瑪先驗放在方差項上，則此數據的完整貝葉斯模型能夠寫爲：

​

假設超參數​是已知的，後驗能夠被寫入到一個比例常數，

​

括號中的項是是數據的聯合分佈或機率。其餘項包括參數的聯合先驗分佈。

​

R代碼從指定的「真實」參數模型生成數據。咱們稍後將用這個數據估計貝葉斯迴歸模型來檢查是否能夠恢復這些真實的參數。

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採樣算法。吉布斯採樣是一種迭代算法，從每一個感興趣的參數的後驗分佈生成樣本。它經過如下方式從每一個參數的條件後驗分佈依次得出的：

​

能夠看出，剩下的1,000個樣本是從後驗分佈中抽取的。這些樣本不是獨立的。每一個步驟都取決於先前的位置。

條件後驗分佈

要使用Gibbs，咱們須要肯定每一個參數的條件後驗。

​

爲了找到參數的條件後驗，咱們簡單地刪除不包含參數後驗的全部項。例如，常數項條件後驗：

​

類似地，

​

條件後驗能夠被認爲是另外一個逆伽馬分佈。

​

條件後驗​不那麼容易識別。可是若是使用網格方法，咱們不須要經過代數方法。

​考慮網格方法。網格方法是暴力方法，從其條件後驗分佈進行抽樣。這個條件分佈只是一個函數​。因此咱們能夠評估密度​值。在R中，能夠是grid = seq（-10，10，by = .001）。這個序列是點的「網格」。

那麼在每一個網格點評估的條件後驗分佈告訴咱們這個抽樣的相對機率。

而後，咱們能夠使用R中的sample（）函數從這些網格中抽取，抽樣機率與網格處的密度評估成比例。

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代碼的第一部分的函數rb0cond（）和rb1cond（）中實現。

使用網格方法時遇到數值問題是很常見的。因爲咱們在評估網格中未標準化的後驗，所以結果可能會變得至關大或很小。可能會在R中產生Inf和-Inf值。

例如，在函數rb0cond（）和rb1cond（）中，我實際上評估了條件後驗分佈的對數。而後，我進行歸一化和對數化。

咱們不須要使用網格方法來從後驗條件抽樣，由於它來自已知的分佈​。

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

首先，這在計算上是複雜的。經過代數，但願獲得一個已知的後驗分佈，從而在計算上更有效率。

其次，網格方法須要指定網格點的區域。若是條件後驗在咱們指定的[-10,10]的網格區間以外具備顯着的密度，在這種狀況下，咱們不會從條件後驗獲得準確的樣本。而且須要普遍的網格區間進行實驗。因此，咱們須要靈活地處理數字問題，例如在R中接近Inf和-Inf值的數字。

仿真結果

如今咱們能夠從每一個參數的條件後驗進行採樣，咱們能夠實現Gibbs採樣器。

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

結果很好。下圖顯示了1000個吉布斯（Gibbs）樣本的序列。紅線表示咱們模擬數據的真實參數值。第四幅圖顯示了截距和斜率項的後驗聯合分佈，紅線表示等高線。

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

​

總結一下，咱們首先推導了一個表達式，用於參數的聯合分佈。而後咱們概述了從後驗抽取樣本的Gibbs算法。在這個過程當中，咱們認識到Gibbs方法依賴於每一個參數的條件後驗分佈，這是容易識別的已知的分佈。對於斜率和截距項，咱們決定用網格方法來規避代數方法。這樣作的好處是咱們能夠繞開不少代數運算。代價是增長了計算複雜性。

參考文獻

1.matlab使用貝葉斯優化的深度學習

2.matlab貝葉斯隱馬爾可夫hmm模型實現

3.R語言Gibbs抽樣的貝葉斯簡單線性迴歸仿真

4.R語言中的block Gibbs吉布斯採樣貝葉斯多元線性迴歸

5.R語言中的Stan機率編程MCMC採樣的貝葉斯模型

6.Python用PyMC3實現貝葉斯線性迴歸模型

7.R語言使用貝葉斯 層次模型進行空間數據分析

8.R語言隨機搜索變量選擇SSVS估計貝葉斯向量自迴歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實現