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

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

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

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

簡介

我的 Garmin Vivosmart 手表可以根據心率和運動情況追蹤我的睡眠和起床狀況。它并非 100％準確，不過真實數據從不完美，我們仍然可以借助正確的模型從噪聲數據中提取有用的知識！

典型的睡眠數據

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

選擇一個概率分布

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

睡眠數據

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

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

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

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

馬爾科夫鏈蒙特卡羅

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

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

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

馬爾科夫鏈

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

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

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

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

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

為 logistic 函數選擇一組初始參數 α 和 β。

根據當前狀態，把新的隨機值分配給 α 和 β。

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

對指定的迭代次數重復執行步驟 2 和 3。

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

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

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

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

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

在 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 左右。

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

這些結果更好地反映了 MCMC 模型真正做了什么。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.

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

睡眠時長

出于好奇心和練習目的，我最終想創造的是關于我睡眠時長的模型。首先，我們需要找到一個函數來模擬數據的分布。我猜想結果應該會是正態分布的形式，但是我們只有通過檢查數據才能得到最終結果。

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

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

小結

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