將相同的策略應用於通過從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算法包含以下內容：

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獨立採樣

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 是我們鏈的初始/起始值。我們可以更改它，以查看收斂是否發生了變化。

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算法

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

上面定義的後驗。

標準差σ是固定的。

所以接受機率等於

######## 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採樣用於回歸的貝葉斯估計》。

點擊標題查閱往期內容

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採樣算法的實現

Metropolis Hastings採樣和貝葉斯泊松回歸Poisson模型

Matlab用BUGS馬爾可夫區制轉換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H採樣分析時間序列R語言RSTAN MCMC：NUTS採樣算法用LASSO 構建貝葉斯線性回歸模型分析職業聲望數據

R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉換隨機波動率SV模型、粒子濾波、Metropolis Hasting採樣時間序列分析

R語言Metropolis Hastings採樣和貝葉斯泊松回歸Poisson模型

R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數據和可視化診斷

R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs採樣算法實例

R語言貝葉斯Poisson泊松-正態分布模型分析職業足球比賽進球數

R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數

R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病

R語言中貝葉斯網絡（BN）、動態貝葉斯網絡、線性模型分析錯頜畸形數據

R語言中的block Gibbs吉布斯採樣貝葉斯多元線性回歸

Python貝葉斯回歸分析住房負擔能力數據集

R語言實現貝葉斯分位數回歸、lasso和自適應lasso貝葉斯分位數回歸分析

Python用PyMC3實現貝葉斯線性回歸模型

R語言用WinBUGS 軟體對學術能力測驗建立層次（分層）貝葉斯模型

R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析

R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數據

R語言基於copula的貝葉斯分層混合模型的診斷準確性研究

R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型

R語言貝葉斯推斷與MCMC：實現Metropolis-Hastings 採樣算法示例

R語言stan進行基於貝葉斯推斷的回歸模型

R語言中RStan貝葉斯層次模型分析示例

R語言使用Metropolis-Hastings採樣算法自適應貝葉斯估計與可視化

R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型

WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較

R語言實現MCMC中的Metropolis–Hastings算法與吉布斯採樣

R語言貝葉斯推斷與MCMC：實現Metropolis-Hastings 採樣算法示例

R語言使用Metropolis-Hastings採樣算法自適應貝葉斯估計與可視化

視頻：R語言中的Stan機率編程MCMC採樣的貝葉斯模型

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