过去的几个月里，我在数据科学领域反复听到一个术语：马尔可夫链蒙特卡罗（MCMC）。工作中每次听到这个词时，我都会点头并觉得这听起来很酷，但实际上对其含义只有模糊的理解。我多次尝试学习MCMC和贝叶斯推理，但每次开始阅读相关书籍时，我很快就放弃了。感到沮丧后，我转向了学习任何新技能的最佳方法：将其应用到一个实际问题中。
利用一些我一直想探索的睡眠数据和一本基于动手实践的书籍（《Bayesian Methods for Hackers》），我最终通过一个现实世界的项目学习了马尔可夫链蒙特卡罗。通常，当我将这些技术概念应用到实际问题时，比单纯阅读书本上的抽象概念要容易得多（而且更有趣）。这篇文章将带你了解如何在Python中实现马尔可夫链蒙特卡罗，这也是我最终掌握这一强大建模和分析工具的方法。
该项目的完整代码和数据已上传到GitHub。我鼓励任何有兴趣的人查看并在自己的数据上使用它。
项目地址：https://github.com/WillKoehrsen/ai-projects/tree/master/markov_chain_monte_carlo
介绍
我的Garmin Vivosmart手表通过心率和动作追踪我的入睡和醒来时间。尽管它并不是100%准确，但现实世界中的数据从来都不完美，然而，我们依然可以通过合适的模型从噪声数据中提取出有用的信息！
这个项目的目标是利用睡眠数据创建一个模型，该模型能够根据时间来指定睡眠的后验概率。由于时间是一个连续变量，直接指定整个后验分布是不可行的，因此我们转向了近似分布的方法，比如马尔可夫链蒙特卡罗（MCMC）。
概率分布选择
在开始使用MCMC之前，我们需要确定一个合适的函数来建模睡眠的后验概率分布。一个简单的方法是通过视觉检查数据。下面显示了我入睡时间作为时间函数的观测数据。
每个数据点都以一个点表示，点的强度显示在特定时间的观测次数。我的手表只记录了我入睡的具体分钟时间，因此为了扩展数据，我在精确时间的两侧每一分钟都添加了数据点。如果我的手表显示我在晚上10:05入睡，那么在这之前的每一分钟都标记为0（清醒），而之后的每一分钟都标记为1（入睡）。通过这种方式，我将大约60个晚上的观测数据扩展到了11340个数据点。
从图中可以看到，我通常在晚上10:00之后不久入睡，但我们希望创建一个模型，以概率的形式捕捉从清醒到入睡的过渡。我们可以使用一个简单的阶跃函数作为模型，该函数在某个精确时间从清醒（0）变为入睡（1），但这无法体现数据中的不确定性。我并不是每晚都在同一时间入睡，因此我们需要一个函数来将这个过渡建模为一个渐进的过程，以显示其中的可变性。鉴于数据的特点，最好的选择是使用一个逻辑函数（Logistic Function），它可以平滑地从0过渡到1。以下是一个用于描述睡眠概率随时间变化的逻辑方程。
在这里，β（贝塔）和α（阿尔法）是模型的参数，我们需要在MCMC过程中学习这些参数。下面显示的是一个具有不同参数的逻辑函数。
逻辑函数适合这组数据，因为入睡的概率是逐渐过渡的，这很好地捕捉了我睡眠模式的可变性。我们希望能够将时间 t 输入到函数中，并输出睡眠的概率，这个概率必须在0和1之间。这样一来，我们可以获得一个概率值，而不是一个简单的“我是否在晚上10:00入睡”的两种答案。为了创建这个模型，我们使用数据通过马尔可夫链蒙特卡罗（MCMC）技术之一来找到最佳的 α 和 β 参数。
马尔科夫蒙特卡罗
马尔可夫链蒙特卡罗（MCMC）指的是一类用于从概率分布中抽样以构建最可能分布的方法。由于我们无法直接计算逻辑分布，因此我们生成数千个函数参数（α 和 β）的值，称为样本，以近似该分布。MCMC的思想是，随着我们生成更多的样本，我们的近似值会越来越接近实际的真实分布。
马尔可夫链蒙特卡罗方法有两个部分。蒙特卡罗（Monte Carlo）指的是一种通过重复随机抽样来获得数值解的一般技术。蒙特卡罗可以被认为是进行许多实验，每次改变模型中的变量并观察响应。通过选择随机值，我们可以探索参数空间的大部分区域，即变量的可能值范围。下图展示了在我们的问题中使用正态先验（关于这一点稍后会详细介绍）的参数空间。
显然，我们无法尝试这些图中的每一个点，但通过从更高概率区域（红色区域）随机抽样，我们可以为我们的问题创建出最可能的模型。
马尔科夫链
马尔可夫链是一个过程，其中下一个状态仅依赖于当前状态。（在这个上下文中，状态指的是参数的值分配）。马尔可夫链是无记忆的，因为只有当前状态是重要的，而与如何到达该状态无关。如果这有点难以理解，可以考虑一个日常现象：天气。如果我们想预测明天的天气，只需根据今天的天气就能得到一个合理的估计。如果今天下雪了，我们可以查看历史数据，了解雪后第二天天气的分布，以估计明天的天气概率。马尔可夫链的概念是，我们不需要知道整个过程的历史才能预测下一个输出，这种近似在许多实际情况下都很有效。
将马尔可夫链和蒙特卡罗结合起来，MCMC是一种方法，通过当前值反复抽取分布参数的随机值。每个样本值都是随机的，但这些值的选择受到当前状态和参数的先验分布的限制。MCMC可以看作是一种随机游走，它逐渐收敛到真实的分布。
为了抽取α和β的随机值，我们需要假设这些值的先验分布。由于我们对参数没有事先假设，我们可以使用正态分布。正态分布（或高斯分布）由均值和方差定义，均值表示数据的位置，方差表示数据的分散程度。以下是几个均值和方差不同的正态分布图。
我们使用的特定MCMC算法称为Metropolis Hastings。为了将我们的观测数据与模型连接起来，每次抽取一组随机值时，算法都会将这些值与数据进行比较。如果这些值与数据不符（这里我稍作简化），则这些值会被拒绝，模型保持在当前状态。如果随机值与数据一致，则这些值会被分配给参数，并成为当前状态。这个过程会继续进行指定数量的步骤，随着步骤的增加，模型的准确性会逐渐提高。
将所有步骤综合起来，我们在这个问题中的马尔可夫链蒙特卡罗基本程序如下：
选择逻辑函数参数α和β的初始值。 根据当前状态随机分配新的α和β值。 检查新的随机值是否与观测数据一致。如果不一致，拒绝这些值并返回到之前的状态。如果一致，接受这些值作为新的当前状态。 对指定的迭代次数重复步骤2和步骤3。
算法返回所有生成的α和β值。我们可以使用这些值的平均值作为逻辑函数中α和β的最可能的最终值。MCMC不能返回“真实”值，而是提供对分布的近似。根据数据得到的最终睡眠概率模型将是具有α和β平均值的逻辑函数。
Python实现
上述细节我多次听不太懂，直到我在Python中实际应用它们！亲眼看到结果比阅读他人描述要有帮助得多。为了在Python中实现MCMC，我们将使用PyMC3贝叶斯推理库。它抽象化了大部分细节，使我们可以在不迷失于理论中的情况下创建模型。
以下代码创建了一个完整的模型，包括参数α和β、概率p和观测数据observed。step变量指代具体的算法，而sleep_trace保存了模型生成的所有参数值。
with pm.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);
为了了解运行这段代码时发生了什么，我们可以查看模型运行期间生成的所有α和β值。
这些图被称为轨迹图（trace plots）。我们可以看到，每个状态与前一个状态相关联——这就是马尔可夫链（Markov Chain）——但这些值显著波动——这就是蒙特卡罗抽样（Monte Carlo sampling）。
在MCMC中，通常会丢弃多达90%的轨迹。算法不会立即收敛到真实分布，初始值往往不准确。参数的后期值通常更好，这意味着它们是我们应该用于构建模型的值。我们使用了10000个样本，并丢弃了前50%的样本，但在实际应用中，可能会使用成千上万甚至数百万个样本。
MCMC在足够多的步骤下会收敛到真实值，但评估收敛性可能会比较困难。我不会在这篇文章中深入探讨这个话题（其中一种方法是测量轨迹的自相关性），但如果我们想要最准确的结果，这是一个重要的考虑因素。PyMC3内置了评估模型质量的函数，包括轨迹图和自相关图。
pm.traceplot(sleep_trace, ['alpha', 'beta'])
pm.autocorrplot(sleep_trace, ['alpha', 'beta'])
入睡模型
在最终构建并运行模型之后，现在是时候使用结果了。我们将使用最后5000个α和β样本的平均值作为参数的最可能值，这使我们能够创建一条单一的曲线来建模后验睡眠概率。
该模型很好地表示了数据。此外，它还捕捉到了我睡眠模式中的内在变异性。与直接给出一个“是”或“否”的答案不同，模型提供了一个概率。例如，我们可以查询模型，找出我在某个特定时间入睡的概率，并找出睡眠概率超过50%的时间点。
9:30 PM probability of being asleep: 4.80%.
10:00 PM probability of being asleep: 27.44%.
10:30 PM probability of being asleep: 73.91%.
The probability of sleep increases to above 50% at 10:14 PM.
虽然我尽量在晚上10点睡觉，但很显然大多数晚上并没有做到！我们可以看到，我的平均入睡时间大约是晚上10点14分。
这些值是根据数据得出的最可能估计值。然而，由于模型是近似的，所以这些概率存在不确定性。为了表示这种不确定性，我们可以使用所有的 alpha 和 beta 样本来预测某个时间点的睡眠概率，而不是使用平均值，然后绘制结果的直方图。
这些结果更好地说明了 MCMC 模型的真正作用。这种方法并不是找到单一的答案，而是得到一组可能值的样本。贝叶斯推断在现实世界中非常有用，因为它以概率的形式表达预测。我们可以说有一个最可能的答案，但更准确的说法是，任何预测都有一系列可能的值。
清醒模型
我可以利用清醒数据来建立一个类似的模型，用于预测我早上醒来的时间。我尽量总是在早上6:00被闹钟叫醒，但我们可以看到，这并不总是发生！下图展示了从睡眠到清醒的转变的最终模型以及相关的观测结果。
我们可以查询模型，以找出我在特定时间处于睡眠状态的概率，以及我最有可能醒来的时间。
Probability of being awake at 5:30 AM: 14.10%.
Probability of being awake at 6:00 AM: 37.94%.
Probability of being awake at 6:30 AM: 69.49%.
The probability of being awake passes 50% at 6:11 AM.
看起来我需要对那个闹钟做一些调整！
睡眠持续时间
我想创建的最后一个模型——既出于好奇，也为了练习——是我的睡眠时长。首先，我们需要找到一个函数来建模数据的分布。事先我认为这可能是正态分布，但我们只能通过检查数据来确定！
正态分布可以有效，但无法捕捉右侧的离群点（例如严重睡过头的情况）。我们可以使用两个独立的正态分布来表示两个模式，但我选择使用偏态正态分布。偏态正态分布有三个参数：均值、方差和偏度 alpha。这三个参数都必须通过 MCMC 算法来学习。以下代码创建了模型并实现了 Metropolis Hastings 采样。
with pm.Model() as duration_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)
现在，我们可以使用这三个参数的平均值来构建最可能的分布。以下是最终的偏态正态分布，叠加在数据上。
这看起来非常契合！我们可以查询模型，以找出我获得至少一定睡眠时间的可能性，以及最可能的睡眠时长。
Probability of at least 6.5 hours of sleep = 99.16%.
Probability of at least 8.0 hours of sleep = 44.53%.
Probability of at least 9.0 hours of sleep = 10.94%.
The most likely duration of sleep is 7.67 hours.
结论
再次完成这个项目让我体会到解决问题的重要性，特别是那些具有现实应用的问题！在构建基于马尔可夫链蒙特卡罗（MCMC）的贝叶斯推断的完整实现过程中，我掌握了许多基本知识，并在这个过程中获得了乐趣。数据科学就是不断为你的工具箱添加新工具，而最有效的方法就是找到一个问题并开始动手！
欢迎扫码关注：
