拓端tecdat|R語(yǔ)言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
原文鏈接：http://tecdat.cn/?p=23524
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
在本文中，我想向你展示如何使用R的Metropolis采樣從貝葉斯Poisson回歸模型中采樣。
Metropolis-Hastings算法
Metropolis-Hastings抽樣算法是一類馬爾科夫鏈蒙特卡洛（MCMC）方法，其主要思想是生成一個(gè)馬爾科夫鏈
使其平穩(wěn)分布為目標(biāo)分布。這種算法最常見(jiàn)的應(yīng)用之一是在貝葉斯統(tǒng)計(jì)中從后驗(yàn)密度中取樣，這也是本文的目標(biāo)。
該算法規(guī)定對(duì)于一個(gè)給定的狀態(tài)Xt，如何生成下一個(gè)狀態(tài)?
?有一個(gè)候選點(diǎn)Y，它是從一個(gè)提議分布?
,中生成的，根據(jù)決策標(biāo)準(zhǔn)被接受，所以鏈條在時(shí)間t+1時(shí)移動(dòng)到狀態(tài)Y，即Xt+1=Y或被拒絕，所以鏈條在時(shí)間t+1時(shí)保持在狀態(tài)Xt，即Xt+1=Xt。
Metropolis 采樣
在Metropolis算法中，提議分布是對(duì)稱的，也就是說(shuō)，提議分布??
滿足
?
，所以Metropolis采樣器產(chǎn)生馬爾科夫鏈的過(guò)程如下。
選擇一個(gè)提議分布
. 在選擇它之前，了解這個(gè)函數(shù)中的理想特征。
從提議分布g中生成X0。
重復(fù)進(jìn)行，直到鏈?zhǔn)諗康揭粋€(gè)平穩(wěn)的分布。
從?
生成Y.
從Uniform(0, 1)中生成U。
如果?
, 接受Y并設(shè)置Xt+1=Y，否則設(shè)置Xt+1=Xt。這意味著候選點(diǎn)Y被大概率地接受
.
遞增t.
貝葉斯方法
正如我之前提到的，我們要從定義為泊松回歸模型的貝葉斯中取樣。
對(duì)于貝葉斯分析中的參數(shù)估計(jì)，我們需要找到感興趣的模型的似然函數(shù)，在這種情況下，從泊松回歸模型中找到。
現(xiàn)在我們必須為每個(gè)參數(shù)β0和β1指定一個(gè)先驗(yàn)分布。我們將對(duì)這兩個(gè)參數(shù)使用無(wú)信息的正態(tài)分布，β0～N(0,100)和β1～N(0,100) 。
最后，我們將后驗(yàn)分布定義為先驗(yàn)分布和似然分布的乘積。
使用Metropolis采樣器時(shí)，后驗(yàn)分布將是目標(biāo)分布。
計(jì)算方法
這里你將學(xué)習(xí)如何使用R語(yǔ)言的Metropolis采樣器從參數(shù)β0和β1的后驗(yàn)分布中采樣。
數(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ù)。在這種情況下，我們將使用這個(gè)函數(shù)的對(duì)數(shù)，這是強(qiáng)烈建議的，以避免在運(yùn)行算法時(shí)出現(xiàn)數(shù)字問(wèn)題。
LikelihoodFunction <- function(param){beta0 <- param[1]beta1 <- param[2]lambda <- exp(beta1*X + beta0)# 對(duì)數(shù)似然函數(shù)loglikelihoods <- sum(dpois(y, lambda = lambda, log=T))return(loglikelihoods)}
先驗(yàn)分布
接下來(lái)我們定義參數(shù)β0和β1的先驗(yàn)分布。與似然函數(shù)一樣，我們將使用先驗(yàn)分布的對(duì)數(shù)。
beta0prior <- dnorm(beta0, 0, sqrt(100), log=TRUE)beta1prior <- dnorm(beta1, 0, sqrt(100), log=TRUE)return(beta0prior + beta1prior) #先驗(yàn)分布的對(duì)數(shù)
后驗(yàn)分布
由于我們是用對(duì)數(shù)工作的，我們把后驗(yàn)分布定義為似然函數(shù)的對(duì)數(shù)與先驗(yàn)分布的對(duì)數(shù)之和。記住，這個(gè)函數(shù)是我們的目標(biāo)函數(shù)f(.)，我們要從中取樣。
提議函數(shù)
最后，我們定義提議分布g(.|Xt)。由于我們將使用Metropolis采樣器，提議分布必須是對(duì)稱的，并且取決于鏈的當(dāng)前狀態(tài)，因此我們將使用正態(tài)分布，其平均值等于當(dāng)前狀態(tài)下的參數(shù)值。
Metropolis 采樣器
最后，我們編寫代碼，幫助我們執(zhí)行Metropolis采樣器。在這種情況下，由于我們使用的是對(duì)數(shù)，我們必須將候選點(diǎn)Y被接受的概率定義為。
# 創(chuàng)建一個(gè)數(shù)組來(lái)保存鏈的值chain[1, ] <- startvalue # 定義鏈的起始值for (i in 1:iterations){# 從提議函數(shù)生成YY <- ProposalFunction(chain[i, ])# 候選點(diǎn)被接受的概率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)無(wú)法代表真實(shí)的基礎(chǔ)后驗(yàn)分布。那么，為了減少自相關(guān)，我們可以只使用鏈上的每一個(gè)n個(gè)值來(lái)稀釋樣本。在這種情況下，我們將在算法的每20次迭代中為我們的最終鏈選擇一個(gè)值。
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個(gè)值burnIn <- 5000
在這里，你可以看到ACF圖，它給我們提供了任何序列與其滯后值的自相關(guān)值。在這種情況下，我們展示了初始MCMC鏈的ACF圖和對(duì)兩個(gè)參數(shù)的樣本進(jìn)行稀釋后的最終鏈。從圖中我們可以得出結(jié)論，所使用的程序?qū)嶋H上能夠大大減少自相關(guān)。
結(jié)果
在這一節(jié)中，我們介紹了由Metropolis采樣器產(chǎn)生的鏈以及它對(duì)參數(shù)β0和β1的分布。參數(shù)的真實(shí)值由紅線表示。
與glm()的比較
現(xiàn)在我們必須將使用Metropolis采樣得到的結(jié)果與glm()函數(shù)進(jìn)行比較，glm()函數(shù)用于擬合廣義linera模型。
下表列出了參數(shù)的實(shí)際值和使用Metropolis采樣器得到的估計(jì)值的平均值。
## ? ? ? True value Mean MCMC ? ? ? glm## beta0 ?1.0578047 1.0769213 1.0769789## beta1 ?0.8113144 0.8007347 0.8009269
結(jié)論
從結(jié)果來(lái)看，我們可以得出結(jié)論，使用Metropolis采樣器和glm()函數(shù)得到的泊松回歸模型的參數(shù)β0和β1的估計(jì)值非常相似，并且接近于參數(shù)的實(shí)際值。另外，必須認(rèn)識(shí)到先驗(yàn)分布、建議分布和鏈的初始值的選擇對(duì)結(jié)果有很大的影響，因此這種選擇必須正確進(jìn)行。
最受歡迎的見(jiàn)解
1.R語(yǔ)言多元Logistic邏輯回歸 應(yīng)用案例
2.面板平滑轉(zhuǎn)移回歸(PSTR)分析案例實(shí)現(xiàn)
3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）
4.R語(yǔ)言泊松Poisson回歸模型分析案例
5.R語(yǔ)言回歸中的Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)
6.r語(yǔ)言中對(duì)LASSO回歸，Ridge嶺回歸和Elastic Net模型實(shí)現(xiàn)
7.在R語(yǔ)言中實(shí)現(xiàn)Logistic邏輯回歸
8.python用線性回歸預(yù)測(cè)股票價(jià)格
9.R語(yǔ)言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)
