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

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

在貝葉斯方法中，馬爾可夫鏈蒙特卡羅方法尤其神秘。它們肯定是數(shù)學(xué)繁重且計(jì)算量大的過程，但它們背后的基本推理，就像數(shù)據(jù)科學(xué)中的許多其他東西一樣，可以變得直觀。這就是我的目標(biāo)。

相關(guān)視頻：馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實(shí)現(xiàn)|數(shù)據(jù)分享

馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實(shí)現(xiàn)

，時(shí)長08:47

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

MCMC 方法用于通過概率空間中的隨機(jī)抽樣來近似感興趣參數(shù)的后驗(yàn)分布。

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

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

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

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

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

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

在貝葉斯統(tǒng)計(jì)中，代表我們對參數(shù)的信念的分布稱為先驗(yàn)分布，因?yàn)樗诳吹饺魏螖?shù)據(jù)之前捕獲了我們的信念。似然分布通過表示一系列參數(shù)值以及每個(gè)參數(shù)解釋我們正在觀察的數(shù)據(jù)的可能性來總結(jié)觀察到的數(shù)據(jù)告訴我們的內(nèi)容。估計(jì)最大化似然分布的參數(shù)值只是回答了這個(gè)問題：什么參數(shù)值最有可能觀察到我們觀察到的數(shù)據(jù)？在沒有先驗(yàn)信念的情況下，我們可能會(huì)止步于此。

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

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

當(dāng)先驗(yàn)可能性結(jié)合起來時(shí)，數(shù)據(jù)（由可能性表示）支配了在巨人中長大的假設(shè)個(gè)人的弱先驗(yàn)信念。盡管那個(gè)人仍然認(rèn)為人類的平均身高比數(shù)據(jù)告訴他的要高一些，但他基本上相信數(shù)據(jù)。

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

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

MCMC 方法

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

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

假設(shè)我們想估計(jì)圓的面積：

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

由于 20 個(gè)點(diǎn)中有 15 個(gè)位于圓內(nèi)，因此該圓看起來約為 0.75 。對于只有 20 個(gè)隨機(jī)點(diǎn)的蒙特卡洛模擬來說還不錯(cuò)。

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

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

回想一下，我們正在嘗試估計(jì)我們感興趣的參數(shù)的后驗(yàn)分布，即人類平均身高：

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

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

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

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

紅點(diǎn)是隨機(jī)參數(shù)樣本

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

藍(lán)點(diǎn)僅代表任意時(shí)間點(diǎn)之后的隨機(jī)樣本，此時(shí)預(yù)計(jì)會(huì)發(fā)生收斂。注意：垂直堆疊點(diǎn)純粹是為了說明目的。

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

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

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

什么是MCMC，什么時(shí)候使用它？

MCMC只是一個(gè)從分布抽樣的算法。

這只是眾多算法之一。這個(gè)術(shù)語代表“馬爾可夫鏈蒙特卡洛”，因?yàn)樗且环N使用“馬爾可夫鏈”（我們將在后面討論）的“蒙特卡羅”（即隨機(jī)）方法。MCMC只是蒙特卡洛方法的一種，盡管可以將許多其他常用方法看作是MCMC的簡單特例。

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

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

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

如果這些樣本（查看文末了解數(shù)據(jù)獲取方式）是來自分布的獨(dú)立樣本，則 估計(jì)均值將會(huì)收斂在真實(shí)均值上。

假設(shè)我們的目標(biāo)分布是一個(gè)具有均值m和標(biāo)準(zhǔn)差的正態(tài)分布s。

作為一個(gè)例子，考慮用均值m和標(biāo)準(zhǔn)偏差s來估計(jì)正態(tài)分布的均值（在這里，我將使用對應(yīng)于標(biāo)準(zhǔn)正態(tài)分布的參數(shù)）：

我們可以很容易地使用這個(gè)rnorm 函數(shù)從這個(gè)分布中抽樣

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

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

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

事實(shí)上，在這種情況下，$ n $樣本估計(jì)的預(yù)期方差是$ 1 / n $，所以我們預(yù)計(jì)大部分值在$ \ 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

這個(gè)函數(shù)計(jì)算累積平均值之和。

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

將x軸轉(zhuǎn)換為對數(shù)坐標(biāo)并顯示另外30個(gè)隨機(jī)方法：

可以從您的一系列采樣點(diǎn)中抽取樣本分位數(shù)。

這是分析計(jì)算的點(diǎn)，其概率密度的2.5％低于：

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

我們可以通過在這種情況下的直接整合來估計(jì)這個(gè)

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積分估計(jì)點(diǎn)：

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

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

?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

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

在一個(gè)等級(jí)模型中，你可能會(huì)有大量的隨機(jī)效應(yīng)項(xiàng)被擬合，但是你最想對一個(gè)參數(shù)做出推論。在

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

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

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

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

馬爾可夫鏈蒙特卡羅

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

定義

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

?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

這個(gè)函數(shù)采用一個(gè)狀態(tài)向量x（其中x[i]是處于狀態(tài)的概率i），并通過將其與轉(zhuǎn)移矩陣相乘來迭代它P，使系統(tǒng)前進(jìn)到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函數(shù)來提取系統(tǒng)的主要特征向量（t()這里轉(zhuǎn)置矩陣以便得到左特征向量）。

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

然后在之前的數(shù)字上加上點(diǎn)，表明我們有多接近收斂：

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

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

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

這鏈條運(yùn)行了100個(gè)步驟：

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

繪制我們在每個(gè)狀態(tài)隨時(shí)間變化的時(shí)間分?jǐn)?shù)，而不是繪制狀態(tài)：

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

再運(yùn)行一下（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)

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

Metropolis算法

這是最簡單的MCMC算法。

MCMC采樣1d（單參數(shù)）問題

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

這里是一些參數(shù)和目標(biāo)密度的定義。

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

概率密度繪制

我們來定義一個(gè)非常簡單的算法，該算法從以當(dāng)前點(diǎn)為中心的標(biāo)準(zhǔn)偏差為4的正態(tài)分布中抽樣

而這只需要運(yùn)行MCMC的幾個(gè)步驟。它將從點(diǎn)x返回一個(gè)矩陣，其nsteps行數(shù)和列數(shù)與x元素的列數(shù)相同。如果在標(biāo)量上運(yùn)行， x它將返回一個(gè)向量。

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

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

?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)

運(yùn)行更長時(shí)間，結(jié)果開始看起來更好：

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

現(xiàn)在，運(yùn)行不同的方案 - 一個(gè)標(biāo)準(zhǔn)差很大（33），另一個(gè)標(biāo)準(zhǔn)差很?。?）。

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

注意三條軌跡正在移動(dòng)的不同方式。

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

藍(lán)色的蹤跡提出了傾向于被接受的小動(dòng)作，但是它隨著大部分的軌跡隨機(jī)行走。它需要數(shù)百次迭代才能達(dá)到概率密度的大部分。

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

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

由此可以計(jì)算獨(dú)立樣本的有效數(shù)量：

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

這更清楚地顯示了鏈條運(yùn)行時(shí)間更長的情況：

?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在兩個(gè)維度

給出了一個(gè)多元正態(tài)密度，給定一個(gè)均值向量（分布的中心）和方差 - 協(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)}}

如上所述，將目標(biāo)密度定義為兩個(gè)mvns的總和（這次未加權(quán)）：

?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)分布取樣也相當(dāng)簡單，但我們將使用MCMC從中抽取樣本。

這里有一些不同的策略 - 我們可以同時(shí)在兩個(gè)維度上提出動(dòng)作，或者我們可以獨(dú)立地沿著每個(gè)軸進(jìn)行采樣。這兩種策略都能奏效，雖然它們的混合速度會(huì)有所不同。

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

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

例如，參數(shù)1 的邊際分布是多少？

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

我們需要整合第一個(gè)參數(shù)的第二個(gè)參數(shù)的所有可能值。那么，因?yàn)槟繕?biāo)函數(shù)本身并不是標(biāo)準(zhǔn)化的，所以我們必須將其分解為一維積分值 。

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

數(shù)據(jù)獲取

在下面公眾號(hào)后臺(tái)回復(fù)“MCMC數(shù)據(jù)”，可獲取完整數(shù)據(jù)。

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

獲取全文完整資料。

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

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

R語言實(shí)現(xiàn)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ù)測心臟病

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 軟件對學(xué)術(shù)能力測驗(yàn)建立層次（分層）貝葉斯模型

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

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

R語言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究

R語言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測模型

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對多元隨機(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ì)