基於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"$eta$ = %.1f" % beta)
plt.legend();
plt.title(r"Logistic Function with Different $eta$ 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"$eta = 1, alpha = 0$", ls="--", lw=2)
plt.plot(x, logistic(x, beta=-1), label=r"$eta = -1 alpha = 0$", ls="--", lw=2)
plt.plot(x, logistic(x, -1, 1),
label=r"$eta = -1, alpha = 1$", color="darkblue")
plt.plot(x, logistic(x, -1, -1),
label=r"$eta = -1, alpha = -1$",color="skyblue")
plt.plot(x, logistic(x, -2, 5),
label=r"$eta = -2, alpha = 5$", color="orangered")
plt.plot(x, logistic(x, -2, -5),
label=r"$eta = -2, alpha = -5$", color="darkred")
plt.legend(); plt.ylabel("Probability"); plt.xlabel("t")
plt.title(r"Logistic Function with Varying $eta$ 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 $eta$ 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
probability 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 $eta$")
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
probability 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
TAG:全球大搜羅 |