大資料文摘作品

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

在過去幾個月裡，我在資料科學的世界裡反覆遇到一個詞：馬爾可夫鏈蒙特卡洛（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）和蒙特卡洛（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中實施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