用Python实现马尔可夫链蒙特卡罗

摘要： 本文经过用Python中的马尔可夫链蒙特卡罗实现了睡眠模型项目，并教会如何使用MCMC。

在过去的几个月里，我在数据科学领域里遇到一个术语：马尔可夫链蒙特卡罗（MCMC）。在博客或文章里，每次看到这个语，我都会摇摇头，有几回我试着学习MCMC和贝叶斯推理，但每次一开始，就很快放弃了。我学习新技术的方式都是把它应用到一个实际问题上。

经过使用一些数据和一本应用实战的书（Bayesian Methods for Hackers），我终于经过一个实际项目弄懂了MCMC。像往常同样，当把这些技术概念应用到实际问题中时，理解它们要比阅读书上的抽象概念更容易。本文经过介绍Python中的MCMC实现过程，最终教会了我使用这个强大的建模和分析工具。

本项目的完整代码和相关数据在GitHub上能够找到。本文重点讨论了应用程序和结果，涵盖了不少有深度的内容。

介绍

实际生活中的数据永远不是完美的，但咱们仍然能够经过正确的模型从噪音数据中提取有价值的信息。

本项目的目的是使用睡眠数据建立一个模型，该模型指定了睡眠的后验几率做为一个时间函数。因为时间是一个连续变量，由于指定整个后验分布是比较困难的，因此咱们转向接近于分布的方法-MCMC。

几率分布的选择

在使用MCMC以前，咱们须要肯定适当的函数来给睡眠的后验几率分布建模。最简单方法是经过可视化的方式进行数据检查，入睡的时间函数的观察结果以下图所示：

在图上，每一个数据节点被标记为一个点，点的强度代表了在指定时间的观测次数。个人表只记录入睡的时间，因此为了扩大数据量，我在时间点两边的每一分钟都加上了节点。若是手表记录我在晚上10：05分睡着了，那么以前的每分钟都被表示为0（清醒），以后的每分钟都表示为1（入睡）。这将会把大约60个晚上的观测扩展到11340个数据节点。

能够看到，我老是于晚上10:00以后入睡，咱们想要建立一个获取从醒到睡的几率转换的模型。当模型在一个精确的时间从清醒（0）到入睡（1）转换的时候，咱们能够给模型用一个简单的阶梯函数，由于我不会每晚都在同一时间入睡，因此咱们须要一个函数来把转换过程建模为渐变过程。给定数据的最佳选择是一个在0和1之间平稳转换的逻辑函数。如下是入睡几率做为时间函数的逻辑方程。

β（beta）和α（alpha）是咱们在MCMC过程当中必须学习的模型参数。以下所示具备不一样参数的逻辑函数：

逻辑函数适合这些数据，由于入睡转换的几率是逐渐的，同时获取了个人睡眠样本的变化。咱们但愿能为函数增长一个时间变量t，并找出入睡的几率，它必须在0和1之间。为了建立这个模型，咱们经过一个分类的技术做为MCMC，用数据来找到最佳的α和β参数。

马尔可夫链蒙特卡罗（MCMC）

MCMC是从几率分布中抽样以构造最可能分布的一类方法。咱们不能直接计算逻辑分布，因此咱们为函数的参数（α和β）生成数千个称为样本的数值，以建立分布的近似值。MCMC背后的思想是，随着生成更多的样本，近似值会愈来愈接近实际的分布。

MCMC方法有两个部分，蒙特卡罗是指使用重复随机样原本得到数值结果的通用技术。蒙特卡罗能够被认为是进行不少次的实验，每次改变模型中的变量并观察反应。经过选择随机值，咱们能够研究参数空间的一大部分，可能的变量值的范围。参数空间，以下图所示：

很显然，咱们不能尝试图中的每一个点，可是经过从较高几率（红色）的区域随机抽样，咱们能够建立最可能的模型。

马尔可夫链（Markov Chain）

马尔可夫链是当前状态被下一个状态所依赖的过程。马尔可夫链是不记录状态的，由于只有当前的状态才是重要的，而它如何到达那个状态并不重要。若是这有点难以理解，那么考虑一个平常现象，天气。若是咱们想预测明天的天气，只能根据今天的天气来进行一个合理的预测。若是今天下雪了，咱们查看一下下雪后次日的天气分布的历史数据，以预测明天是什么天气的几率。马尔可夫链的概念是，咱们不须要必定知道一个过程的整个历史数据来预测下一个输出，类似的现象在许多现实状况中都能很好地预测。

MCMC结合了马尔可夫链和蒙特卡罗的思想，是一种基于当前值对分布的参数重复抽取随机值的方法。每一个值的样本都是随机的，可是这些值的选择受到当前状态和参数的假定先验分布的限制。MCMC能够看做是逐渐收敛到实际分布的随机游动。

为了抽取α和β的随机值，咱们须要假设这些值的先验分布。因为咱们没有预先给参数假设，因此可使用一个正态分布。正态分布或高斯分布由平均值定义，代表了数据的位置、方差和分布范围。下图是具备不一样平均值和范围的几个正态分布:

咱们使用的特定MCMC算法称为Metropolis Hastings。为了将观测数据与模型联系起来，每次绘制一组随机值时，该算法都会根据数据进行评估。若是这些随机值与数据不一致，则会被拒绝，而且模型保持在当前状态。反之，若是与数据一致，则将这些随机值分配给参数并变成当前状态。这个过程持续必定数量的迭代，那么模型的精确度就会随着迭代数量的增长而提升。

综上所述，用MCMC解决问题的基本步骤以下:

（1）选择α和β以及逻辑函数的参数初始值集合；

（2）根据当前状态给α和β随机分配新的值；

（3）检查新的随机值是否符合观测值。若是不符合，则拒绝这些随机值，并返回到先前状态，反之，要是符合，则接受，并更新为新的当前状态；

（4）如需迭代，则重复步骤2和3；

算法返回其生成的全部α和β的值。而后，咱们可使用这些值的平均值做为逻辑函数中最有可能的α和β最终值。MCMC不能返回“True”值，而是返回一个分布的近似值。由已有数据获得的睡眠几率，其最终模型是具备α和β平均值的逻辑函数。

Python的实现

为了在Python中实现MCMC，咱们将使用贝叶斯推理库PyMC3，它对大部分细节进行了抽象，使咱们可以建立模型而不是空有理论。

下面的代码建立了具备参数α和β、几率p和观察值的完整模型。步骤变量引用特定的算法，而且变量sleep_trace保存了由模型生成的全部参数值。

withpm.Model() as sleep_model:

# Create the alpha and beta parameters
# Assume a normal distribution
alpha=pm.Normal('alpha', mu=0.0, tau=0.05, testval=0.0)
beta=pm.Normal('beta', mu=0.0, tau=0.05, testval=0.0)

# The sleep probability is modeled as a logistic function
p=pm.Deterministic('p', 1. / (1. +tt.exp(beta * time + alpha)))

# Create the bernoulli parameter which uses observed data to inform the algorithm
observed=pm.Bernoulli('obs', p, observed=sleep_obs)

# Using Metropolis Hastings Sampling
step=pm.Metropolis()

# Draw the specified number of samples
sleep_trace=pm.sample(N_SAMPLES, step=step);

为了更好地了解代码的运行，咱们能够看下全部的模型运行过程当中所产生的α和β值。

上图被称为轨迹图。咱们能够看到，每一个状态都与先前的马尔可夫链相关，但其值在蒙特卡罗抽样中振荡明显。

在MCMC中，很常见的作法是丢弃多达90%的轨迹。算法不当即收敛到真实分布，并且初始值每每也不许确。以后的参数值一般会更好一些，这意味着应该用它们来构建模型。咱们使用了10000个样本，并丢弃了前面的50%，可是在企业级应用中可能会使用成百上千个或数百万个样本。

MCMC收敛到真实值，但评估收敛多是很困难的。若是咱们想要最精确的结果，这是一个重要的因素。PyMC3库已经建立了用于评估模型质量的函数，包括轨迹图和自相关图。

pm.traceplot(sleep_trace, ['alpha', 'beta'])pm.autocorrplot(sleep_trace, ['alpha', 'beta'])

睡眠模型

在最终创建和运行模型以后，如今应该使用了。咱们将最后5000个α和β样本的平均值做为最有可能的参数值，这容许咱们建立单个曲线图来模拟指定时间以后的入睡几率:

该模型能很好地表示数据，还获取了我入睡模式固有的变化，它不是给出一个结果，而是给出一个几率。例如，能够查询该模型以找出我在指定时间入睡的几率，并找到几率超过50%的时间点:

晚上9点半入睡几率: 4.80%.

晚上9点半入睡几率: 27.44%.

晚上10点入睡几率: 73.91%.

在晚上10点14分，入睡几率提升到了50%以上。

能够看到我入睡的平均时间是晚上10点14分左右。

这些值是给定数据的最可能的估计值。然而，会有这些几率相关的不肯定性，由于模型是近似的。为了表示这种不肯定性，咱们能够在一个给定的时间点使用全部α和β的样原本进行入睡几率预测，而后根据结果绘制柱状：

上述结果给出了一个更好的MCMC模型能实际作到的指标。这个方法并无找到一个正确答案，而是可能值的一个样本。贝叶斯推理是有实际用处的，由于它以几率的形式表示预测的结果。咱们能够说获得一个最有可能的答案，可是更准确的说法应该是任何预测都是一系列值的范围。

睡醒模型

能够用我早上睡醒的相关数据来找到一个相似的模型。我一般在早上6点起床，下图根据观察结果显示了从睡觉到睡醒的最终模型。

能够经过模型来发现我在某一指定时间入睡的几率和最有可能睡醒的时间。

早上5点半睡醒的几率: 14.10%.

早上6点睡醒的几率: 37.94%.

早上6点半睡醒的几率: 69.49%.

在早上6点11分睡醒的几率超过50%。

入睡持续时间

最后一个我想建立的是入睡时间模型。首先，咱们须要找到一个函数来模拟数据的分布，但只能经过检查数据来找到。

一个普通的分布便可实现，但它不会获取右边的离群点。可使用两个相互独立的正态分布来表示这两种模式，可是，我会使用一个偏态分布。偏态分布有三个参数，平均值、方差、误差，而且这三个参数必须从MCMC算法中进行学习。下面的代码建立模型并实现Metropolis Hastings抽样：

withpm.Model() asduration_model:
# Three parameters to sample
alpha_skew=pm.Normal('alpha_skew', mu=0, tau=0.5, testval=3.0)
mu_ =pm.Normal('mu', mu=0, tau=0.5, testval=7.4)
tau_ =pm.Normal('tau', mu=0, tau=0.5, testval=1.0)

# Duration is a deterministic variable
duration_ =pm.SkewNormal('duration', alpha=alpha_skew, mu= mu_, 
sd=1/tau_, observed= duration)

# Metropolis Hastings for sampling
step=pm.Metropolis()
duration_trace=pm.sample(N_SAMPLES, step=step)

如今，咱们可使用三个参数的平均值来构造最有可能的分布。下图表示数据的最终偏态分布:

上图看起来很合乎实际状况，经过查询模型能够发现我至少得到必定的入睡持续时间，和最可能的入睡时间范围的几率：

入睡至少持续6.5小时的几率= 99.16%.

入睡至少持续8小时的几率= 44.53%.

入睡至少持续9小时的几率= 10.94%.

入睡最可能持续的时间是7.67小时

结论

经过完成这个项目，告诉咱们解决问题的重要性是最好解决实际生活中的问题。数据科学须要不断地进行知识积累，最有效的方法就是有了问题尽快开始。

---

本文做者：【方向】

阅读原文

本文为云栖社区原创内容，未经容许不得转载。