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

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

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

这只是众多算法之一。这个术语代表“马尔可夫链蒙特卡洛”，因为它是一种使用“马尔可夫链”（我们将在后面讨论）的“蒙特卡罗”（即随机）方法。MCMC只是蒙特卡洛方法的一种，尽管可以将许多其他常用方法看作是MCMC的简单特例。

最近我们被客户要求撰写关于MCMC的研究报告，包括一些图形和统计输出。

相关视频：马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

，时长07:25

 相关视频：马尔可夫链蒙特卡罗方法MCMC原理与R语言实现|数据分享

，时长08:47

我为什么要从分布中抽样？

从分布中抽取样本是解决一些问题的最简单的方法。

可能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

我们可以通过在这种情况