用Python入門不明覺厲的馬爾可夫鏈蒙特卡羅（附案例代碼）

大數據文摘作品

編譯：Niki、張南星、Shan LIU、Aileen

這篇文章讓小白也能讀懂什麼是人們常說的Markov Chain Monte Carlo。

在過去幾個月里，我在數據科學的世界裡反覆遇到一個詞：馬爾可夫鏈蒙特卡洛（Markov Chain Monte Carlo , MCMC）。在我的研究室、podcast和文章里，每每遇到這個詞我都會「不明覺厲」地點點頭，覺得這個演算法聽起來很酷，但每次聽人提起也只是有個模模糊糊的概念。

我屢次嘗試學習MCMC和貝葉斯推論，而一拿起書，又很快就放棄了。無奈之下，我選擇了學習任何新東西最佳的方法：應用到一個實際問題中。

通過使用一些我曾試圖分析的睡眠數據和一本實操類的、基於應用教學的書（《寫給開發者的貝葉斯方法》，我最終通過一個實際項目搞明白了MCMC。

《寫給開發者的貝葉斯方法》

https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers

和學習其他東西一樣，當我把這些技術性的概念應用於一個實際問題中而不是單純地通過看書去了解這些抽象概念，我更容易理解這些知識，並且更享受學習的過程。

這篇文章介紹了馬爾可夫鏈蒙特卡洛在Python中入門級的應用操作，這個實際應用最終也使我學會使用這個強大的建模分析工具。

此項目全部的代碼和數據：

https://github.com/WillKoehrsen/ai-projects/blob/master/bayesian/bayesian_inference.ipynb

這篇文章側重於應用和結果，因此很多知識點只會粗淺的介紹，但對於那些想了解更多知識的讀者，在文章也嘗試提供了一些學習鏈接。

案例簡介

我的智能手環在我入睡和起床時會根據心率和運動進行記錄。它不是100%準確的，但現實世界中的數據永遠不可能是完美的，不過我們依然可以運用正確的模型從這些雜訊數據提取出有價值的信息。

典型睡眠數據

這個項目的目的在於運用睡眠數據建立一個能夠確立睡眠相對於時間的後驗分布模型。由於時間是個連續變數，我們無法知道後驗分布的具體表達式，因此我們轉向能夠近似後驗分布的演算法，比如馬爾可夫鏈蒙特卡洛（MCMC）。

選擇一個概率分布

在我們開始MCMC之前，我們需要為睡眠的後驗分布模型選擇一個合適的函數。一種簡單的做法是觀察數據所呈現的圖像。下圖呈現了當我入睡時時間函數的數據分布。

睡眠數據

每個數據點用一個點來表示，點的密度展現了在固定時刻的觀測個數。我的智能手錶只記錄我入睡的那個時刻，因此為了擴展數據，我在每分鐘的兩端添加了數據點。如果我的手錶記錄我在晚上10：05入睡，那麼所有在此之前的時間點被記為0（醒著），所有在此之後的時間點記為1（睡著）。這樣一來，原本大約60夜的觀察量被擴展成11340個數據點。

可以看到我趨向於在10：00後幾分鐘入睡，但我們希望建立一個把從醒到入睡的轉變用概率進行表達的模型。我們可以用一個簡單的階梯函數作為模型，在一個精確時間點從醒著（0）變到入睡（1），但這不能代表數據中的不確定性。

我不會每天在同一時間入睡，因此我們需要一個能夠模擬出這個個漸變過程的函數來展現變化當中的差異性。在現有數據下最好的選擇是logistic函數，在0到1之前平滑地移動。下面這個公式是睡眠狀態相對時間的概率分布，也是一個logistic公式。

在這裡，β (beta) 和 α (alpha) 是模型的參數，我們只能通過MCMC去模擬它們的值。下面展示了一個參數變化的logistic函數。

一個logistic函數能夠很好的擬合數據，因為在logistic函數中入睡的概率在逐漸改變，捕捉了我睡眠模式的變化性。我們希望能夠帶入一個具體的時間t到函數中，從而得到一個在0到1之間的睡眠狀態的概率分布。我們並不會直接得到我是否在10：00睡著了的準確答案，而是一個概率。創建這個模型，我們通過數據和馬爾可夫鏈蒙特卡洛去尋找最優的alpha和beta係數估計。

馬爾可夫鏈蒙特卡洛

馬爾可夫鏈蒙特卡羅是一組從概率分布中抽樣，從而建立最近似原分布的函數的方法。因為我們不能直接去計算logistic分布，所以我們為係數(alpha 和 beta)生成成千上萬的數值-被稱為樣本-去建立分布的一個模擬。MCMC背後的基本思想就是當我們生成越多的樣本，我們的模擬就更近似於真實的分布。

馬爾可夫鏈蒙特卡洛由兩部分組成。蒙特卡洛代表運用重複隨機的樣本來獲取一個準確答案的一種模擬方法。蒙特卡洛可以被看做大量重複性的實驗，每次更改變數的值並觀察結果。通過選擇隨機的數值，我們可以在係數的範圍空間，也就是變數可取值的範圍，更大比例地探索。下圖展示了在我們的問題中，一個使用高斯分布作為先驗的係數空間。

能夠清楚地看到我們不能在這些圖中一一找出單個的點，但通過在更高概率的區域（紅色）進行隨機抽樣，我們就能夠建立最近似的模型。

馬爾可夫鏈（Markov Chain）

馬爾可夫鏈是一個「下個狀態值只取決於當前狀態」的過程。（在這裡，一個狀態指代當前時間係數的數值分配）。一個馬爾可夫鏈是「健忘」的，因為如何到達當前狀態並不要緊，只有當前的狀態值是關鍵。如果這有些難以理解的話，讓我們來設想一個每天都會經歷的情景--天氣。

如果我們希望預測明天的天氣，那麼僅僅使用今天的天氣狀況我們就能夠得到一個較為合理的預測。如果今天下雪，我們可以觀測有關下雪後第二天天氣的歷史數據去預測明天各種天氣狀況的可能性。馬爾可夫鏈的定義就是我們不需要知道一個過程中的全部歷史狀態去預測下一節點的狀態，這種近似在許多現實問題中都很有用。

把馬爾可夫鏈（Markov Chain）和蒙特卡洛（Monte Carlo），兩者放到一起，就有了MCMC。MCMC是一種基於當前值，重複為概率分布係數抽取隨機數值的方法。每個樣本都是隨機的，但是數值的選擇也受當前值和係數先驗分布的影響。MCMC可以被看做是一個最終趨於真實分布的隨機遊走。

為了能夠抽取alpha 和 beta的隨機值，我們需要為每個係數假設一個先驗分布。因為我們沒有對於這兩個係數的任何假設，我們可以使用正太分布作為先驗。正太分布，也稱高斯分布，是由均值（展示數據分布），和方差（展示離散性）來定義的。下圖展示了多個不同均值和離散型的正態分布。

具體的MCMC演算法被稱作Metropolis Hastings。為了連接我們的觀察數據到模型中，每次一組隨機值被抽取，演算法將把它們與數據進行比較。一旦它們與數據不吻合（在這裡我簡化了一部分內容），這些值就會被捨棄，模型將停留在當前的狀態值。

如果這些隨機值與數據吻合，那麼這些值就被接納為各個係數新的值，成為當前的狀態值。這個過程會有一個提前設置好的迭代次數，次數越多，模型的精確度也就越高。

把上面介紹的整合到一起，就能得到在我們的問題中所需進行的最基本的MCMC步驟：

• 為logistic函數的係數alpha 和beta選擇初始值。

• 基於當前狀態值隨機分配給alpha 和beta新的值。

• 檢查新的隨機值是否與觀察數據吻合。如果不是，捨棄掉這個值，並回到上一狀態值。如果吻合，接受這個新的值作為當前狀態值。

• 重複步驟2和3（重複次數提前設置好）。

• 這個演算法會給出所有它所生成的alpha 和beta值。我們可以用這些值的平均數作為alpha 和beta在logistic函數中可能性最大的終值。MCMC不會返回「真實」的數值，而是函數分布的近似值。睡眠狀態概率分布的最終模型將會是以alph和beta均值作為係數的logistic函數。

Python實施

我再三思考模擬上面提到的細節，最終我開始用Python將它們變成現實。觀察一手的結果會比閱讀別人的經驗貼有幫助得多。想要在Python中實施MCMC，我們需要用到PyMC3貝葉斯庫，它省略了很多細節，方便我們創建模型，避免迷失在理論之中。

通過下面的這些代碼可以創建完整的模型，其中包含了參數alpha 、beta、概率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);

為了更直觀地展現代碼運行的效果，我們可以看一下模型運行時alpha和beta生成的值。

這些圖叫做軌跡圖，可以看到每個狀態都與其歷史狀態相關，即馬爾可夫鏈；同時每個值劇烈波動，即蒙特卡洛抽樣。

使用MCMC時，常常需要放棄軌跡圖中90%的值。這個演算法並不能立即展現真實的分布情況，最初生成的值往往是不準確的。隨著演算法的運行，後產生的參數值才是我們真正需要用來建模的值。我使用了一萬個樣本，放棄了前50%的值，但真正在行業中應用時，樣本量可達成千上萬個、甚至上百萬個。

通過足夠多的迭代，MCMC逐漸趨近於真實的值，但是估算收斂性並不容易。這篇文章中並不會涉及到具體的估算方法（方法之一就是計算軌跡的自我相關性），但是這是得到最準確結果的必要條件。PyMC3的函數能夠評估模型的質量，包括對軌跡、自相關圖的評估。

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

軌跡圖（左）和自相關性圖（右）

睡眠模型

建模、模型運行完成後，該最終結果上場了。我們將使用最終的5000個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。

這些值是基於樣本數據的最有可能值，但這些概率值都有一定的不確定性，因為模型本身就是近似的。為了展現這種不確定性，我們可以使用所有的alpha、beta值來估計某個時間點的睡覺概率，而不是使用平均值，並且把這些概率值展現在圖中。

晚上10:00睡覺的概率分布

這些結果能夠更好地展現MCMC模型真正在做的事情，即它並不是在尋找單一的答案，而是一系列可能值。貝葉斯推論在現實世界中非常有用，因為它是對概率進行了預測。我們可以說存在一個最可能的答案，但其實更準確的回復應當是：每個預測都有一系列的可能值。

起床模型

同樣我可以用我的起床數據創建類似的模型。我希望能夠在鬧鐘的幫助下總能在早上6:00起床，但實際上並不如此。下面這張圖展現了基於觀測值我起床的最終模型：

基於5000個樣本的起床事後概率

可以通過模型得出我在某個特定時間起床的概率，以及我最有可能起床的時間：

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.

看來我需要一個更生猛的鬧鐘了….

睡眠的時間

出於好奇以及實踐需求，最後我想創建的模型是我的睡眠時間模型。首先，我們需要尋找到一個描述數據分布的函數。事先，我認為應該是正態函數，但無論如何我們需要用數據來證明。

睡眠時間長短分布

正態分布的確能夠解釋大部分數據，但是圖中右側的異常值卻無法得到解釋（當我睡懶覺的時候）。我們可以用兩個單獨的正態分布來代表兩種模式，但我要用偏態分布。偏態分布有三個參數：平均值、偏離值，以及alpha傾斜值。這三個參數的值都需要從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.

結論

我想再次強調，完成這個項目讓我體會到解決問題的重要性，尤其是有現實應用意義的項目！在我嘗試使用馬爾可夫鏈蒙特卡洛來端到端建立貝葉斯推論的時候，我重新熟悉了許多基礎知識，並且非常享受這個過程。

我不僅了解到自身需要改進的習慣，而且當別人在談論MCMC和貝葉斯推論時，我終於真的明白他們在談論什麼了！數據科學正是關於持續不斷地在你自己的知識庫中輸入新的工具，而這最有效的辦法就是發現一個問題，然後應用你所學的去解決問題！

原文鏈接：

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

推薦閱讀：