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

什么是 MCMC，什么时候应用它？

MCMC 只是一个从散布抽样的算法。

这只是泛滥算法之一。这个术语代表“马尔可夫链蒙特卡洛”，因为它是一种应用“马尔可夫链”（咱们将在前面探讨）的“蒙特卡罗”（即随机）办法。MCMC 只是蒙特卡洛办法的一种，只管能够将许多其余罕用办法看作是 MCMC 的简略特例。

正如下面的段落所示，这个话题有一个疏导问题，咱们会缓缓解决。

我为什么要从调配中抽样？

你可能没有意识到你想（实际上，你可能并不想）。然而，从散布中抽取样本是解决一些问题的最简略的办法。

可能 MCMC 最罕用的办法是从贝叶斯推理中的某个模型的后验概率分布中抽取样本。通过这些样本，你能够问一些问题：“参数的平均值和可信度是多少？”。

如果这些样本是来自散布的独立样本，则 预计均值将会收敛在实在均值上。

假如咱们的指标散布是一个具备均值 m 和标准差的正态分布 s。显然，这种散布的意思是 m，但咱们试图通过从散布中抽取样本来展现。

作为一个例子，思考用均值 m 和标准偏差 s 来预计正态分布的均值（在这里，我将应用对应于规范正态分布的参数）：

咱们能够很容易地应用这个 rnorm 函数从这个散布中抽样

 seasamples<-rn 000,m,s)

样本的平均值十分靠近实在平均值（零）：

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

事实上，在这种状况下，$ n $ 样本预计的预期方差是 $ 1 / n $，所以咱们预计大部分值在 $ \\ pm 2 \\，/ \\ sqrt {n} = 0.02 $ 10000 分的实在意思。

 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

这个函数计算累积平均值（即元素 $ k $，元素 $ 1,2，\\ ldots，k $ 除以 $ k $）之和。

 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 积分预计点：

1 2a.mc<-unnasamples,p))a.mc1## [1] -2.0231a.true-a.mc1## [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

这种事件真的很常见。在大多数贝叶斯推理中，后验散布是一些（可能很大的）参数向量的函数，您想对这些参数的子集进行推理。在一个等级模型中，你可能会有大量的随机效应项被拟合，然而你最想对一个参数做出推论。在贝叶斯框架中，您能够计算您感兴趣的参数在所有其余参数上的边际散布（这是咱们下面要做的）。

为了阐明这个问题，

思考在边长为 $ 2r $ 的方格内半径为 $ r $ 的圆; 空间的“乏味”区域是一个随机抉择的点位于圆圈内的一个很好的机会。

对于半径为 $ 2r $ 的立方体中半径为 $ r $ 的球体，球体的体积为 $ 4 /（3 \\ pi r ^ 3）$，立方体的体积为 $（2d）^ 3 $，作为问题的维数，$ d $ 减少（应用超立方体中的超球面）

 d<-2:10plot(d,pi^(d/2)/(d*2^(d-1)*gamma(d/2)),logahere")

所以咱们不须要减少很多维度来次要对潜在空间的一小部分感兴趣。

 a<-function(d,n)mean(repa00)) ## dimension p.interesting## 1 1 0.5219## 2 2 0.2173## 3 3 0.0739## 4 4 0.0218## 5 5 0.0070## 6 6 0.0025## 7 7 0.0006## 8 8 0.0000## 9 9 0.0000## 10 10 0.0000

即便只看 4 - 5 个维度，如果咱们试图对参数空间进行彻底整合

为什么“失常统计”不应用蒙特卡洛办法？

对于传统教学统计中的许多问题，而不是从散布中抽样，能够使函数最大化或最大化。所以咱们须要一些函数来形容可能性并使其最大化（最大似然推理），或者一些计算平方和并使其最小化的函数。

然而，蒙特卡罗办法在贝叶斯统计中的作用与频率统计中的优化程序雷同，这只是执行推理的算法。所以，一旦你根本晓得 MCMC 正在做什么，你能够像大多数人把他们的优化程序当作黑匣子一样看待它，像一个黑匣子。

马尔可夫链蒙特卡罗

假如咱们想要抽取一些指标散布，然而咱们不能像从前那样抽取独立样本。有一个应用马尔科夫链蒙特卡洛（MCMC）来做这个的解决方案。首先，咱们必须定义一些事件，以便下一句话是有情理的：咱们要做的是试图结构一个马尔科夫链，它的难以抽样的指标散布作为它的安稳散布。

定义

让 $ X\_t $ 示意在工夫 $ t $ 时的一些随机变量的值。马尔可夫链从某个点 $ X\_0 $ 开始，生成一系列样本 $ \[X0，X1，X2，\\ ldots，Xt\] $，而后遵循一系列随机步骤。

 假如咱们有一个三态马尔科夫过程。让咱们 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] 给出了从状态 i 到状态的概率 j（所以这是下面提到的 $ P（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,lena<-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)# normalise eigenvector

而后在之前的数字上加上一点，表明咱们有多靠近收敛：

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 采样器，而只是为了看到他们的确在工作。

算法进行如下。

从 $ x\_t $ 开始。

倡议一个新的状态 $ x ^ \\ prime $

计算“承受概率”

从 $ \[0,1\] $ 中抽取一些均匀分布的随机数 $ u $; 如果 $ u <\\ alpha $ 承受该点，则设置 $ x {t + 1} = x ^ \\ prime $。否则回绝它并设置 $ x {t + 1} = x\_t $。

请留神，在下面的步骤 3 中，未知归一化常数因为而退出

这将产生一系列样本 $ {x 0，x 1，\\ ldots} $。请留神，如果倡议的样本被回绝，雷同的值将呈现在间断的样本中。

还要留神，这些不是来自指标散布的独立样本; 他们是依赖样本 ; 也就是说，示例 $ x\_t $ 取决于 $ x\_ {t-1} $ 等等。然而，因为链条靠近安稳散布，只有咱们抽取足够的点数，这种依赖性就不会有问题。

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

咱们抉择一个中央开始（如何 -10，只是抉择一个十分蹩脚的一点）

这里是马尔可夫链的前 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="x",ylabaadd=TRUE,col="red",n=200)

当初，运行不同的提案机制 – 一个标准差很大（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,main="Large steps")

由此能够计算独立样本的无效数量：

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 步：

 par(mfrow=c(2,2),mar=rep(.5,4),oma=c(4,4,0,0))for(hinhh){plot(h,main="",freq=a=300)}

MCMC 在两个维度

给出了一个多元正态密度，给定一个均值向量（散布的核心）和方差 – 协方差矩阵。

 make.mvn<-function(mean,vcv){logdet<-as.numeric(detea+logdetvcv.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<-sapply(xx,m)z<-integrate(splinefun(xx,yy),min(xx),max(xx))$valuehist(samples[,1],freq=FALSE,ma,0.25))lines(xx,yy/z,col="red")

相干文章

R 语言用 Backfitting _MCMC_抽样算法进行贝叶斯推理案例

 R 语言中实现马尔可夫链蒙特卡罗_MCMC_模型 

R 语言中的 Stan 概率编程_MCMC_采样的贝叶斯模型

R 语言实现_MCMC_中的 Metropolis–Hastings 算法与吉布斯 …

matlab 实现_MCMC_的马尔可夫切换 ARMA – GARCH 模型预计