原文鏈接：http://tecdat.cn/?p=2687

在貝葉斯方法中，馬爾可夫鏈蒙特卡羅方法尤其神秘

它們肯定是數學繁重且計算量大的過程，但它們背后的基本推理，就像數據科學中的許多其他東西一樣，可以變得直觀。這就是我的目標。

那么，什么是馬爾可夫鏈蒙特卡羅（MCMC）方法？簡短的回答是：

MCMC 方法用于通過概率空間中的隨機抽樣來近似感興趣參數的后驗分布。

在這篇文章中，我將解釋這個簡短的答案。

首先，一些術語。感興趣的參數只是總結我們感興趣的現(xiàn)象的一些數字。通常我們使用統(tǒng)計數據來估計參數。例如，如果我們想了解成年人的身高，我們感興趣的參數可能是平均身高。分布是我們參數的每個可能值的數學表示，以及我們觀察每個值的可能性。最著名的例子是鐘形曲線：

在貝葉斯統(tǒng)計方法中，分布有額外的解釋。貝葉斯不僅表示參數的值以及每個參數成為真實值的可能性，而是將分布視為描述我們對參數的信念。因此，上面的鐘形曲線表明我們非常確定參數的值非常接近于零，但我們認為真實值高于或低于該值的可能性相同，直到某個點。

碰巧的是，人類身高確實遵循正態(tài)曲線，所以假設我們相信人類平均身高的真實值遵循如下鐘形曲線：

顯然，這張圖所代表的有信仰的人多年來一直生活在巨人中間，因為據他們所知，最有可能的成人平均身高是1米8（但他們并不是特別自信）。

讓我們想象一下這個人去收集了一些數據，他們觀察了1米6到 1米8之間的人群。我們可以在下面表示該數據，以及另一條顯示人類平均身高值最能解釋數據的正態(tài)曲線：

在貝葉斯統(tǒng)計中，代表我們對參數的信念的分布稱為先驗分布，因為它在看到任何數據之前捕獲了我們的信念。似然分布通過表示一系列參數值以及每個參數解釋我們正在觀察的數據的可能性來總結觀察到的數據告訴我們的內容。估計最大化似然分布的參數值只是回答了這個問題：什么參數值最有可能觀察到我們觀察到的數據？在沒有先驗信念的情況下，我們可能會止步于此。

然而，貝葉斯分析的關鍵是結合先驗分布和似然分布來確定后驗分布。這告訴我們，考慮到我們先前的信念，哪些參數值可以最大限度地提高觀察我們所做的特定數據的機會。在我們的例子中，后驗分布如下所示：

上圖，紅線代表后驗分布。您可以將其視為先驗分布和似然分布的一種平均值。由于先驗分布更短且分布更廣，因此它代表了一組對人類平均身高的真實值“不太確定”的信念。同時，似然總結了相對狹窄范圍內的數據，因此它代表了對真實參數值的“更確定”的猜測。

當先驗可能性結合起來時，數據（由可能性表示）支配了在巨人中長大的假設個人的弱先驗信念。盡管那個人仍然認為人類的平均身高比數據告訴他的要高一些，但他基本上相信數據。

在兩條鐘形曲線的情況下，求解后驗分布非常容易。有一個簡單的公式可以將兩者結合起來。但是，如果我們的先驗分布和似然分布沒有那么好怎么辦？有時，使用沒有方便形狀的分布對我們的數據或我們的先驗信念進行建模是最準確的。如果我們的概率最好由具有兩個峰值的分布來表示，并且出于某種原因我們想要解釋一些非常古怪的先驗分布怎么辦？我通過手繪一個丑陋的先驗分布來可視化下面的場景：

和以前一樣，存在一些后驗分布，它給出了每個參數值的可能性。但它有點難以看出它可能是什么樣子，并且不可能通過分析來解決。

MCMC 方法

MCMC 方法允許我們估計后驗分布的形狀，以防我們無法直接計算它?；叵胍幌拢琈CMC 代表馬爾可夫鏈蒙特卡羅方法。為了理解它們是如何工作的，我將介紹蒙特卡羅模擬。

蒙特卡羅模擬只是一種通過重復生成隨機數來估計固定參數的方法。通過獲取生成的隨機數并對它們進行一些計算，蒙特卡洛模擬提供了一個參數的近似值。

假設我們想估計圓的面積：

由于圓在邊長為 1的正方形內，因此面積可以很容易地計算為 0.785 。但是，我們可以在正方形內隨機放置 20 個點。然后我們計算落在圓圈內的點的比例，并將其乘以正方形的面積。這個數字是圓面積的一個很好的近似值。

由于 20 個點中有 15 個位于圓內，因此該圓看起來約為 0.75 。對于只有 20 個隨機點的蒙特卡洛模擬來說還不錯。

蒙特卡羅模擬不僅用于估計困難形狀的區(qū)域。通過生成大量隨機數，它們可用于對非常復雜的過程進行建模。

有了蒙特卡羅模擬和馬爾可夫鏈的一些知識，我希望對 MCMC 方法如何工作的無數學解釋非常直觀。

回想一下，我們正在嘗試估計我們感興趣的參數的后驗分布，即人類平均身高：

我不是可視化專家，顯然我也不擅長將我的示例保持在常識范圍內：我的后驗分布示例嚴重高估了人類的平均身高。

我們知道后驗分布在我們的先驗分布和似然分布的范圍內，但無論出于何種原因，我們都無法直接計算它。使用 MCMC 方法，我們將有效地從后驗分布中抽取樣本，然后計算統(tǒng)計數據，例如抽取樣本的平均值。

首先，MCMC 方法選擇一個隨機參數值來考慮。模擬將繼續(xù)生成隨機值（這是蒙特卡洛部分），但要遵守一些規(guī)則來確定什么是好的參數值。訣竅是，對于一對參數值，可以通過計算每個值解釋數據的可能性來計算哪個是更好的參數值，給定我們的先驗信念。如果隨機生成的參數值比上一個更好，則以一定的概率將其添加到參數值鏈中，該概率取決于它的好壞程度（這是馬爾可夫鏈部分）。

為了直觀地解釋這一點，讓我們回想一下某個值的分布高度代表觀察該值的概率。因此，我們可以認為我們的參數值（x 軸）展示了高概率和低概率的區(qū)域，顯示在 y 軸上。對于單個參數，MCMC 方法從沿 x 軸隨機采樣開始：

紅點是隨機參數樣本

由于隨機樣本受固定概率的影響，它們往往會在一段時間后收斂于我們感興趣的參數的最高概率區(qū)域：

藍點僅代表任意時間點之后的隨機樣本，此時預計會發(fā)生收斂。注意：垂直堆疊點純粹是為了說明目的。

收斂后，MCMC 采樣會產生一組點，這些點是來自后驗分布的樣本。圍繞這些點繪制直方圖，并計算您喜歡的任何統(tǒng)計數據：

在 MCMC 模擬生成的樣本集上計算的任何統(tǒng)計量都是我們對真實后驗分布統(tǒng)計量的最佳猜測。

MCMC 方法也可用于估計多個參數（例如人的身高和體重）的后驗分布。對于n 個參數，在 n 維空間中存在高概率區(qū)域，其中某些參數值集可以更好地解釋觀察到的數據。因此，我認為 MCMC 方法是在概率空間內隨機抽樣以近似后驗分布。

點擊標題查閱往期內容

R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

左右滑動查看更多

01

02

03

04

什么是MCMC，什么時候使用它？

MCMC只是一個從分布抽樣的算法。

這只是眾多算法之一。這個術語代表“馬爾可夫鏈蒙特卡洛”，因為它是一種使用“馬爾可夫鏈”（我們將在后面討論）的“蒙特卡羅”（即隨機）方法。MCMC只是蒙特卡洛方法的一種，盡管可以將許多其他常用方法看作是MCMC的簡單特例。

我為什么要從分布中抽樣？

從分布中抽取樣本是解決一些問題的最簡單的方法。

可能MCMC最常用的方法是從貝葉斯推理中的某個模型的后驗概率分布中抽取樣本。通過這些樣本，你可以問一些問題：“參數的平均值和可信度是多少？”。

如果這些樣本?（?查看文末了解數據獲取方式?）?是來自分布的獨立樣本，則 估計均值將會收斂在真實均值上。

假設我們的目標分布是一個具有均值m和標準差的正態(tài)分布s。

作為一個例子，考慮用均值m和標準偏差s來估計正態(tài)分布的均值（在這里，我將使用對應于標準正態(tài)分布的參數）：

我們可以很容易地使用這個rnorm 函數從這個分布中抽樣

?seasamples<-rn?000,m,s)

樣本的平均值非常接近真實平均值（零）：

?mean(sa?es) ?##?[1]?-0.?537

事實上，在這種情況下，n樣本估計的預期方差是1/n，所以我們預計大部分值在$ \ pm 2 \，/ \ sqrt {n} = 0.02 。

?summary(re?0,mean(rnorm(10000,m,s))))?##?Min.?1st?Qu.?Median?Mean?3rd?Qu.?Max.?##?-0.03250?-0.00580?0.00046?0.00042?0.00673?0.03550

這個函數計算累積平均值之和。

?cummean<-fun?msum(x)/seq_along(x) ?plot(cummaaSample",ylab="Cumulative?mean",panel.aabline(h=0,col="red"),las=1)

將x軸轉換為對數坐標并顯示另外30個隨機方法：

可以從您的一系列采樣點中抽取樣本分位數。

這是分析計算的點，其概率密度的2.5％低于：

?p<-0.025a.true<-qnorm(p,m,s)a.true1##?[1]?-1.96

我們可以通過在這種情況下的直接整合來估計這個

aion(x)dnorm(x,m,s) g<-function(a)integrate(f,-Inf,a)$valuea.int<-uniroot(function(x)g(a10,0))$roota.int1##?[1]?-1.96

并用Monte Carlo積分估計點：

a.mc<-unnasamples,p))a.mc##?[1]?-2.023a.true-a.mc##?[1]?0.06329

但是，在樣本量趨于無窮大的極限內，這將會收斂。此外，有可能就錯誤的性質作出陳述; 如果我們重復采樣過程100次，那么我們得到一系列與均值附近的誤差相同幅度的誤差的估計：

?a.mc<-replicate(anorm(10000,m,s),p))summary(a.true-a.mc) ?##?Min.?1st?Qu.?Median?Mean?3rd?Qu.?Max.?##?-0.05840?-0.01640?-0.00572?-0.00024?0.01400?0.07880

這種事情真的很常見。在大多數貝葉斯推理中，后驗分布是一些（可能很大的）參數向量的函數，您想對這些參數的子集進行推理。

在一個等級模型中，你可能會有大量的隨機效應項被擬合，但是你最想對一個參數做出推論。在

貝葉斯框架中，您可以計算您感興趣的參數在所有其他參數上的邊際分布（這是我們上面要做的）。

為什么“傳統(tǒng)統(tǒng)計”不使用蒙特卡洛方法？

對于傳統(tǒng)教學統(tǒng)計中的許多問題，不是從分布中抽樣，可以使函數最大化或最大化。所以我們需要一些函數來描述可能性并使其最大化（最大似然推理），或者一些計算平方和并使其最小化的函數。

然而，蒙特卡羅方法在貝葉斯統(tǒng)計中的作用與頻率統(tǒng)計中的優(yōu)化程序相同，這只是執(zhí)行推理的算法。所以，一旦你基本知道MCMC正在做什么，你可以像大多數人把他們的優(yōu)化程序當作黑匣子一樣對待它，像一個黑匣子。

馬爾可夫鏈蒙特卡羅

假設我們想要抽取一些目標分布，但是我們不能像從前那樣抽取獨立樣本。有一個使用馬爾科夫鏈蒙特卡洛（MCMC）來做這個的解決方案。首先，我們必須定義一些事情，以便下一句話是有道理的：我們要做的是試圖構造一個馬爾科夫鏈，它抽樣的目標分布作為它的平穩(wěn)分布。

定義

假設我們有一個三態(tài)馬爾科夫過程。讓我們P為鏈中的轉移概率矩陣：

?P<-rbind(a(.2,.1,.7),c(.25,.25,.5))P?##?[,1]?[,2]?[,3]##?[1,]?0.50?0.25?0.25##?[2,]?0.20?0.10?0.70##?[3,]?0.25?0.25?0.50?rowSums(P) ?##?[1]?1?1?1

P[i,j]給出了從狀態(tài)i到狀態(tài)的概率j。

請注意，與行不同，列不一定總和為1：

?colSums(P) ?##?[1]?0.95?0.60?1.45

這個函數采用一個狀態(tài)向量x（其中x[i]是處于狀態(tài)的概率i），并通過將其與轉移矩陣相乘來迭代它P，使系統(tǒng)前進到n步驟。

?iterate.P<-function(x,P,n){ res<-matrix(NA,n+1,len a<-xfor(iinseq_len(n)) res[i+1,]<-x<-x%*%P? res}

從處于狀態(tài)1的系統(tǒng)開始（x向量 [1,0,0] 也是如此，表示處于狀態(tài)1的概率為100％，不處于任何其他狀態(tài)）

同樣，對于另外兩種可能的起始狀態(tài)：

?y2<-iterate.P(c(0,1,0),P,n)y3<-iterate.P(c(0,0,1),P,n)

這表明了平穩(wěn)分布的收斂性。

ma=1,xlab="Step",ylab="y",las=1)matlines(0:n,y2,lty=2)matlines(0:n,y3,lty=3)

我們可以使用R的eigen函數來提取系統(tǒng)的主要特征向量（t()這里轉置矩陣以便得到左特征向量）。

?v<-eigen(t(P) ars[,1]v<-v/sum(v)#?歸一化特征向量

然后在之前的數字上加上點，表明我們有多接近收斂：

matplot(0:n,y1a3,lty=3)points(rep(10,3),v,col=1:3)

上面的過程迭代了不同狀態(tài)的總體概率; 而不是通過系統(tǒng)的實際轉換。所以，讓我們迭代系統(tǒng)，而不是概率向量。

?run<-function(i,P,n){ res<-integer(n)for(a(n)) res[[t]]<-i<-sample(nrow(P),1,pr=P[i,])? res}

這鏈條運行了100個步驟：

?samples<-run(1,P,100)ploaes,type="s",xlab="Step",ylab="State",las=1)

繪制我們在每個狀態(tài)隨時間變化的時間分數，而不是繪制狀態(tài)：

?plot(cummean(samplesa2)lines(cummean(samples==3),col=3)

再運行一下（5000步）

?n<-5000set.seed(1) samples<-run(1,P,n)plot(cummeanasamples==2),col=2)lines(cummean(samples==3),col=3)abline(h=v,lty=2,col=1:3)

所以這里的關鍵是：馬爾可夫鏈有一些不錯的屬性。馬爾可夫鏈有固定的分布，如果我們運行它們足夠長的時間，我們可以看看鏈條在哪里花費時間，并對該平穩(wěn)分布進行合理的估計。

Metropolis算法

這是最簡單的MCMC算法。

MCMC采樣1d（單參數）問題

這是兩個正態(tài)分布的加權和。這種分布相當簡單，可以從MCMC中抽取樣本。

這里是一些參數和目標密度的定義。

?p<-0.4ma1,2)sd<-c(.5,2)f<-function(x)p*dnora],sd[1])+(1-p)*dnorm(x,mu[2],sd[2])

概率密度繪制

我們來定義一個非常簡單的算法，該算法從以當前點為中心的標準偏差為4的正態(tài)分布中抽樣

而這只需要運行MCMC的幾個步驟。它將從點x返回一個矩陣，其nsteps行數和列數與x元素的列數相同。如果在標量上運行， x它將返回一個向量。

?run<-funagth(x))for(iinseq_len(nsteps)) res[i,]<-x<-step(x,f,q)drop(res)}

這里是馬爾可夫鏈的前1000步，目標密度在右邊：

?layout(matrix(ca,type="s",xpd=NA,ylab="Parameter",xlab="Sample",las=1) usr<-par("usr") xx<-seq(usr[a4],length=301)plot(f(xx),xx,type="l",yaxs="i",axes=FALSE,xlab="")

hist(res,5aALSE,main="",ylim=c(0,.4),las=1,xlab="x",ylab="Probability?density") z<-integrate(f,-Inf,Inf)$valuecurve(f(x)/z,add=TRUE,col="red",n=200)

運行更長時間，結果開始看起來更好：

res.long<-run(-10,f,q,50000)hist(res.long,100,freq=FALSE,main="",ylim=c(0,.4),las=1,xlab

現(xiàn)在，運行不同的方案 - 一個標準差很大（33），另一個標準差很小（3）。

res.fast<-run(-10action(x)rnorm(1,x,33),1000)res.slow<-run(-10,f,functanorm(1,x,.3),1000)

注意三條軌跡正在移動的不同方式。

相反，紅色的痕跡拒絕其中的大部分空間。

藍色的蹤跡提出了傾向于被接受的小動作，但是它隨著大部分的軌跡隨機行走。它需要數百次迭代才能達到概率密度的大部分。

您可以在隨后的參數中看到不同方案步驟在自相關中的效果 - 這些圖顯示了不同滯后步驟之間自相關系數的衰減，藍線表示統(tǒng)計獨立性。

?par(mfrow=c(1,3ain="Intermediate")acf(res.fast,las=1,m

由此可以計算獨立樣本的有效數量：

1coda::effectiveSize(res)1?2##?var1?##?1871coda::effectiveSize(res.fast)1?2##?var1?##?33.191coda::effectiveSize(res.slow)1?2##?var1?##?5.378

這更清楚地顯示了鏈條運行時間更長的情況：

?naun(-10,f,q,n)) xlim<-range(sapply(saa100) hh<-lapply(samples,function(x)hist(x,br,plot=FALSE)) ylim<-c(0,max(f(xx)))

顯示100,1,000,10,000和100,000步：

for(hinhh){plot(h,main="",freq=a=300)}

MCMC在兩個維度

給出了一個多元正態(tài)密度，給定一個均值向量（分布的中心）和方差 - 協(xié)方差矩陣。

?make.mvn<-function(mean,vcv){ logdet<-as.numeric(detea+logdet vcv.i<-solve(vcv)function(x){ dx<-x-meanexp(-(tmp+rowSums((dx%*%vcv.i)*dx))/2)}}

如上所述，將目標密度定義為兩個mvns的總和（這次未加權）：

?mu1<-c(-1,1)mu2<-c(2,-2)vCV1<-ma5,.25,1.5),2,2)vCV2<-matrix(c(2,-.5,-.5,2aunctioax)+f2(x)x<-seq(-5,6,length=71)y<-seq(-7,6,lena-expand.grid(x=x,y=y)z<-matrix(aaTRUE)

從多元正態(tài)分布取樣也相當簡單，但我們將使用MCMC從中抽取樣本。

這里有一些不同的策略 - 我們可以同時在兩個維度上提出動作，或者我們可以獨立地沿著每個軸進行采樣。這兩種策略都能奏效，雖然它們的混合速度會有所不同。

假設我們實際上并不知道如何從mvn中抽樣 ，讓我們提出一個在兩個維度上一致的提案分布，從每邊的寬度為“d”的正方形取樣。

比較抽樣分布與已知分布：

例如，參數1 的邊際分布是多少？

?hisales[,1],freq=FALSa",xlab="x",ylab="Probability?density")

我們需要整合第一個參數的第二個參數的所有可能值。那么，因為目標函數本身并不是標準化的，所以我們必須將其分解為一維積分值 。

?m<-function(x1){ g<-Vectorize(function(x2)f(c(x1,ae(g,-Inf,Inf)$value} xx<-seq(mina]),max(sales[,1]),length=201) yy<-suehist(samples[,1],freq=FALSE,ma,0.25))lines(xx,yy/z,col="red")

數據獲取

在下面公眾號后臺回復“MCMC數****據”，可獲取完整數據。

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言中實現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型》。

點擊標題查閱往期內容

R語言實現(xiàn)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采樣用于回歸的貝葉斯估計