基于Markov Chain Monte Carlo的智能手表睡眠数据分析

为了帮助数友们提升数据科学实战能力以及加深对数据科学理论的认识水平，中国人民大学朝乐门老师团队策划并推出【Python数据科学实战系列】，为您全景详解数据科学领域的最佳实践。目前，已公布的课程有：

9.基于Markov Chain Monte Carlo的智能手表睡眠数据分析

8.按主题抓取多个网页及建立自己的数据集

7.Jupyter Notebook/Lab中添加R Kernel的详细步骤

6.Web信息爬取 | 详解 + Reddit等2个案例实践

5.基于MovieLens的影评趋势分析|详解

4.Windows和PC机上搭建Spark+Python开发环境的详细步骤

3.Jupyter Notebook/Lab中添加R Kernel的详细步骤

2.盘点数据科学领域常用的Python库

1.如何用Python学习数据科学

本文基于马尔可夫链蒙特卡罗（Markov Chain Monte Carlo简称MCMC）分析智能手表睡眠数据，为您解读基于Python与MCMC的数据科学实践。

1

研究问题及前期准备

智能手表可以记录下睡眠时间以及活动时间，如上图。本案例的目标是利用睡眠数据构建一个模型，模型指定睡眠后验概率作为时间函数，即返回在给定的时间进入睡眠状态的概率，数学表达式为P（sleeptime）。由于时间是一个连续变量，确定整个后验分布比较棘手，因此可以尝试研究近似分布，比如马尔可夫链蒙特卡罗（MCMC）。

1.1 数据准备

数据集包含两个月的活动和入睡信息，首先导入数据，并进行可视化。

import pandas as pd

import numpy as np

import scipy

from scipy import stats

import pymc3 as pm

import theano.tensor as tt

import scipy

from scipy import optimize

import matplotlib.pyplot as plt

%matplotlib inline

from IPython.core.pylabtools import figsize

import matplotlib

import json

s = json.load(open('../style/bmh_matplotlibrc.json'))

matplotlib.rcParams.update(s)

matplotlib.rcParams['figure.figsize'] = (10, 3)

matplotlib.rcParams['font.size'] = 14

matplotlib.rcParams['ytick.major.size'] = 20

# Markov Chain Monte Carlo设置初始样本数量

N_SAMPLES = 5000

# 读入数据

sleep_data = pd.read_csv('data/sleep_data.csv')

wake_data = pd.read_csv('data/wake_data.csv')

# 设置绘图坐标信息

sleep_labels = ['9:00', '9:30', '10:00', '10:30', '11:00', '11:30', '12:00']

wake_labels = ['5:00', '5:30', '6:00', '6:30', '7:00', '7:30', '8:00']

1.2入睡数据

首先对入睡数据进行可视化

print('Number of sleep observations %d' % len(sleep_data))

得到

Number of sleep observations 11340

figsize(16, 6)

# 入睡数据

plt.scatter(sleep_data['time_offset'], sleep_data['indicator'],

s= 60, alpha=0.01, facecolor = 'b', edgecolors='b')

plt.yticks([0, 1], ['Awake', 'Asleep']); plt.xlabel('PM Time');

plt.title('Falling Asleep Data', size = 18)

plt.xticks([-60, -30, 0, 30, 60, 90, 120], sleep_labels);

每个点表示特定时间的观察记录，颜色强度对应当时的点数。从图中可以看出，一般在晚上10点以后入睡。

1.3 活动数据

对活动数据进行可视化

# Wake data

plt.scatter(wake_data['time_offset'], wake_data['indicator'],s= 50, alpha = 0.01, facecolor='r', edgecolors = 'r');

plt.yticks([0, 1], ['Awake', 'Asleep']);

plt.xlabel('AM Time');

plt.title('Waking Up Data')

plt.xticks([-60, -30, 0, 30, 60, 90, 120], wake_labels);

从图中可以看出，活动数据的临界点在6点左右，可以判断出醒来时间在6点左右。

1.4 选择概率分布

在使用MCMC之前，需要确定一个适当的函数对睡眠的后验概率分布进行建模。有许多模型可以选择，首先可以尝试使用逻辑回归模型。 随着t→−∞,p(s|t)→0；t→+∞,p(s|t)→1。如果将睡眠的逻辑概率分布作为一个时间函数，那么表达式如下：

参数β未知，可以通过马尔可夫链蒙特卡罗采样来设定。MCMC从每个参数的先验样本中采样，尝试最大化给定数据的参数概率。以下展示了几个参数β值的逻辑回归函数：

figsize(16, 6)

def logistic(x, beta):

return 1. / (1. + np.exp(beta * x))

x = np.linspace(-5, 5, 1000)

for beta in [-5, -1, 0.5, 1, 5]:

plt.plot(x, logistic(x, beta), label = r"$\beta$ = %.1f" % beta)

plt.legend();

plt.title(r'Logistic Function with Different $\beta$ values');

像上面展示的基本逻辑回归模型，存在一个问题，转换的中心为0。然而在睡眠数据中，转换点为10:00pm和6:00am左右。因此需要添加偏差来调整逻辑函数的位置，最终表达式为：

因此又得到另外一个未知参数α，这也可以通过马尔可夫链蒙特卡罗得到。

对于参数α、β，设定不同的值，得到结果如下：

figsize(20, 8)

def logistic(x, beta, alpha=0):

return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha))

x = np.linspace(-4, 6, 1000)

plt.plot(x, logistic(x, beta=1), label=r"$\beta = 1, \alpha = 0$", ls="--", lw=2)

plt.plot(x, logistic(x, beta=-1), label=r"$\beta = -1 \alpha = 0$", ls="--", lw=2)

plt.plot(x, logistic(x, -1, 1),

label=r"$\beta = -1, \alpha = 1$", color="darkblue")

plt.plot(x, logistic(x, -1, -1),

label=r"$\beta = -1, \alpha = -1$",color="skyblue")

plt.plot(x, logistic(x, -2, 5),

label=r"$\beta = -2, \alpha = 5$", color="orangered")

plt.plot(x, logistic(x, -2, -5),

label=r"$\beta = -2, \alpha = -5$", color="darkred")

plt.legend(); plt.ylabel('Probability'); plt.xlabel('t')

plt.title(r'Logistic Function with Varying $\beta$ and $\alpha$');

从上图可以看出，参数β影响曲线方向和陡峭程度，参数α影响位置。下面将通过MCMC发现最可能的参数值。

2

Markov Chain Monte Carlo（MCMC）

马尔可夫链蒙特卡罗是一种概率分布抽样的方法，目的是构建最可能的分布模型。以本例来说，无法直接计算逻辑分布，所以可以为函数参数alpha和beta生成数千个值，称为样本，进而从中得到分布的近似值。MCMC背后的观点是，当生成的样本越多时，得到到近似值越接近真实分布。

马尔可夫链蒙特卡罗方法有两部分。蒙特卡罗指的是使用重复随机样本来获得数值答案的一般技术。蒙特卡罗可以认为是进行了许多实验，每次都改变模型中的变量并观察响应。通过选择随机值，可以探索大部分参数空间，即变量可能值的范围。

马尔可夫链是下一个状态仅依赖于当前状态的过程（在此上下文中的状态是指赋值给参数）。马尔可夫链是无记忆的，因为只有当前状态重要，而不是它到达当前状态的过程。如果这有点难以理解，请考虑日常现象即天气。如果我们想要预测明天的天气，我们可以仅使用今天的天气来获得合理的估计。如果今天下雪，我们会查看历史数据，显示下雪后第二天的天气分布情况，以估计明天天气的可能性。马尔可夫链的概念是我们不需要知道一个过程的整个历史来预测下一个输出，这种近似在许多现实情况中都能较好地工作。

综合马尔可夫链和蒙特卡罗的思想，MCMC是一种基于当前值重复绘制分布参数随机值的方法。每个值的样本都是随机的，但是值的选择受当前状态和假定的参数先验分布限制。 MCMC可以被认为是随机漫步，逐渐收敛到真正的分布。

为了绘制alpha和beta的随机值，需要为这些值假设一个先验分布。由于事先对参数没有任何假设，可以采用正态分布。

马尔可夫链蒙特卡罗将从两个正态分布中对β和α进行采样以找到参数。每次迭代，对β和α的估计都是从前面得出的。如果参数增加了数据的概率，那么该状态被接受，但如果参数与数据不一致，则状态被拒绝。蒙特卡罗指的是算法的采样部分。马尔可夫链意味着下一个状态仅依赖于一阶过程中的当前状态（二阶取决于当前和前一个步骤，当前和前两个步骤中的三阶等等）。 MCMC将返回指定步骤数量的每个参数样本，这被称为模型轨迹。为了找到最可能的参数，可以取迹线中样本的平均值。 MCMC没有给出确切的答案，而是试图找到基于数据的最大似然状态。

2.1 参数β和α的先验分布

没有证据表明模型参数参数β和α的先验分布是超前的。因此可以对他们建模，并假设来自于正态分布。正态分布，或者高斯分布，由均值μ和精度τ定义。精度是标准偏差σ的倒数。均值决定了分布的位置，精度显示分布情况。精度越高，表示数据越集中，因此变化越小。均值可以是正值或者负值，但精度是正值。此处定义的正态分布表示为：

下面显示了三种正态分布的概率密度函数。

figsize(20, 8)

nor = stats.norm

x= np.linspace(-10, 10, 1000)

mu = (-5, 0, 4)

tau = (0.5, 1, 2.5)

colors = ("forestgreen", "navy", "darkred")

params = zip(mu, tau, colors)

for param in params:

y = nor.pdf(x, loc = param[0], scale = 1 / param[1])

plt.plot(x, y,

label="$\mu = %d,\;\\tau = %.1f$" % (param[0], param[1]),

color = param[2])

plt.fill_between(x, y, color = param[2], alpha = 0.3)

plt.legend(prop={'size':18});

plt.xlabel("$x$")

plt.ylabel("Probability Density", size = 18)

plt.title("Normal Distributions", size = 20);

2.2 参数搜索空间

参数搜索空间，即变量可能值的范围。对于参数正态先验分布的参数空间可视化如下：

%matplotlib inline

import scipy.stats as stats

from IPython.core.pylabtools import figsize

import numpy as np

figsize(16, 8)

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

jet = plt.cm.jet

fig = plt.figure()

x = y = np.linspace(0, 5, 100)

X, Y = np.meshgrid(x, y)

plt.subplot(121)

norm_x = stats.norm.pdf(x, loc=0, scale=1)

norm_y = stats.norm.pdf(y, loc=0, scale=1)

M = np.dot(norm_x[:, None], norm_y[None, :])

im = plt.imshow(M, interpolation='none', origin='lower',

cmap=jet)

plt.xlim(0, 40)

plt.ylim(0, 40)

plt.title("Parameter Search Space for Normal Priors.")

ax = fig.add_subplot(122, projection='3d')

ax.plot_surface(X, Y, M, cmap=plt.cm.jet)

ax.view_init(azim=390)

plt.title("Parameter Search Space for Normal Priors");

正态分布的期望值就是均值。

正态分布的方差等同于：

再一次，对于参数β和α的先验分布中的μ或τ值，并没有任何假设。当初始建模时，可以使用μ=0以及相对较大的方差，比如τ=0.05。马尔可夫链蒙特卡罗将对μ和τ的值进行采样，试图使α和β的可能性最大化。

显然，并不能尝试参数空间中的每个点，而是通过对从较高概率地区（红色）随机抽样，进而创建最可能的模型。

2.3 具体算法

此处使用的MCMC具体算法为Metropolis–Hastings algorithm，译为梅特罗波利斯－黑斯廷斯算法（源于维基百科）。为了将观察数据连接到模型，每次绘制一组随机值时，算法会根据数据对他们进行评估。如果与数据不一致，参数值将会被拒绝，模型将保持在当前状态。如果随机值与数据一致，则将这些值分配给参数并成为当前状态。该过程继续进行指定的步骤数，模型的准确性随着步骤数量不断改进。

综合来说，在本案例中使用MCMC的基本过程如下：

为alpha和beta选择一组初始值，作为逻辑回归函数的参数。

根据当前状态随机地将新值分配给alpha和beta。

检查新的随机值是否与观察结果一致。若不一致，则拒绝这些值并返回到之前的状态。若一致，则接受这些值作为新的当前状态。

对指定的迭代次数重复步骤2和3。

该算法返回它为alpha和beta生成的所有值。然后，我们可以将这些值的平均值用作逻辑函数中alpha和beta的最可能的最终值。 MCMC无法返回“真正”的值，而是返回分布的近似值。给定数据下睡眠概率的最终模型将是具有α和β平均值的逻辑回归函数。

3

构建模型

基于以上准备工作，逻辑回归函数可以表达睡眠状态的转换，但是并不能确定模型参数alpha和beta的值。因此，构建模型的目标在于确定参数的值，以最大可能地揭示出观察数据。假定参数来自于由均值μ和方差τ定义的正态分布，MCMC算法将对参数alpha和beta采样均值μ和方差τ的值，以得到给定数据下逻辑回归函数最大可能性的参数值。数据通过伯努利变量连接到参数。

伯努利变量是一个离散的随机变量，可以时0或者1。在本案例中，可以将睡眠模型或者活动模型设置为一个伯努利变量，其中清醒为0，睡眠为1。睡眠数据的伯努利变量取决于时间，通过逻辑回归函数定义。

p(ti)是具有时间自变量的逻辑回归函数，因此以上公式变为：

MCMC的目标就是在使用数据集、假定正态先验分布的基础上找到参数alpha和beta的值。

以下构建三个模型：

睡眠模型，即入睡时间的模型；

活动模型，即醒来状态的模型；

睡眠时长模型，即睡眠时间长度的模型。

3.1 睡眠模型

本文通过Python中PyMC3库实现目的，这是一个功能强大的贝叶斯推断库，具有运算马尔可夫链蒙特卡罗和其他推断算法的功能。以下代码创建模型并执行MCMC，为β和α绘制N_SAMPLES个样本。具体的采样算法是Metropolic Hastings。我们提供数据并告诉模型它是伯努利变量的观测值。该模型在数据基础上尝试最大化参数。

# 根据时间偏移进行排序

sleep_data.sort_values('time_offset', inplace=True)

# 设置时间

time = np.array(sleep_data.loc[:, 'time_offset'])

# 设置时间观测点

sleep_obs = np.array(sleep_data.loc[:, 'indicator'])

with pm.Model() as sleep_model:

# 创建参数alpha和beta

alpha = pm.Normal('alpha', mu=0.0, tau=0.01, testval=0.0)

beta = pm.Normal('beta', mu=0.0, tau=0.01, testval=0.0)

# 根据逻辑函数创建概率

p = pm.Deterministic('p', 1. / (1. + tt.exp(beta * time + alpha)))

# 通过观测数据创建伯努利参数

observed = pm.Bernoulli('obs', p, observed=sleep_obs)

# 通过最大化后验估计设置初始值

# start = pm.find_MAP()

# 设置采样算法

step = pm.Metropolis()

# 进行采样

sleep_trace = pm.sample(N_SAMPLES, step=step, njobs=2);

注：采样过程需要花费大约3-5小时，具体时间由机器性能决定

变量轨迹包含从参数β和α后验估计中得到的所有样本，可以进一步进行可视化来探索他们在采样过程中的变化。MCMC的观点在于，随着算法的推进，给定数据下样本的可能性越高。也就说，当样本增加时，MCMC算法收敛于最可能的值。我们期望后面得到的值比前面的值更准确。在马尔可夫链蒙特卡罗中，通常做法是丢掉一部分样本，一般是50%，这些被称为老化样本。在本文中，没有丢弃任何样本，但在实际应用中，通常会多次运行该模型，而且丢弃初始样本。

以下对参数β和α后验估计中得到的所有样本进行可视化。

# 提取样本

alpha_samples = sleep_trace["alpha"][5000:, None]

beta_samples = sleep_trace["beta"][5000:, None]

figsize(16, 10)

plt.subplot(211)

plt.title(r"""Distribution of $\alpha$ with %d samples""" % N_SAMPLES)

plt.hist(alpha_samples, histtype='stepfilled', color = 'darkred', bins=30, alpha=0.8);

plt.ylabel('Probability Density')

plt.subplot(212)

plt.title(r"""Distribution of $\beta$ with %d samples""" % N_SAMPLES)

plt.hist(beta_samples, histtype='stepfilled', color = 'darkblue', bins=30, alpha=0.8)

plt.ylabel('Probability Density');

如果β值以0为中心，则表明时间对入睡概率没有影响。 α值也不为0，表明在睡眠时间从晚上10点开始有偏移。此处选择将时间表示为从下午10:00起的偏移量，以尽可能避免处理数据时间。

数据的传播导致测量时存在不确定性。由于睡眠和活动的观察值存在相当大的重叠，不确定性也会更大。为了找到最可能的睡眠后验分布，此处取取α和β样本的平均值。

# 设置概率预测的时间值

time_est = np.linspace(time.min()- 15, time.max() + 15, 1e3)[:, None]

# 取样本平均值

alpha_est = alpha_samples.mean()

beta_est = beta_samples.mean()

# 计算概率

sleep_est = logistic(time_est, beta=beta_est, alpha=alpha_est)

确定参数值后，进一步，对睡眠后验分布进行可视化。

figsize(16, 6)

plt.plot(time_est, sleep_est, color = 'navy',

lw=3, label="Most Likely Logistic Model")

plt.scatter(time, sleep_obs, edgecolor = 'slateblue',

s=50, alpha=0.2, label='obs')

plt.title('Probability Distribution for Sleep with %d Samples' % N_SAMPLES);

plt.legend(prop={'size':18})

plt.ylabel('Probability')

plt.xlabel('PM Time');

plt.xticks([-60, -30, 0, 30, 60, 90, 120], sleep_labels);

如上图所示，随着时间向后变化，后验概率从0增长到1。由于数据存在噪音，该模型并不完美，但它是基于观测数据的充分近似，可以提供有用的估计。

print('The probability of sleep increases to above 50% at 10:{} PM.'.format(int(time_est[np.where(sleep_est > 0.5)[0][0]][0])))

得到

The probability of sleep increases to above 50% at 10:14 PM.

进一步可视化

colors = ["#348ABD", "#A60628", "#7A68A6"]

cmap = matplotlib.colors.LinearSegmentedColormap.from_list("BMH", colors)

figsize(12, 6)

probs = sleep_trace['p']

plt.scatter(time, probs.mean(axis=0), cmap = cmap,

c = probs.mean(axis=0), s = 50);

plt.title('Probability of Sleep as Function of Time')

plt.xlabel('PM Time');

plt.ylabel('Probability');

plt.xticks([-60, -30, 0, 30, 60, 90, 120], sleep_labels);

从晚上10点开始的任何偏移点上，可以查询后验分布得到入睡的概率。

print('10:00 PM probability of being asleep: {:.2f}%.'.

format(100 * logistic(0, beta_est, alpha_est)))

print('9:30 PM probability of being asleep: {:.2f}%.'.

format(100 * logistic(-30, beta_est, alpha_est)))

print('10:30 PM probability of being asleep: {:.2f}%.'.

format(100 * logistic(30, beta_est, alpha_est)))

得到

10:00 PM probability of being asleep: 27.41%.

9:30 PM probability of being asleep: 4.80%.

10:30 PM probability of being asleep: 73.90%.

得到模型后，有多种模型诊断方法。比如，我们知道对α和β的估计存在相当大的不确定性。为了在图表中反映这一点，我们可以基于所有样本在每个时间得到95％的置信区间。

plt.fill_between(time_est[:, 0], *quantiles, alpha=0.6,

color='slateblue', label = '95% CI')

plt.plot(time_est, sleep_est, lw=2, ls='--',

color='black', label="average posterior \nprobability of sleep")

plt.xticks([-60, -30, 0, 30, 60, 90, 120], sleep_labels);

plt.scatter(time, sleep_obs, edgecolor = 'skyblue', s=50, alpha=0.1);

plt.legend(prop={'size':14})

plt.xlabel('PM Time'); plt.ylabel('Probability');

plt.title('Posterior Probabilty with 95% CI');

可以看出，几乎每个时间点上是否入睡都存在不确定性。这也表明MCMC得到的只是对事实数据的最近似估计，而无法得到真正的参数。

同时，也可以基于参数所有样本，对某个时间点上的后验分布绘制直方图，从而看一看模型中的不确定性。

def sleep_posterior(time_offset, time):

figsize(16, 8)

prob = logistic(time_offset, beta_samples, alpha_samples)

plt.hist(prob, bins=100, histtype='step', lw=4)

plt.title('Probability Distribution for Sleep at %s PM' % time)

plt.xlabel('Probability of Sleep'); plt.ylabel('Samples')

plt.show();

sleep_posterior(0, '10:00')

sleep_posterior(-30, '9:30')

sleep_posterior(30, '10:30')

print('Most likely alpha parameter: {:.6f}.'.format(alpha_est))

print('Most likely beta parameter: {:.6f}.'.format(beta_est))

得到

Most likely alpha parameter: 0.973708.

Most likely beta parameter: -0.067142.

我们如何知道模型是否收敛？可以看看轨迹，或取样过程中值的变化路径。另一种选择是查看样本的自相关性。在马尔可夫链建模过程中，样本与自身相关，因为下一个值取决于当前状态（或基于顺序的当前状态和过去状态）。最初，该算法倾向于在搜索空间中漫步，并且具有较高的自相关性。随着算法收敛，样本将稳定在一个值附近，并且收敛度是低自相关。本文中不对收敛进行严格的研究，但可以绘制所有样本的变化轨迹。

以下对参数α和β的所有样本进行可视化

figsize(12, 6)

# alpha的轨迹

plt.subplot(211)

plt.title(r'Trace of $\alpha$')

plt.plot(alpha_samples, color = 'darkred')

plt.xlabel('Samples'); plt.ylabel('Parameter');

# beta的轨迹

plt.subplot(212)

plt.title(r'Trace of $\beta$')

plt.plot(beta_samples, color='b')

plt.xlabel('Samples'); plt.ylabel('Parameter');

plt.tight_layout(h_pad=0.8)

此外，PyMC3库有许多用于模型评估的内置诊断函数。以下是轨迹图和自相关图。

figsize(20, 12)

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

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

3.2 活动模型

以下重复上述过程，对活动数据即清醒时候的数据进行建模。

首先进行采样

wake_data.sort_values('time_offset', inplace=True)

time = np.array(wake_data.loc[:, 'time_offset'])

wake_obs = np.array(wake_data.loc[:, 'indicator'])

with pm.Model() as wake_model:

alpha = pm.Normal('alpha', mu=0.0, tau=0.01, testval=0.0)

beta = pm.Normal('beta', mu=0.0, tau=0.01, testval=0.0)

p = pm.Deterministic('p', 1. / (1. + tt.exp(beta * time + alpha)))

observed = pm.Bernoulli('obs', p, observed=wake_obs)

step = pm.Metropolis()

wake_trace = pm.sample(N_SAMPLES, step=step, njobs=2);

注：采样过程需要花费大约3-5小时，具体时间由机器性能决定

进一步，绘制活动数据的后验分布

alpha_samples = wake_trace["alpha"][5000:, None]

beta_samples = wake_trace["beta"][5000:, None]

time_est = np.linspace(time.min()- 15, time.max() + 15, 1e3)[:, None]

alpha_est = alpha_samples.mean()

beta_est = beta_samples.mean()

wake_est = logistic(time_est, beta=beta_est, alpha=alpha_est)

figsize(16, 6)

plt.plot(time_est, wake_est, color = 'darkred', lw=3, label="average posterior \nprobability of wake")

plt.scatter(time, wake_obs, edgecolor = 'r', facecolor = 'r', s=50, alpha=0.05, label='obs')

plt.title('Posterior Probability of Wake with %d Samples' % N_SAMPLES);

plt.legend(prop={'size':14})

plt.ylabel('Probability')

plt.xlabel('AM Time');

plt.xticks([-60, -30, 0, 30, 60, 90, 120], wake_labels);

该模型没有陡峭的过渡，符合预期，因为研究对象倾向于早上6点左右醒来。几个小时后醒来多次，使得曲线向右移动。

print('The probability of being awake passes 50% at 6:{} AM.'.format(int(time_est[np.where(wake_est

得到

The probability of being awake passes 50% at 6:11 AM.

colors = ["#348ABD", "#A60628", "#7A68A6"]

cmap = matplotlib.colors.LinearSegmentedColormap.from_list("BMH", colors)

figsize(12, 6)

probs = wake_trace['p']

plt.scatter(time, probs.mean(axis=0), cmap = cmap, c = probs.mean(axis=0), s = 50);

plt.title('Probability of Sleep as Function of Time')

plt.xlabel('AM Time');

plt.ylabel('Probability');

plt.xticks([-60, -30, 0, 30, 60, 90, 120], wake_labels);

进一步查看某个时间点上活动的概率

print('Probability of being awake at 5:30 AM: {:.2f}%.'.

format(100 - (100 * logistic(-30, beta=beta_est, alpha=alpha_est))))

print('Probability of being awake at 6:00 AM: {:.2f}%.'.

format(100 - (100 * logistic(0, beta=beta_est, alpha=alpha_est))))

print('Probability of being awake at 6:30 AM: {:.2f}%.'.

format(100 - (100 * logistic(30, beta=beta_est, alpha=alpha_est))))

得到

Probability of being awake at 5:30 AM: 14.07%.

Probability of being awake at 6:00 AM: 37.89%.

Probability of being awake at 6:30 AM: 69.46%.

3.3 睡眠时长模型

基于给定数据，也可以得到一个模型来估计睡眠时长。首先查看数据，然后确定哪个函数符合概率分布。

raw_data = pd.read_csv('data/sleep_wake.csv')

raw_data['length'] = 8 - (raw_data['Sleep'] / 60) + (raw_data['Wake'] / 60)

duration = raw_data['length']

figsize(10, 8)

plt.hist(duration, bins = 20, color = 'darkred')

plt.xlabel('Hours'); plt.title('Length of Sleep Distribution');

plt.ylabel('Observations');

从图中可以看出，呈右偏分布。因此可以一个偏斜的分布对睡眠时长建模。也可以使用双模分布，因为似乎存在两种模式。现在，决定将睡眠时长分布表示为偏斜正态分布。偏斜的正态分布如下所示。

对模型参数进行采样

with pm.Model() as duration_model:

# 三个参数进行采样

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_ = pm.SkewNormal('duration', alpha = alpha_skew, mu = mu_,sd = 1/tau_, observed = duration)

# 进行采样

step = pm.Metropolis()

duration_trace = pm.sample(N_SAMPLES, step=step)

进一步，取平均值处理

# 从采样中提取最大可能性估计

alpha_skew_samples = duration_trace['alpha_skew'][5000:]

mu_samples = duration_trace['mu'][5000:]

tau_samples = duration_trace['tau'][5000:]

alpha_skew_est = alpha_skew_samples.mean()

mu_est = mu_samples.mean()

tau_est = tau_samples.mean()

对睡眠时长的后验分布进行可视化

x = np.linspace(6, 12, 1000)

y = stats.skewnorm.pdf(x, a = alpha_skew_est, loc=mu_est, scale=1/tau_est)

plt.plot(x, y, color = 'forestgreen')

plt.fill_between(x, y, color = 'forestgreen', alpha = 0.2);

plt.xlabel('Hours'); plt.ylabel('Probability');

plt.title('Posterior Distribution for Duration of Sleep');

plt.vlines(x = x[np.argmax(y)], ymin=0, ymax=y.max(),

linestyles='--', linewidth=2, color='red',

label = 'Most Likely Duration');

print('The most likely duration of sleep is {:.2f} hours.'.format(x[np.argmax(y)]))

查询某些时长的概率

print('Probability of at least 6.5 hours of sleep = {:.2f}%.'.format(100 * (1 - stats.skewnorm

.cdf(6.5, a = alpha_skew_est, loc = mu_est, scale = 1/tau_est))))

print('Probability of at least 8.0 hours of sleep = {:.2f}%.'.format(100 * (1 - stats.skewnorm

.cdf(8.0, a = alpha_skew_est, loc = mu_est, scale = 1/tau_est))))

print('Probability of at least 9.0 hours of sleep = {:.2f}%.'.format(100 * (1 - stats.skewnorm

.cdf(9.0, a = alpha_skew_est, loc = mu_est, scale = 1/tau_est))))

得到

Probability of at least 6.5 hours of sleep = 99.07%.

Probability of at least 8.0 hours of sleep = 44.12%.

Probability of at least 9.0 hours of sleep = 11.03%.

对原始数据分布及模型进行可视化

x = np.linspace(6, 12, 1000)

y = stats.skewnorm.pdf(x, a = alpha_skew_est, loc=mu_est, scale=1/tau_est)

figsize(10, 8)

# Plot the posterior distribution

plt.plot(x, y, color = 'forestgreen',

label = 'Model', lw = 3)

plt.fill_between(x, y, color = 'forestgreen', alpha = 0.2);

# Plot the observed values

plt.hist(duration, bins=10, color = 'red', alpha=0.8,

label='Observed', normed=True)

plt.xlabel('Hours'); plt.ylabel('Probability');

plt.title('Duration Model');

plt.vlines(x = x[np.argmax(y)], ymin=0, ymax=y.max(),

linestyles='--', linewidth=2, color='k',

label = 'Most Likely Duration');

plt.legend(prop={'size':12});

后验偏斜正态分布看起来对数据的拟合比较好。然而，从右边分布来看，数据拟合为两个分布似乎更好，第二种模式并不是以单个偏斜正态分布捕捉的。但总体而言，这个模型仍然为睡眠时长的概率提供了合理估计。

4

结 论

基于以上建模估计，可以得到一些结论：

平均入睡时间为10:14PM

平均醒来时间为6:11AM

平均睡眠时长为7.67个小时

通过以上模型可以得到关于研究对象睡眠模式的知识，更多的数据可以提高适用性。虽然本文提出了几个假设，并没有对模型进行全面的调查，但是用贝叶斯方法分析实际数据是一个很好的开始。项目，尤其使用实际的应用数据，是最好的学习方式。

END

本文原作者为 William Koehrsen，文章发布于 Towards Data Science，刘岩、朝乐门对其进行翻译、整理及实践后发布。

原文：https://towardsdatascience.com/markov-chain-monte-carlo-in-python-44f7e609be98