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

最近我們被客戶要求撰寫關(guān)于采樣的研究報告，包括一些圖形和統(tǒng)計輸出。

在本文中，我想向你展示如何使用R的Metropolis采樣從貝葉斯Poisson回歸模型中采樣。

Metropolis-Hastings算法

Metropolis-Hastings抽樣算法是一類馬爾科夫鏈蒙特卡洛（MCMC）方法，其主要思想是生成一個馬爾科夫鏈

使其平穩(wěn)分布為目標(biāo)分布。這種算法最常見的應(yīng)用之一是在貝葉斯統(tǒng)計中從后驗密度中取樣，這也是本文的目標(biāo)。

該算法規(guī)定對于一個給定的狀態(tài)Xt，如何生成下一個狀態(tài)?

?有一個候選點Y，它是從一個提議分布?

,中生成的，根據(jù)決策標(biāo)準(zhǔn)被接受，所以鏈條在時間t+1時移動到狀態(tài)Y，即Xt+1=Y或被拒絕，所以鏈條在時間t+1時保持在狀態(tài)Xt，即Xt+1=Xt。

Metropolis 采樣

在Metropolis算法中，提議分布是對稱的，也就是說，提議分布??

滿足

?

，所以Metropolis采樣器產(chǎn)生馬爾科夫鏈的過程如下。

1. 選擇一個提議分布

1. . 在選擇它之前，了解這個函數(shù)中的理想特征。

2. 從提議分布g中生成X0。

3. 重復(fù)進(jìn)行，直到鏈?zhǔn)諗康揭粋€平穩(wěn)的分布。

• 從?

• 生成Y.

• 從Uniform(0, 1)中生成U。

• 如果?

• , 接受Y并設(shè)置Xt+1=Y，否則設(shè)置Xt+1=Xt。這意味著候選點Y被大概率地接受

• .

• 遞增t.

貝葉斯方法

正如我之前提到的，我們要從定義為泊松回歸模型的貝葉斯中取樣。

對于貝葉斯分析中的參數(shù)估計，我們需要找到感興趣的模型的似然函數(shù)，在這種情況下，從泊松回歸模型中找到。

現(xiàn)在我們必須為每個參數(shù)β0和β1指定一個先驗分布。我們將對這兩個參數(shù)使用無信息的正態(tài)分布，β0～N(0,100)和β1～N(0,100) 。

最后，我們將后驗分布定義為先驗分布和似然分布的乘積。

使用Metropolis采樣器時，后驗分布將是目標(biāo)分布。

計算方法

這里你將學(xué)習(xí)如何使用R語言的Metropolis采樣器從參數(shù)β0和β1的后驗分布中采樣。

數(shù)據(jù)

首先，我們從上面介紹的泊松回歸模型生成數(shù)據(jù)。

n?<-?1000?#??樣本大小J?<-?2?#?參數(shù)的數(shù)量X?<-?runif(n,-2,2)?#?生成自變量的值beta?<-?runif(J,-2,2)?#生成參數(shù)的值y?<-?rpois(n,?lambda?=?lambda)?#?生成因變量的值

似然函數(shù)

現(xiàn)在我們定義似然函數(shù)。在這種情況下，我們將使用這個函數(shù)的對數(shù)，這是強(qiáng)烈建議的，以避免在運行算法時出現(xiàn)數(shù)字問題。

LikelihoodFunction?<-?function(param){????????beta0?<-?param[1]? ????????beta1?<-?param[2]? ????????lambda?<-?exp(beta1*X?+?beta0)????????#?對數(shù)似然函數(shù)????????loglikelihoods?<-?sum(dpois(y,?lambda?=?lambda,?log=T))? ????????return(loglikelihoods)}

先驗分布

接下來我們定義參數(shù)β0和β1的先驗分布。與似然函數(shù)一樣，我們將使用先驗分布的對數(shù)。

????????beta0prior?<-?dnorm(beta0,?0,?sqrt(100),?log=TRUE)????????beta1prior?<-?dnorm(beta1,?0,?sqrt(100),?log=TRUE)????????return(beta0prior?+?beta1prior)?#先驗分布的對數(shù)

后驗分布

由于我們是用對數(shù)工作的，我們把后驗分布定義為似然函數(shù)的對數(shù)與先驗分布的對數(shù)之和。記住，這個函數(shù)是我們的目標(biāo)函數(shù)f(.)，我們要從中取樣。

提議函數(shù)

最后，我們定義提議分布g(.|Xt)。由于我們將使用Metropolis采樣器，提議分布必須是對稱的，并且取決于鏈的當(dāng)前狀態(tài)，因此我們將使用正態(tài)分布，其平均值等于當(dāng)前狀態(tài)下的參數(shù)值。

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

R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

左右滑動查看更多

01

02

03

04

Metropolis 采樣器

最后，我們編寫代碼，幫助我們執(zhí)行Metropolis采樣器。在這種情況下，由于我們使用的是對數(shù)，我們必須將候選點Y被接受的概率定義為。

????????#?創(chuàng)建一個數(shù)組來保存鏈的值????????chain[1,?]?<-?startvalue?#?定義鏈的起始值????????for?(i?in?1:iterations){ ????????????????#?從提議函數(shù)生成Y????????????????Y?<-?ProposalFunction(chain[i,?])? ????????????????#?候選點被接受的概率???????????????????????????????????????????PosteriorFunction(chain[i,?])) ????????????????#?接受或拒絕Y的決策標(biāo)準(zhǔn)?????????????????if?(runif(1)?<?probability)?{ ????????????????????????chain[i+1,?]?<-?Y ????????????????}else{? ????????????????????????chain[i+1,?]?<-?chain[i,?]

由于MCMC鏈具有很強(qiáng)的自相關(guān)，它可能產(chǎn)生的樣本在短期內(nèi)無法代表真實的基礎(chǔ)后驗分布。那么，為了減少自相關(guān)，我們可以只使用鏈上的每一個n個值來稀釋樣本。在這種情況下，我們將在算法的每20次迭代中為我們的最終鏈選擇一個值。

startvalue?<-?c(0,?0)?#?定義鏈條的起始值#每20次迭代選擇最終鏈的值for?(i?in?1:10000){ ????????if?(i?==?1){ ????????????????cfinal[i,?]?<-?chain[i*20,]????????}?else?{ ????????????????cfinal[i,?]?<-?chain[i*20,]#?刪除鏈上的前5000個值burnIn?<-?5000

在這里，你可以看到ACF圖，它給我們提供了任何序列與其滯后值的自相關(guān)值。在這種情況下，我們展示了初始MCMC鏈的ACF圖和對兩個參數(shù)的樣本進(jìn)行稀釋后的最終鏈。從圖中我們可以得出結(jié)論，所使用的程序?qū)嶋H上能夠大大減少自相關(guān)。

結(jié)果

在這一節(jié)中，我們介紹了由Metropolis采樣器產(chǎn)生的鏈以及它對參數(shù)β0和β1的分布。參數(shù)的真實值由紅線表示。

與glm()的比較

現(xiàn)在我們必須將使用Metropolis采樣得到的結(jié)果與glm()函數(shù)進(jìn)行比較，glm()函數(shù)用于擬合廣義linera模型。

下表列出了參數(shù)的實際值和使用Metropolis采樣器得到的估計值的平均值。

##???????True?value?Mean?MCMC???????glm##?beta0??1.0578047?1.0769213?1.0769789##?beta1??0.8113144?0.8007347?0.8009269

結(jié)論

從結(jié)果來看，我們可以得出結(jié)論，使用Metropolis采樣器和glm()函數(shù)得到的泊松回歸模型的參數(shù)β0和β1的估計值非常相似，并且接近于參數(shù)的實際值。另外，必須認(rèn)識到先驗分布、建議分布和鏈的初始值的選擇對結(jié)果有很大的影響，因此這種選擇必須正確進(jìn)行。

本文摘選?《?R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型?》?，點擊“閱讀原文”獲取全文完整資料。

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

Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(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的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
R語言隨機(jī)搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機(jī)波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計