共计 9543 个字符，预计需要花费 24 分钟才能阅读完成。
原文链接：http://tecdat.cn/?p=2687
在贝叶斯办法中，马尔可夫链蒙特卡罗办法尤其神秘
它们必定是数学沉重且计算量大的过程，但它们背地的根本推理，就像数据迷信中的许多其余货色一样，能够变得直观。这就是我的指标。
那么，什么是马尔可夫链蒙特卡罗（MCMC）办法？简短的答复是：
MCMC 办法用于通过概率空间中的随机抽样来近似感兴趣参数的后验散布。
在这篇文章中，我将解释这个简短的答案。
首先，一些术语。感兴趣的参数只是总结咱们感兴趣的景象的一些数字。通常咱们应用统计数据来预计参数。例如，如果咱们想理解成年人的身高，咱们感兴趣的参数可能是均匀身高。散布是咱们参数的每个可能值的数学示意，以及咱们察看每个值的可能性。最驰名的例子是钟形曲线：
在贝叶斯统计办法中，散布有额定的解释。贝叶斯不仅示意参数的值以及每个参数成为实在值的可能性，而是将散布视为形容咱们对参数的信念。因而，下面的钟形曲线表明咱们十分确定参数的值十分接近于零，但咱们认为实在值高于或低于该值的可能性雷同，直到某个点。
碰巧的是，人类身高的确遵循正态曲线，所以假如咱们置信人类均匀身高的实在值遵循如下钟形曲线：
显然，这张图所代表的有信奉的人多年来始终生存在伟人两头，因为据他们所知，最有可能的成人均匀身高是 1 米 8（但他们并不是特地自信）。
让咱们设想一下这个人去收集了一些数据，他们察看了 1 米 6 到 1 米 8 之间的人群。咱们能够在上面示意该数据，以及另一条显示人类均匀身高值最能解释数据的正态曲线：
在贝叶斯统计中，代表咱们对参数的信念的散布称为先验散布，因为它在看到任何数据之前捕捉了咱们的信念。似然散布通过示意一系列参数值以及每个参数解释咱们正在察看的数据的可能性来总结察看到的数据通知咱们的内容。预计最大化似然散布的参数值只是答复了这个问题：什么参数值最有可能察看到咱们察看到的数据？在没有先验信念的状况下，咱们可能会止步于此。
然而，贝叶斯剖析的要害是联合先验散布和似然散布来确定后验散布。这通知咱们，思考到咱们先前的信念，哪些参数值能够最大限度地进步察看咱们所做的特定数据的机会。在咱们的例子中，后验散布如下所示：
上图，红线代表后验散布。您能够将其视为先验散布和似然散布的一种平均值。因为先验散布更短且散布更广，因而它代表了一组对人类均匀身高的实在值“不太确定”的信念。同时，似然总结了绝对狭隘范畴内的数据，因而它代表了对实在参数值的“更确定”的猜想。
当先验可能性联合起来时，数据（由可能性示意）摆布了在伟人中长大的假如集体的弱先验信念。只管那个人依然认为人类的均匀身高比数据通知他的要高一些，但他基本上置信数据。
在两条钟形曲线的状况下，求解后验散布非常容易。有一个简略的公式能够将两者联合起来。然而，如果咱们的先验散布和似然散布没有那么好怎么办？有时，应用没有不便形态的散布对咱们的数据或咱们的先验信念进行建模是最精确的。如果咱们的概率最好由具备两个峰值的散布来示意，并且出于某种原因咱们想要解释一些十分乖僻的先验散布怎么办？我通过手绘一个俊俏的先验散布来可视化上面的场景：
和以前一样，存在一些后验散布，它给出了每个参数值的可能性。但它有点难以看出它可能是什么样子，并且不可能通过剖析来解决。
MCMC 办法
MCMC 办法容许咱们预计后验散布的形态，以防咱们无奈间接计算它。回忆一下，MCMC 代表马尔可夫链蒙特卡罗办法。为了了解它们是如何工作的，我将介绍蒙特卡罗模仿。
蒙特卡罗模仿只是一种通过反复生成随机数来预计固定参数的办法。通过获取生成的随机数并对它们进行一些计算，蒙特卡洛模仿提供了一个参数的近似值。
假如咱们想预计圆的面积：
因为圆在边长为 1 的正方形内，因而面积能够很容易地计算为 0.785。然而，咱们能够在正方形内随机搁置 20 个点。而后咱们计算落在圆圈内的点的比例，并将其乘以正方形的面积。这个数字是圆面积的一个很好的近似值。
因为 20 个点中有 15 个位于圆内，因而该圆看起来约为 0.75。对于只有 20 个随机点的蒙特卡洛模仿来说还不错。
蒙特卡罗模仿不仅用于预计艰难形态的区域。通过生成大量随机数，它们可用于对非常复杂的过程进行建模。
有了蒙特卡罗模仿和马尔可夫链的一些常识，我心愿对 MCMC 办法如何工作的无数学解释十分直观。
回忆一下，咱们正在尝试预计咱们感兴趣的参数的后验散布，即人类均匀身高：
我不是可视化专家，显然我也不善于将我的示例放弃在常识范畴内：我的后验散布示例重大高估了人类的均匀身高。
咱们晓得后验散布在咱们的先验散布和似然散布的范畴内，但无论出于何种起因，咱们都无奈间接计算它。应用 MCMC 办法，咱们将无效地从后验散布中抽取样本，而后计算统计数据，例如抽取样本的平均值。
首先，MCMC 办法抉择一个随机参数值来思考。模仿将持续生成随机值（这是蒙特卡洛局部），但要恪守一些规定来确定什么是好的参数值。窍门是，对于一对参数值，能够通过计算每个值解释数据的可能性来计算哪个是更好的参数值，给定咱们的先验信念。如果随机生成的参数值比上一个更好，则以肯定的概率将其增加到参数值链中，该概率取决于它的好坏水平（这是马尔可夫链局部）。
为了直观地解释这一点，让咱们回忆一下某个值的散布高度代表察看该值的概率。因而，咱们能够认为咱们的参数值（x 轴）展现了高概率和低概率的区域，显示在 y 轴上。对于单个参数，MCMC 办法从沿 x 轴随机采样开始：
红点是随机参数样本
因为随机样本受固定概率的影响，它们往往会在一段时间后收敛于咱们感兴趣的参数的最高概率区域：
蓝点仅代表任意工夫点之后的随机样本，此时预计会产生收敛。留神：垂直重叠点纯正是为了阐明目标。
收敛后，MCMC 采样会产生一组点，这些点是来自后验散布的样本。围绕这些点绘制直方图，并计算您喜爱的任何统计数据：
在 MCMC 模仿生成的样本集上计算的任何统计量都是咱们对实在后验散布统计量的最佳猜想。
MCMC 办法也可用于预计多个参数（例如人的身高和体重）的后验散布。对于 n 个参数，在 n 维空间中存在高概率区域，其中某些参数值集能够更好地解释察看到的数据。因而，我认为 MCMC 办法是在概率空间内随机抽样以近似后验散布。
点击题目查阅往期内容
R 语言实现 MCMC 中的 Metropolis–Hastings 算法与吉布斯采样
左右滑动查看更多
01
02
03
04
什么是 MCMC，什么时候应用它？
MCMC 只是一个从散布抽样的算法。
这只是泛滥算法之一。这个术语代表“马尔可夫链蒙特卡洛”，因为它是一种应用“马尔可夫链”（咱们将在前面探讨）的“蒙特卡罗”（即随机）办法。MCMC 只是蒙特卡洛办法的一种，只管能够将许多其余罕用办法看作是 MCMC 的简略特例。
我为什么要从散布中抽样？
从散布中抽取样本是解决一些问题的最简略的办法。
可能 MCMC 最罕用的办法是从贝叶斯推理中的某个模型的后验概率分布中抽取样本。通过这些样本，你能够问一些问题：“参数的平均值和可信度是多少？”。
如果这些样本 （ 查看文末理解数据获取形式 ） 是来自散布的独立样本，则 预计均值将会收敛在实在均值上。
假如咱们的指标散布是一个具备均值 m 和标准差的正态分布 s。
作为一个例子，思考用均值 m 和标准偏差 s 来预计正态分布的均值（在这里，我将应用对应于规范正态分布的参数）：
咱们能够很容易地应用这个 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这种事件真的很常见。在大多数贝叶斯推理中，后验散布是一些（可能很大的）参数向量的函数，您想对这些参数的子集进行推理。
在一个等级模型中，你可能会有大量的随机效应项被拟合，然而你最想对一个参数做出推论。在
贝叶斯框架中，您能够计算您感兴趣的参数在所有其余参数上的边际散布（这是咱们下面要做的）。
为什么“传统统计”不应用蒙特卡洛办法？
对于传统教学统计中的许多问题，不是从散布中抽样，能够使函数最大化或最大化。所以咱们须要一些函数来形容可能性并使其最大化（最大似然推理），或者一些计算平方和并使其最小化的函数。
然而，蒙特卡罗办法在贝叶斯统计中的作用与频率统计中的优化程序雷同，这只是执行推理的算法。所以，一旦你根本晓得 MCMC 正在做什么，你能够像大多数人把他们的优化程序当作黑匣子一样看待它，像一个黑匣子。
马尔可夫链蒙特卡罗
假如咱们想要抽取一些指标散布，然而咱们不能像从前那样抽取独立样本。有一个应用马尔科夫链蒙特卡洛（MCMC）来做这个的解决方案。首先，咱们必须定义一些事件，以便下一句话是有情理的：咱们要做的是试图结构一个马尔科夫链，它抽样的指标散布作为它的安稳散布。
定义
假如咱们有一个三态马尔科夫过程。让咱们 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 1P[i,j]给出了从状态 i 到状态的概率 j。
请留神，与行不同，列不肯定总和为 1：
colSums(P)
## [1] 0.95 0.60 1.45这个函数采纳一个状态向量 x（其中 x[i]是处于状态的概率 i），并通过将其与转移矩阵相乘来迭代它 P，使零碎后退到 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}从处于状态 1 的零碎开始（x 向量 [1,0,0] 也是如此，示意处于状态 1 的概率为 100%，不处于任何其余状态）
同样，对于另外两种可能的起始状态：
y2<-iterate.P(c(0,1,0),P,n)y3<-iterate.P(c(0,0,1),P,n)这表明了安稳散布的收敛性。
ma=1,xlab="Step",ylab="y",las=1)matlines(0:n,y2,lty=2)matlines(0:n,y3,lty=3)咱们能够应用 R 的 eigen 函数来提取零碎的次要特征向量（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)下面的过程迭代了不同状态的总体概率; 而不是通过零碎的理论转换。所以，让咱们迭代零碎，而不是概率向量。
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)绘制咱们在每个状态随工夫变动的工夫分数，而不是绘制状态：
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)所以这里的要害是：马尔可夫链有一些不错的属性。马尔可夫链有固定的散布，如果咱们运行它们足够长的工夫，咱们能够看看链条在哪里破费工夫，并对该安稳散布进行正当的预计。
Metropolis 算法
这是最简略的 MCMC 算法。
MCMC 采样 1d（单参数）问题
这是两个正态分布的加权和。这种散布相当简略，能够从 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 的正态分布中抽样
而这只须要运行 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当初，运行不同的计划 – 一个标准差很大（33），另一个标准差很小（3）。
res.fast<-run(-10action(x)rnorm(1,x,33),1000)res.slow<-run(-10,f,functanorm(1,x,.3),1000)留神三条轨迹正在挪动的不同形式。
相同，红色的痕迹回绝其中的大部分空间。
蓝色的形迹提出了偏向于被承受的小动作，然而它随着大部分的轨迹随机行走。它须要数百次迭代能力达到概率密度的大部分。
您能够在随后的参数中看到不同计划步骤在自相干中的成果 – 这些图显示了不同滞后步骤之间自相关系数的衰减，蓝线示意统计独立性。
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 在两个维度
给出了一个多元正态密度，给定一个均值向量（散布的核心）和方差 – 协方差矩阵。
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)从多元正态分布取样也相当简略，但咱们将应用 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<-s
uehist(samples[,1],freq=FALSE,ma,0.25))lines(xx,yy/z,col="red")数据获取
在上面公众号后盾回复“MCMC 数 据”，可获取残缺数据。
点击文末 “浏览原文”
获取全文残缺材料。
本文选自《R 语言中实现马尔可夫链蒙特卡罗 MCMC 模型》。
点击题目查阅往期内容
R 语言实现 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 采样用于回归的贝叶斯预计
