詳解R語言MCMC:Metropolis-Hastings采樣用於回歸的貝葉斯估計

MCMC是從復雜概率模型中采樣的通用技術。

  1. 蒙特卡洛
  2. 馬爾可夫鏈
  3. 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算法

該算法是 後驗密度采樣最常見的貝葉斯統計應用之一 。

  1. 上面定義的後驗。
  2. 從隨機參數值開始
  3. 根據某個候選函數的概率密度,選擇一個接近舊值的新參數值
  4. 以概率p(new)/ p(old)跳到這個新點,其中p是目標函數,並且p> 1也意味著跳躍
  5. 請註意,我們有一個 對稱的跳躍/候選分佈 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!