关于数据挖掘:R语言中实现马尔可夫链蒙特卡罗MCMC模型

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

原文出处：拓端数据部落公众号

什么是 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.025



a.true<-qnorm(p,m,s)



a.true



1## \[1\] -1.96

咱们能够通过在这种状况下的间接整合来预计这个

aion(x)

dnorm(x,m,s)

g<-function(a)

integrate(f,-Inf,a)$value

a.int<-uniroot(function(x)g(a10,0))$roota.int



1## \[1\] -1.96

并用 Monte Carlo 积分预计点：

a.mc<-unnasamples,p))

a.mc



## \[1\] -2.023



a.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 1

P[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<-5000

set.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 ## 187



1coda::effectiveSize(res.fast)



1 2## var1 ## 33.19



1coda::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

ue

hist(samples\[,1\],freq=FALSE,ma,0.25))

lines(xx,yy/z,col="red")

最受欢迎的见解

1.R 语言中的马尔科夫机制转换 (Markov regime switching) 模型模型 ”)

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

3.matlab 贝叶斯隐马尔可夫 hmm 模型

4. 用机器学习辨认一直变动的股市情况—隐马尔科夫模型 (HMM) 的利用的利用 ”)

5.R 语言隐马尔可夫模型 hmm 辨认一直变动的市场条件

6.matlab 中的隐马尔可夫模型 (HMM) 实现实现 ”)

7.matlab 实现 MCMC 的马尔可夫切换 ARMA – GARCH 模型

8.R 语言如何做马尔科夫转换模型 markov switching model

9.R 语言马尔可夫转换模型钻研交通伤亡人数