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

最近我們被客戶要求撰寫關于Metropolis-Hastings采樣的研究報告，包括一些圖形和統(tǒng)計輸出。

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）。過渡核是從狀態(tài)xi遷移到狀態(tài)xj的概率。

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

Metropolis算法

?對于一個Markov鏈是平穩(wěn)的?；旧媳硎?/p>

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

或者

方法是將轉換分為兩個子步驟；候選和接受拒絕。

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

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

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

該?轉移概率?可以寫成：

插入上一個方程式，我們有

Metropolis-Hastings算法?

A的選擇遵循以下邏輯。

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

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

即，當接受度大于1時，我們總是接受，而當接受度小于1時，我們將相應地拒絕。因此，Metropolis-Hastings算法包含以下內容：

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

2. 根據q（x'| x）隨機選擇一個新狀態(tài)x';

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

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

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

原則上，我們從分布P（x）提取保存的狀態(tài)，因為步驟4保證它們是不相關的。必須根據候選分布等不同因素來選擇T的值。 重要的是，尚不清楚應該使用哪種分布q（x'| x）；必須針對當前的特定問題進行調整。

屬性

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

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

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

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

Metropolis采樣

一個簡單的Metropolis-Hastings采樣

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

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

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

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

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

3. 計算“接受概率”

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

點擊標題查閱往期內容

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

左右滑動查看更多

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

初始值

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

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

選擇方案

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

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

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

示例2：回歸的貝葉斯估計

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

設定參數

DGP和圖

# 創(chuàng)建獨立的x值，大約為零 x <- (-(Size-1)/2):((Size-1)/2) # 根據ax + b + N（0，sd）創(chuàng)建相關值 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)

為什么使用對數

似然函數中概率的對數，這也是我求和所有數據點的概率（乘積的對數等于對數之和）的原因。

我們?yōu)槭裁匆鲞@個？強烈建議這樣做，因為許多小概率相乘的概率會變得很小。在某個階段，計算機程序會陷入數值四舍五入或下溢問題。

因此，?當您編寫概率時，請始終使用對數

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

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

先驗分布

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

# 先驗分布 # 更改優(yōu)先級，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算法

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

上面定義的后驗。

1. 從隨機參數值開始

2. 根據某個候選函數的概率密度，選擇一個接近舊值的新參數值

3. 以概率p（new）/ p（old）跳到這個新點，其中p是目標函數，并且p> 1也意味著跳躍

4. 請注意，我們有一個?對稱的跳躍/?候選分布?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-16summary(lm(y~x))$sigma[1] 9.780494coefficients(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采樣用于回歸的貝葉斯估計》。

點擊標題查閱往期內容

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