MCMC是從復(fù)雜概率模型中采樣的通用技術(shù)。

1. 蒙特卡洛
2. 馬爾可夫鏈
3. Metropolis-Hastings算法

問題

如果需要計(jì)算有復(fù)雜后驗(yàn)pdf p（θ| y）的隨機(jī)變量θ的函數(shù)f（θ）的平均值或期望值。

您可能需要計(jì)算后驗(yàn)概率分布p（θ）的最大值。

解決期望值的一種方法是從p（θ）繪制N個(gè)隨機(jī)樣本，當(dāng)N足夠大時(shí)，我們可以通過以下公式逼近期望值或最大值

將相同的策略應(yīng)用于通過從p（θ| y）采樣并取樣本集中的最大值來找到argmaxp（θ| y）。

解決方法

1.1直接模擬

1.2逆CDF

1.3拒絕/接受抽樣

如果我們不知道精確/標(biāo)準(zhǔn)化的pdf或非常復(fù)雜，則MCMC會(huì)派上用場(chǎng)。

馬爾可夫鏈

為了模擬馬爾可夫鏈，我們必須制定一個(gè) 過渡核T（xi，xj）。過渡核是從狀態(tài)xi遷移到狀態(tài)xj的概率。

 馬爾可夫鏈的收斂性意味著它具有平穩(wěn)分布π。馬爾可夫鏈的統(tǒng)計(jì)分布是平穩(wěn)的,那么它意味著分布不會(huì)隨著時(shí)間的推移而改變。

Metropolis算法

 對(duì)于一個(gè)Markov鏈?zhǔn)瞧椒€(wěn)的。基本上表示

處于狀態(tài)x并轉(zhuǎn)換為狀態(tài)x'的概率必須等于處于狀態(tài)x'并轉(zhuǎn)換為狀態(tài)x的概率

或者

方法是將轉(zhuǎn)換分為兩個(gè)子步驟；候選和接受拒絕。

令q（x'| x）表示 候選密度，我們可以使用概率 α（x'| x）來調(diào)整q  。

候選分布 Q（X'| X）是給定的候選X的狀態(tài)X'的條件概率，

和 接受分布 α（x'| x）的條件概率接受候選的狀態(tài)X'-X'。我們?cè)O(shè)計(jì)了接受概率函數(shù)，以滿足詳細(xì)的平衡。

該 轉(zhuǎn)移概率 可以寫成：

插入上一個(gè)方程式，我們有

Metropolis-Hastings算法 

A的選擇遵循以下邏輯。

在q下從x到x'的轉(zhuǎn)移太頻繁了。因此，我們應(yīng)該選擇α（x | x'）=1。但是，為了滿足 細(xì)致平穩(wěn)，我們有

下一步是選擇滿足上述條件的接受。Metropolis-Hastings是一種常見的 選擇：

即，當(dāng)接受度大于1時(shí)，我們總是接受，而當(dāng)接受度小于1時(shí)，我們將相應(yīng)地拒絕。因此，Metropolis-Hastings算法包含以下內(nèi)容：

初始化：隨機(jī)選擇一個(gè)初始狀態(tài)x；

根據(jù)q（x'| x）隨機(jī)選擇一個(gè)新狀態(tài)x';

3.接受根據(jù)α（x'| x）的狀態(tài)。如果不接受，則不會(huì)進(jìn)行轉(zhuǎn)移，因此無需更新任何內(nèi)容。否則，轉(zhuǎn)移為x'；

4.轉(zhuǎn)移到2，直到生成T狀態(tài)；

5.保存狀態(tài)x，執(zhí)行2。

原則上，我們從分布P（x）提取保存的狀態(tài)，因?yàn)椴襟E4保證它們是不相關(guān)的。必須根據(jù)候選分布等不同因素來選擇T的值。 重要的是，尚不清楚應(yīng)該使用哪種分布q（x'| x）；必須針對(duì)當(dāng)前的特定問題進(jìn)行調(diào)整。

屬性

Metropolis-Hastings算法的一個(gè)有趣特性是它 僅取決于比率

是候選樣本x'與先前樣本xt之間的概率，

是兩個(gè)方向（從xt到x'，反之亦然）的候選密度之比。如果候選密度對(duì)稱，則等于1。

馬爾可夫鏈從任意初始值x0開始，并且算法運(yùn)行多次迭代，直到“初始狀態(tài)”被“忘記”為止。這些被丟棄的樣本稱為預(yù)燒（burn-in）。其余的x可接受值集代表分布P（x）中的樣本

Metropolis采樣

一個(gè)簡(jiǎn)單的Metropolis-Hastings采樣

讓我們看看從 伽瑪分布 模擬任意形狀和比例參數(shù)，使用具有Metropolis-Hastings采樣算法。

下面給出了Metropolis-Hastings采樣器的函數(shù)。該鏈初始化為零，并在每個(gè)階段都建議使用N（a / b，a /（b * b））個(gè)候選對(duì)象。

基于正態(tài)分布且均值和方差相同gamma的Metropolis-Hastings獨(dú)立采樣

從某種狀態(tài)開始xt。代碼中的x。在代碼中提出一個(gè)新的狀態(tài)x'候選計(jì)算“接受概率”

從[0,1] 得出一些均勻分布的隨機(jī)數(shù)u；如果u <α接受該點(diǎn)，則設(shè)置xt + 1 = x'。否則，拒絕它并設(shè)置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

畫圖

設(shè)置參數(shù)。

nrep<- 54000
burnin<- 4000
shape<- 2.5
rate<-2.6

修改圖，僅包含預(yù)燒期后的鏈

vec=vec[-(1:burnin)]
#vec=vec[burnin:length(vec)]

par(mfrow=c(2,1)) # 更改主框架，在一幀中有多少個(gè)圖形
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

初始值

第一個(gè)樣本 vec 是我們鏈的初始/起始值。我們可以更改它，以查看收斂是否發(fā)生了變化。

  x <- 3*a/b
    vec[1] <- x

選擇方案

如果候選密度與目標(biāo)分布P（x）的形狀匹配，即q（x'| xt）≈P（x'）q（x'|），則該算法效果最佳。 xt）≈P（x'）。如果使用正態(tài)候選密度q，則在預(yù)燒期間必須調(diào)整方差參數(shù)σ2。

通常，這是通過計(jì)算接受率來完成的，接受率是在最后N個(gè)樣本的窗口中接受的候選樣本的比例。

如果σ2太大，則接受率將非常低，因?yàn)楹蜻x可能落在概率密度低得多的區(qū)域中，因此a1將非常小，且鏈將收斂得非常慢。

示例2：回歸的貝葉斯估計(jì)

Metropolis-Hastings采樣用于貝葉斯估計(jì)回歸模型。

設(shè)定參數(shù)

DGP和圖

# 創(chuàng)建獨(dú)立的x值，大約為零
x <- (-(Size-1)/2):((Size-1)/2)
# 根據(jù)ax + b + N（0，sd）創(chuàng)建相關(guān)值
y <- trueA * x + trueB + rnorm(n=Size,mean=0,sd=trueSd)

正態(tài)分布擬然

pred = a*x + b
  singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
  sumll = sum(singlelikelihoods)

為什么使用對(duì)數(shù)

似然函數(shù)中概率的對(duì)數(shù)，這也是我求和所有數(shù)據(jù)點(diǎn)的概率（乘積的對(duì)數(shù)等于對(duì)數(shù)之和）的原因。

我們?yōu)槭裁匆鲞@個(gè)？強(qiáng)烈建議這樣做，因?yàn)樵S多小概率相乘的概率會(huì)變得很小。在某個(gè)階段，計(jì)算機(jī)程序會(huì)陷入數(shù)值四舍五入或下溢問題。

因此， 當(dāng)您編寫概率時(shí)，請(qǐng)始終使用對(duì)數(shù)

示例：繪制斜率a的似然曲線

# 示例：繪制斜率a的似然曲線
plot (seq(3, 7, by=.05), slopelikelihoods , type="l")

先驗(yàn)分布

這三個(gè)參數(shù)的均勻分布和正態(tài)分布。

# 先驗(yàn)分布
 
# 更改優(yōu)先級(jí)，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)

后驗(yàn)

先驗(yàn)和概率的乘積是MCMC將要處理的實(shí)際量。此函數(shù)稱為后驗(yàn)函數(shù)。同樣，這里我們使用和，因?yàn)槲覀兪褂脤?duì)數(shù)。

posterior <- function(param){
  return (likelihood(param) + prior(param))
}

Metropolis算法

該算法是從 后驗(yàn)密度中采樣最常見的貝葉斯統(tǒng)計(jì)應(yīng)用之一 。

1. 上面定義的后驗(yàn)。
2. 從隨機(jī)參數(shù)值開始
3. 根據(jù)某個(gè)候選函數(shù)的概率密度，選擇一個(gè)接近舊值的新參數(shù)值
4. 以概率p（new）/ p（old）跳到這個(gè)新點(diǎn)，其中p是目標(biāo)函數(shù)，并且p> 1也意味著跳躍
5. 請(qǐng)注意，我們有一個(gè) 對(duì)稱的跳躍/候選分布 q（x'| x）。

標(biāo)準(zhǔn)差σ是固定的。

所以接受概率等于

######## 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,]
    }

實(shí)施

（e）輸出接受的值，并解釋。

chain = metrMCMC(startvalue, 5500)
 
burnIn = 5000
accep = 1-mean(duplicated(chain[-(1:burnIn),]))

算法的第一步可能會(huì)因初始值而有偏差，因此通常會(huì)被丟棄來進(jìn)行進(jìn)一步分析（預(yù)燒期）。令人感興趣的輸出是接受率：候選多久被算法接受拒絕一次？候選函數(shù)會(huì)影響接受率：通常，候選越接近，接受率就越大。但是，非常高的接受率通常是無益的：這意味著算法在同一點(diǎn)上“停留”，這導(dǎo)致對(duì)參數(shù)空間（混合）的處理不夠理想。

我們還可以更改初始值，以查看其是否更改結(jié)果/是否收斂。

startvalue = c(4,0,10)

小結(jié)

 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

總結(jié)：

### 總結(jié): #######################
 
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")

到此這篇關(guān)于R語(yǔ)言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)的文章就介紹到這了,更多相關(guān)R語(yǔ)言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論