詳解R語言MCMC:Metropolis-Hastings采樣用於回歸的貝葉斯估計
MCMC是從復雜概率模型中采樣的通用技術。
- 蒙特卡洛
- 馬爾可夫鏈
- Metropolis-Hastings算法
問題
如果需要計算有復雜後驗pdf p（θ| y）的隨機變量θ的函數f（θ）的平均值或期望值。
您可能需要計算後驗概率分佈p（θ）的最大值。
解決期望值的一種方法是從p（θ）繪制N個隨機樣本，當N足夠大時，我們可以通過以下公式逼近期望值或最大值
將相同的策略應用於通過從p（θ| y）采樣並取樣本集中的最大值來找到argmaxp（θ| y）。
解決方法
1.1直接模擬
1.2逆CDF
1.3拒絕/接受抽樣
如果我們不知道精確/標準化的pdf或非常復雜，則MCMC會派上用場。
馬爾可夫鏈
為瞭模擬馬爾可夫鏈，我們必須制定一個 過渡核T（xi，xj）。過渡核是從狀態xi遷移到狀態xj的概率。
馬爾可夫鏈的收斂性意味著它具有平穩分佈π。馬爾可夫鏈的統計分佈是平穩的,那麼它意味著分佈不會隨著時間的推移而改變。
Metropolis算法
對於一個Markov鏈是平穩的。基本上表示
處於狀態x並轉換為狀態x’的概率必須等於處於狀態x’並轉換為狀態x的概率
或者
方法是將轉換分為兩個子步驟；候選和接受拒絕。
令q（x’| x）表示 候選密度，我們可以使用概率 α（x’| x）來調整q 。
候選分佈 Q（X’| X）是給定的候選X的狀態X’的條件概率，
和 接受分佈 α（x’| x）的條件概率接受候選的狀態X’-X’。我們設計瞭接受概率函數，以滿足詳細的平衡。
該 轉移概率 可以寫成：
插入上一個方程式，我們有
Metropolis-Hastings算法
A的選擇遵循以下邏輯。
在q下從x到x’的轉移太頻繁瞭。因此，我們應該選擇α（x | x’）=1。但是，為瞭滿足 細致平穩，我們有
下一步是選擇滿足上述條件的接受。Metropolis-Hastings是一種常見的 選擇：
即，當接受度大於1時，我們總是接受，而當接受度小於1時，我們將相應地拒絕。因此，Metropolis-Hastings算法包含以下內容：
初始化：隨機選擇一個初始狀態x；
根據q（x’| x）隨機選擇一個新狀態x’;
3.接受根據α（x’| x）的狀態。如果不接受，則不會進行轉移，因此無需更新任何內容。否則，轉移為x’；
4.轉移到2，直到生成T狀態；
5.保存狀態x，執行2。
原則上，我們從分佈P（x）提取保存的狀態，因為步驟4保證它們是不相關的。必須根據候選分佈等不同因素來選擇T的值。 重要的是，尚不清楚應該使用哪種分佈q（x’| x）；必須針對當前的特定問題進行調整。
屬性
Metropolis-Hastings算法的一個有趣特性是它 僅取決於比率
是候選樣本x’與先前樣本xt之間的概率，
是兩個方向（從xt到x’，反之亦然）的候選密度之比。如果候選密度對稱，則等於1。
馬爾可夫鏈從任意初始值x0開始，並且算法運行多次迭代，直到“初始狀態”被“忘記”為止。這些被丟棄的樣本稱為預燒（burn-in）。其餘的x可接受值集代表分佈P（x）中的樣本
Metropolis采樣
一個簡單的Metropolis-Hastings采樣
讓我們看看從 伽瑪分佈 模擬任意形狀和比例參數，使用具有Metropolis-Hastings采樣算法。
下面給出瞭Metropolis-Hastings采樣器的函數。該鏈初始化為零，並在每個階段都建議使用N（a / b，a /（b * b））個候選對象。
基於正態分佈且均值和方差相同gamma的Metropolis-Hastings獨立采樣
從某種狀態開始xt。代碼中的x。在代碼中提出一個新的狀態x’候選計算“接受概率”
從[0,1] 得出一些均勻分佈的隨機數u；如果u <α接受該點，則設置xt + 1 = x’。否則，拒絕它並設置xt + 1 = xt。
MH可視化
set.seed(123)
for (i in 2:n) {
can <- rnorm(1, mu, sig)
aprob <- min(1, (dgamma(can, a, b)/dgamma(x,
a, b))/(dnorm(can, mu, sig)/dnorm(x,
mu, sig)))
u <- runif(1)
if (u < aprob)
x <- can
vec[i] <- x
畫圖
設置參數。
nrep<- 54000 burnin<- 4000 shape<- 2.5 rate<-2.6
修改圖，僅包含預燒期後的鏈
vec=vec[-(1:burnin)] #vec=vec[burnin:length(vec)]
par(mfrow=c(2,1)) # 更改主框架，在一幀中有多少個圖形 plot(ts(vec), xlab="Chain", ylab="Draws") abline(h = mean(vec), lwd="2", col="red" )
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.007013 0.435600 0.724800 0.843300 1.133000 3.149000
var(vec[-(1:burnin)])
[1] 0.2976507
初始值
第一個樣本 vec 是我們鏈的初始/起始值。我們可以更改它，以查看收斂是否發生瞭變化。
x <- 3*a/b
vec[1] <- x
選擇方案
如果候選密度與目標分佈P（x）的形狀匹配，即q（x’| xt）≈P（x’）q（x’|），則該算法效果最佳。 xt）≈P（x’）。如果使用正態候選密度q，則在預燒期間必須調整方差參數σ2。
通常，這是通過計算接受率來完成的，接受率是在最後N個樣本的窗口中接受的候選樣本的比例。
如果σ2太大，則接受率將非常低，因為候選可能落在概率密度低得多的區域中，因此a1將非常小，且鏈將收斂得非常慢。
示例2：回歸的貝葉斯估計
Metropolis-Hastings采樣用於貝葉斯估計回歸模型。
設定參數
DGP和圖
# 創建獨立的x值，大約為零 x <- (-(Size-1)/2):((Size-1)/2) # 根據ax + b + N（0，sd）創建相關值 y <- trueA * x + trueB + rnorm(n=Size,mean=0,sd=trueSd)
正態分佈擬然
pred = a*x + b singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T) sumll = sum(singlelikelihoods)
為什麼使用對數
似然函數中概率的對數，這也是我求和所有數據點的概率（乘積的對數等於對數之和）的原因。
我們為什麼要做這個？強烈建議這樣做，因為許多小概率相乘的概率會變得很小。在某個階段，計算機程序會陷入數值四舍五入或下溢問題。
因此， 當您編寫概率時，請始終使用對數
示例：繪制斜率a的似然曲線
# 示例：繪制斜率a的似然曲線 plot (seq(3, 7, by=.05), slopelikelihoods , type="l")
先驗分佈
這三個參數的均勻分佈和正態分佈。
# 先驗分佈 # 更改優先級，log為True，因此這些均為log density/likelihood aprior = dunif(a, min=0, max=10, log = T) bprior = dnorm(b, sd = 2, log = T) sdprior = dunif(sd, min=0, max=30, log = T)
後驗
先驗和概率的乘積是MCMC將要處理的實際量。此函數稱為後驗函數。同樣，這裡我們使用和，因為我們使用對數。
posterior <- function(param){
return (likelihood(param) + prior(param))
}
Metropolis算法
該算法是從 後驗密度中采樣最常見的貝葉斯統計應用之一 。
- 上面定義的後驗。
- 從隨機參數值開始
- 根據某個候選函數的概率密度，選擇一個接近舊值的新參數值
- 以概率p（new）/ p（old）跳到這個新點，其中p是目標函數，並且p> 1也意味著跳躍
- 請註意，我們有一個 對稱的跳躍/候選分佈 q（x’| x）。
標準差σ是固定的。
所以接受概率等於
######## Metropolis 算法 ################
for (i in 1:iterations){
probab = exp(posterior(proposal) - posterior(chain[i,]))
if (runif(1) < probab){
chain[i+1,] = proposal
}else{
chain[i+1,] = chain[i,]
}
實施
（e）輸出接受的值，並解釋。
chain = metrMCMC(startvalue, 5500) burnIn = 5000 accep = 1-mean(duplicated(chain[-(1:burnIn),]))
算法的第一步可能會因初始值而有偏差，因此通常會被丟棄來進行進一步分析（預燒期）。令人感興趣的輸出是接受率：候選多久被算法接受拒絕一次？候選函數會影響接受率：通常，候選越接近，接受率就越大。但是，非常高的接受率通常是無益的：這意味著算法在同一點上“停留”，這導致對參數空間（混合）的處理不夠理想。
我們還可以更改初始值，以查看其是否更改結果/是否收斂。
startvalue = c(4,0,10)
小結
V1 V2 V3 Min. :4.068 Min. :-6.7072 Min. : 6.787 1st Qu.:4.913 1st Qu.:-2.6973 1st Qu.: 9.323 Median :5.052 Median :-1.7551 Median :10.178 Mean :5.052 Mean :-1.7377 Mean :10.385 3rd Qu.:5.193 3rd Qu.:-0.8134 3rd Qu.:11.166 Max. :5.989 Max. : 4.8425 Max. :19.223
#比較: summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-22.259 -6.032 -1.718 6.955 19.892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.1756 1.7566 -1.808 0.081 .
x 5.0469 0.1964 25.697 <2e-16 ***
---
Signif. codes: 0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1
Residual standard error: 9.78 on 29 degrees of freedom
Multiple R-squared: 0.9579, Adjusted R-squared: 0.9565
F-statistic: 660.4 on 1 and 29 DF, p-value: < 2.2e-16
summary(lm(y~x))$sigma
[1] 9.780494
coefficients(lm(y~x))[1]
(Intercept) -3.175555
coefficients(lm(y~x))[2]
x 5.046873
總結：
### 總結: ####################### par(mfrow = c(2,3)) hist(chain[-(1:burnIn),1],prob=TRUE,nclass=30,col="109" abline(v = mean(chain[-(1:burnIn),1]), lwd="2")
到此這篇關於R語言MCMC:Metropolis-Hastings采樣用於回歸的貝葉斯估計的文章就介紹到這瞭,更多相關R語言MCMC:Metropolis-Hastings采樣用於回歸的貝葉斯估計內容請搜索WalkonNet以前的文章或繼續瀏覽下面的相關文章希望大傢以後多多支持WalkonNet！