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

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

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

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

簡介

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

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

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

選擇一個概率分布

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

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

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

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

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

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

馬爾科夫鏈蒙特卡羅

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

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

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

馬爾科夫鏈

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

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

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

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

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

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

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

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

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

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

Python 實現(xiàn)

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

下面的代碼創(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);

（請在 notebook 中查閱完整代碼）

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

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

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

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

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

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

睡眠模型

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

該模型能很好地反映數(shù)據(jù)的結(jié)果。此外，它捕捉了我睡眠模式當中的固有變化。該模型給出的不是一個簡單的是非答案，而是一個概率。例如，我們可以通過該模型找到在給定時間我睡著的概率，并能找到睡眠概率經(jīng)過 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 點上床睡覺，但這顯然不是大多數(shù)下的實際情況。我們可以發(fā)現(xiàn)，我上床的平均時間是晚上 10:14 左右。

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

這些結(jié)果更好地反映了 MCMC 模型真正做了什么。MCMC 找到的不是一個簡單的答案，而是可能值的樣本。貝葉斯推理在現(xiàn)實世界中起到了重要作用，是因為它從概率的角度表示預(yù)測結(jié)果。我們可以說，問題會有一個可能性最大的答案，但是更加準確的回應(yīng)是任何預(yù)測都存在一系列的可能值。

喚醒模型

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

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

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.

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

睡眠時長

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

正態(tài)分布確實可行，但這無法捕捉右側(cè)的偏離點（即我睡眠時間非常長的情況）。我們可以用兩個獨立的正態(tài)分布來表示兩個模型，但是，我想使用偏正態(tài)分布。偏正態(tài)分布有三個參數(shù)：均值、方差、偏斜度 α。以上三個參數(shù)都需要通過 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)

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

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

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.

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

小結(jié)

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