通過把馬爾科夫鏈蒙特卡羅（MCMC）應(yīng)用于一個(gè)具體問題，本文介紹了 Python 中 MCMC 的入門級應(yīng)用。機(jī)器之心對本文進(jìn)行了編譯介紹。

過去幾月中，我總是反復(fù)遇到同一個(gè)數(shù)據(jù)科學(xué)術(shù)語：馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo/MCMC）。每當(dāng)我在實(shí)驗(yàn)室、博客、文章中聽到這個(gè)概念，我常常點(diǎn)頭贊同，覺得它很酷，但實(shí)際上并沒有一個(gè)清晰的認(rèn)知。有幾次我嘗試著學(xué)習(xí) MCMC 和貝葉斯推理，我每次從閱讀書籍開始，結(jié)果卻很快放棄。我感到很惱怒，于是決定轉(zhuǎn)向一種學(xué)習(xí)任何新技能的最佳方法：將它應(yīng)用于一個(gè)具體問題。

使用我的睡眠數(shù)據(jù)（我一直打算對此一探究竟）和一本實(shí)際應(yīng)用的書（Bayesian Methods for Hackers），我終于通過一個(gè)實(shí)際問題學(xué)習(xí)了馬爾科夫鏈蒙特卡羅。像往常一樣，比起閱讀抽象的概念，將這些技術(shù)應(yīng)用到具體問題中能讓學(xué)習(xí)變得更簡單、更愉快。本文介紹了 Python 中的馬爾科夫鏈蒙特卡羅的入門級應(yīng)用，正是它教會了我使用這個(gè)強(qiáng)大的建模分析工具。

我鼓勵(lì)大家參閱 GitHub，并將其用于自己的數(shù)據(jù)當(dāng)中。本文將重點(diǎn)介紹它的應(yīng)用和結(jié)果，所以會產(chǎn)生很多高層次的話題。如果你在閱讀后想了解更多，可以閱讀文中提供的鏈接。

簡介

我的 Garmin Vivosmart 手表可以根據(jù)心率和運(yùn)動(dòng)情況追蹤我的睡眠和起床狀況。它并非 100％準(zhǔn)確，不過真實(shí)數(shù)據(jù)從不完美，我們?nèi)匀豢梢越柚_的模型從噪聲數(shù)據(jù)中提取有用的知識！

典型的睡眠數(shù)據(jù)

本項(xiàng)目的目標(biāo)是借助睡眠數(shù)據(jù)創(chuàng)建一個(gè)模型，通過把睡眠看作時(shí)間函數(shù)，而確定睡眠的后驗(yàn)概率。由于時(shí)間是連續(xù)變量，確定整個(gè)后驗(yàn)分布非常棘手。因此我們轉(zhuǎn)而使用一些可實(shí)現(xiàn)近似分布的方法，比如馬爾可夫鏈蒙特卡羅（MCMC）。

選擇一個(gè)概率分布

在開始使用 MCMC 之前，我們需要確定一個(gè)合適的函數(shù)來對睡眠的后驗(yàn)概率分布進(jìn)行建模。一個(gè)簡單的方法是直觀檢查這些數(shù)據(jù)。對于我的睡眠的時(shí)間函數(shù)的觀察如下圖所示。

睡眠數(shù)據(jù)

上圖中，每個(gè)數(shù)據(jù)點(diǎn)都用點(diǎn)表示，點(diǎn)的強(qiáng)度顯示在特定時(shí)間的觀測數(shù)量。我的手表只能記錄我入睡的那一分鐘，所以為了擴(kuò)大數(shù)據(jù)量，我在精確時(shí)間的兩邊增加以分鐘為單位的數(shù)據(jù)點(diǎn)。舉例而言，如果我的手表顯示我在晚上 10:05 入睡，那么 10:05 之前的每一分鐘都被表示為 0（清醒），10:05 之后的每一分鐘都被表示為 1（睡著）。這將大約 60 個(gè)夜晚的觀測數(shù)據(jù)擴(kuò)展到了 11340 個(gè)數(shù)據(jù)點(diǎn)。

我們可以發(fā)現(xiàn)，我一般在晚上 10 點(diǎn)之后入睡。但是我們想創(chuàng)建一個(gè)模型，以概率的形式捕捉從清醒到入睡的過渡過程。我們可以在模型中使用一個(gè)簡單的階躍函數(shù)，它在一個(gè)精確的時(shí)間從喚醒（0）過渡成入睡（1），但是這無法表現(xiàn)數(shù)據(jù)的不確定性。我不可能在每天晚上的同一時(shí)間睡覺，因此需要一個(gè)能模擬過渡過程的函數(shù)對這一漸進(jìn)過程進(jìn)行建模，顯示變化特性。給定上述數(shù)據(jù)的情況下，我們的最佳選擇是在 0 和 1 的邊界之間平滑過渡的 logistic 函數(shù)。以下是睡眠概率作為時(shí)間函數(shù)的 logistic 方程：

其中，β 和 α 是我們在 MCMC 過程中必須學(xué)習(xí)的模型參數(shù)。具有不同參數(shù)的 logsitic 函數(shù)圖像如下所示。

logsitic 函數(shù)很適合本案例中的數(shù)據(jù)，因?yàn)槿胨目赡苄詴饾u轉(zhuǎn)變，此函數(shù)能捕捉睡眠模式之中的變化情況。我們希望能夠在函數(shù)中插入時(shí)間 t，獲得睡眠概率（其值在 0 和 1 之間）。我們最終得到的不是在晚上 10:00 入睡與否的直接答案，而是一個(gè)概率。為了建立這個(gè)模型，我們使用這些數(shù)據(jù)，通過 MCMC 尋找最佳的 α 和 β 參數(shù)。

馬爾科夫鏈蒙特卡羅

馬爾可夫鏈蒙特卡羅指從概率分布中抽樣以構(gòu)建最大可能分布的一類方法。我們不能直接構(gòu)建 logistic 分布，所以，與之相反，我們?yōu)楹瘮?shù)的參數(shù)（α 和 β）生成了上千個(gè)值——被稱為樣本——從而創(chuàng)造分布的近似值。MCMC 背后的思想是，當(dāng)我們生成更多的樣本時(shí)，我們的近似值越來越接近實(shí)際的真實(shí)分布。

馬爾科夫鏈蒙特卡羅方法分為兩部分。蒙特卡羅指的是使用重復(fù)隨機(jī)樣本獲得數(shù)值解的一般性技術(shù)。蒙特卡羅可以被視為進(jìn)行了若干次實(shí)驗(yàn)，其中每次都對模型中的變量進(jìn)行改變并觀察其響應(yīng)。通過選擇隨機(jī)數(shù)，我們可以探索大部分參數(shù)空間，即變量可能值的范圍。下圖顯示了我們的問題使用正常先驗(yàn)后的參數(shù)空間。

顯然，我們無法一一嘗試圖像中的每一個(gè)點(diǎn)。但是通過對較高概率區(qū)域（紅色區(qū)域）進(jìn)行隨機(jī)抽樣，我們可以為問題建立最可能的模型。

馬爾科夫鏈

馬爾科夫鏈?zhǔn)且粋€(gè)隨機(jī)過程，其中次態(tài)僅依賴于當(dāng)前狀態(tài)（在此語境中，一個(gè)狀態(tài)指的是參數(shù)的一次賦值）。馬爾科夫鏈沒有記憶性，因?yàn)橹挥挟?dāng)前狀態(tài)對下一狀態(tài)起作用，而與到達(dá)當(dāng)前狀態(tài)的方式無關(guān)。如果這種說法還是有些難以理解，我們可以用日?，F(xiàn)象中的天氣來舉例。如果我們想預(yù)測明日天氣，我們可以僅通過今日天氣來得到一個(gè)合理的估計(jì)。如果今天下雪了，我們可以查看下雪次日天氣分布的歷史數(shù)據(jù)，估算明天天氣的概率。馬爾科夫鏈的概念在于，我們無需了解整個(gè)歷史過程就能預(yù)測下一狀態(tài)，這個(gè)近似在許多現(xiàn)實(shí)情況中就能很好地工作。

綜合馬爾科夫鏈和蒙特卡羅的思想，馬爾科夫鏈蒙特卡羅是一種基于當(dāng)前值重復(fù)繪制某一分布參數(shù)隨機(jī)值的方法。每個(gè)值的樣本都是隨機(jī)的，但是值的選擇受限于當(dāng)前狀態(tài)和假定的參數(shù)先驗(yàn)分布。MCMC 可以被認(rèn)為是一種隨機(jī)游走，在這個(gè)過程中逐漸收斂到真實(shí)分布。

為了繪制 α 和 β 的隨機(jī)值，我們需要假設(shè)這些值的先驗(yàn)分布。由于我們對參數(shù)沒有任何提前的假設(shè)，我們可以使用正態(tài)分布。正態(tài)分布也稱高斯分布，它由均值和方差定義，分別顯示數(shù)據(jù)的位置以及擴(kuò)散情況。下圖是具有不同均值和方差的幾種正態(tài)分布：

我們所使用的 MCMC 算法被稱為 Metropolis Hastings。為了將我們觀察的數(shù)據(jù)與模型聯(lián)系起來，每繪制一組隨機(jī)值，算法會根據(jù)數(shù)據(jù)對其進(jìn)行評估。如果隨機(jī)值與數(shù)據(jù)不一致（這里稍微進(jìn)行了一些簡化），這些值將被拒絕，模型保持當(dāng)前狀態(tài)。反之，如果隨機(jī)值與數(shù)據(jù)一致，這些值將會分配給參數(shù)并成為當(dāng)前狀態(tài)。該過程將持續(xù)進(jìn)行指定的步驟數(shù)目，模型的準(zhǔn)確率也隨著步驟數(shù)量的增加而改善。

綜合而言，馬爾科夫鏈蒙特卡羅在我們的問題當(dāng)中基本步驟如下：

為 logistic 函數(shù)選擇一組初始參數(shù) α 和 β。

根據(jù)當(dāng)前狀態(tài)，把新的隨機(jī)值分配給 α 和 β。

檢查新的隨機(jī)值是否與觀察結(jié)果一致。如果不一致，拒絕這些隨機(jī)值并返回前一個(gè)狀態(tài)。如果一致，則接受這些值，將其作為新的當(dāng)前狀態(tài)。

對指定的迭代次數(shù)重復(fù)執(zhí)行步驟 2 和 3。

該算法將返回它為 α 和 β 生成的所有值。然后，我們可以使用這些值的平均值作為 logistc 函數(shù)中 α 和 β 的最終可能值。MCMC 無法返回「真實(shí)」值，它給出的是分布的近似值。給定數(shù)據(jù)的情況下，最終輸出的睡眠概率模型將是具有 α 和 β 均值的 logistic 函數(shù)。

Python 實(shí)現(xiàn)

上述細(xì)節(jié)在我腦海中徘徊已久，最后終于在 Python 中進(jìn)行了實(shí)現(xiàn)！親眼看到第一手的結(jié)果比讀取別人的描述有幫助得多。要在 Python 中實(shí)現(xiàn) MCMC，我們需要使用 PyMC3 貝葉斯推理庫。它將大部分細(xì)節(jié)進(jìn)行了抽象，從而讓我們能不迷失在理論中，并建立我們的模型。

下面的代碼創(chuàng)建的模型帶有參數(shù) α 和 β、概率 p 和觀察結(jié)果 observed。step 變量指的是特定的算法，sleep_trace 則保存了模型生成的所有參數(shù)值。

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);

（請?jiān)?notebook 中查閱完整代碼）

為了了解運(yùn)行此代碼時(shí)發(fā)生的情況，我們可以查看模型運(yùn)行過程中生成的 α 和 β 的所有值。

它們被稱為軌跡圖。我們可以看到，每個(gè)狀態(tài)都與之前的狀態(tài)有關(guān)（馬爾科夫鏈），但是這些值波動(dòng)顯著（蒙特卡羅采樣）。

在 MCMC 中，通常高達(dá) 90% 的軌跡會被拋棄。該算法無法立即收斂到真正的分布，且初始值往往并不準(zhǔn)確。后期的參數(shù)值通常更好，這意味著它們是適用于建立模型的參數(shù)。我們使用了 10000 個(gè)樣本并丟棄了前 50%，但是一個(gè)行業(yè)的應(yīng)用可能會使用數(shù)十萬甚至上百萬個(gè)樣本。

給定足夠多的迭代次數(shù)，MCMC 將收斂于真實(shí)值。但是，對收斂進(jìn)行評估可能比較困難。對此我將不在本文討論（一個(gè)方法是測量軌跡的自相關(guān)），但是，如果我們想要結(jié)果最準(zhǔn)確，這是一個(gè)重要的考慮因素。PyMC3 建立了評估模型好壞的函數(shù)，其中包括軌跡圖和自相關(guān)圖。

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

軌跡圖（左）和自相關(guān)圖（右）

睡眠模型

最終建立并運(yùn)行模型之后，是時(shí)候使用結(jié)果了。我們將最后 5000 個(gè) α 和 β 樣本的平均值作為參數(shù)最可能的值，這就讓我們能夠創(chuàng)建一條曲線，建模睡眠后驗(yàn)概率：

該模型能很好地反映數(shù)據(jù)的結(jié)果。此外，它捕捉了我睡眠模式當(dāng)中的固有變化。該模型給出的不是一個(gè)簡單的是非答案，而是一個(gè)概率。例如，我們可以通過該模型找到在給定時(shí)間我睡著的概率，并能找到睡眠概率經(jīng)過 50% 的時(shí)間：

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 點(diǎn)上床睡覺，但這顯然不是大多數(shù)下的實(shí)際情況。我們可以發(fā)現(xiàn)，我上床的平均時(shí)間是晚上 10:14 左右。

在數(shù)據(jù)給定的情況下，這些值是最有可能的估計(jì)值。然而，因?yàn)槟Ｐ捅旧硎墙频?，所以存在與這些概率相關(guān)的不確定性。為了表示這種不確定性，我們可以使用所有的 α 和 β 樣本（而不是它們的平均值）來預(yù)測某一給定時(shí)間的睡眠概率，然后據(jù)此繪制直方圖。

這些結(jié)果更好地反映了 MCMC 模型真正做了什么。MCMC 找到的不是一個(gè)簡單的答案，而是可能值的樣本。貝葉斯推理在現(xiàn)實(shí)世界中起到了重要作用，是因?yàn)樗鼜母怕实慕嵌缺硎绢A(yù)測結(jié)果。我們可以說，問題會有一個(gè)可能性最大的答案，但是更加準(zhǔn)確的回應(yīng)是任何預(yù)測都存在一系列的可能值。

喚醒模型

我可以使用描述早晨醒來時(shí)間的數(shù)據(jù)建立一個(gè)類似的模型。我定了一個(gè)鬧鐘，努力在早晨 6:00 起床。但是可以看到，我并非每日都是如此。下圖展現(xiàn)了我從入睡到醒來過渡過程的最終模型以及觀察數(shù)據(jù)。

通過查詢模型，我們可以找出在給定時(shí)間我睡著的概率以及最有可能醒來的時(shí)間。

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.

看起來我得處理一下我的鬧鐘了！

睡眠時(shí)長

出于好奇心和練習(xí)目的，我最終想創(chuàng)造的是關(guān)于我睡眠時(shí)長的模型。首先，我們需要找到一個(gè)函數(shù)來模擬數(shù)據(jù)的分布。我猜想結(jié)果應(yīng)該會是正態(tài)分布的形式，但是我們只有通過檢查數(shù)據(jù)才能得到最終結(jié)果。

正態(tài)分布確實(shí)可行，但這無法捕捉右側(cè)的偏離點(diǎn)（即我睡眠時(shí)間非常長的情況）。我們可以用兩個(gè)獨(dú)立的正態(tài)分布來表示兩個(gè)模型，但是，我想使用偏正態(tài)分布。偏正態(tài)分布有三個(gè)參數(shù)：均值、方差、偏斜度 α。以上三個(gè)參數(shù)都需要通過 MCMC 來學(xué)習(xí)。下面的代碼建立了上述模型，并進(jìn)行了 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)

現(xiàn)在，我們可以使用這三個(gè)參數(shù)的均值來構(gòu)建最有可能的分布。以下是根據(jù)數(shù)據(jù)觀察得到的最終偏正態(tài)分布。

看起來擬合得不錯(cuò)！我們可以通過查詢模型，找到我至少可以獲得一定睡眠時(shí)長的概率，同時(shí)也能找到最可能的睡眠時(shí)長：

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.

我對這個(gè)結(jié)果并不是完全滿意，但是對一個(gè)研究生而言，這樣的結(jié)果已經(jīng)不錯(cuò)啦。

小結(jié)

這個(gè)項(xiàng)目的順利完成再次展示了解決問題的重要性，而且我們最好選擇解決真實(shí)存在的應(yīng)用問題（https://towardsdatascience.com/how-to-master-new-skills-656d42d0e09c）。在使用馬爾科夫鏈蒙特卡羅構(gòu)建貝葉斯推理的端對端實(shí)現(xiàn)過程中，我學(xué)習(xí)了許多基礎(chǔ)知識，而且非常享受這個(gè)過程。我不僅更加了解我的習(xí)慣（以及我需要改進(jìn)的方面），而且終于弄明白了 MCMC 和貝葉斯推理到底是什么。在數(shù)據(jù)科學(xué)領(lǐng)域，我們要不斷地給自己的庫存知識增加新的工具，而最有效的學(xué)習(xí)方法就是找到一個(gè)問題并著手開始解決它！