全文鏈接：http://tecdat.cn/?p=19664?

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

1. 蒙特卡洛

2. 馬爾可夫鏈

3. Metropolis-Hastings算法（點(diǎn)擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)）。

問題

如果需要計(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）。

相關(guān)視頻

解決方法

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)容：

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

2. 根據(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ú)立采樣

1. 從某種狀態(tài)開始xt。代碼中的x。

2. 在代碼中提出一個(gè)新的狀態(tài)x'候選

3. 計(jì)算“接受概率”

1. 從[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" )


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

Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法

左右滑動(dòng)查看更多

01

02

03

04

? ?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)用之一?。

上面定義的后驗(yàn)。

1. 從隨機(jī)參數(shù)值開始

2. 根據(jù)某個(gè)候選函數(shù)的概率密度，選擇一個(gè)接近舊值的新參數(shù)值

3. 以概率p（new）/ p（old）跳到這個(gè)新點(diǎn)，其中p是目標(biāo)函數(shù)，并且p> 1也意味著跳躍

4. 請(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")


點(diǎn)擊文末“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)》。

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

Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法
R語言貝葉斯METROPOLIS-HASTINGS GIBBS 吉布斯采樣器估計(jì)變點(diǎn)指數(shù)分布分析泊松過程車站等待時(shí)間
R語言馬爾可夫MCMC中的METROPOLIS HASTINGS，MH算法抽樣（采樣）法可視化實(shí)例
python貝葉斯隨機(jī)過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實(shí)現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動(dòng)率模型、序列蒙特卡羅SMC、M H采樣分析時(shí)間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構(gòu)建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動(dòng)率SV模型、粒子濾波、Metropolis Hasting采樣時(shí)間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)
R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測(cè)心臟病
R語言中貝葉斯網(wǎng)絡(luò)（BN）、動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯(cuò)頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集
R語言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測(cè)選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究
R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測(cè)模型
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進(jìn)行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型
WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較
R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)