摘要： 本文通過用Python中的馬爾可夫鏈蒙特卡羅實現(xiàn)了睡眠模型項目，并教會如何使用MCMC。

在過去的幾個月里，我在數(shù)據(jù)科學(xué)領(lǐng)域里遇到一個術(shù)語：馬爾可夫鏈蒙特卡羅（MCMC）。在博客或文章里，每次看到這個語，我都會搖搖頭，有幾次我試著學(xué)習(xí)MCMC和貝葉斯推理，但每次一開始，就很快放棄了。我學(xué)習(xí)新技術(shù)的方式都是把它應(yīng)用到一個實際問題上。

通過使用一些數(shù)據(jù)和一本應(yīng)用實戰(zhàn)的書（Bayesian Methods for Hackers），我終于通過一個實際項目弄懂了MCMC。像往常一樣，當(dāng)把這些技術(shù)概念應(yīng)用到實際問題中時，理解它們要比閱讀書上的抽象概念更容易。本文通過介紹Python中的MCMC實現(xiàn)過程，最終教會了我使用這個強大的建模和分析工具。

本項目的完整代碼和相關(guān)數(shù)據(jù)在GitHub上可以找到。本文重點討論了應(yīng)用程序和結(jié)果，涵蓋了很多有深度的內(nèi)容。

介紹

實際生活中的數(shù)據(jù)永遠(yuǎn)不是完美的，但我們?nèi)匀豢梢酝ㄟ^正確的模型從噪音數(shù)據(jù)中提取有價值的信息。

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

本項目的目的是使用睡眠數(shù)據(jù)創(chuàng)建一個模型，該模型指定了睡眠的后驗概率作為一個時間函數(shù)。由于時間是一個連續(xù)變量，因為指定整個后驗分布是比較困難的，所以我們轉(zhuǎn)向接近于分布的方法-MCMC。

概率分布的選擇

在使用MCMC之前，我們需要確定適當(dāng)?shù)暮瘮?shù)來給睡眠的后驗概率分布建模。最簡單方法是通過可視化的方式進(jìn)行數(shù)據(jù)檢查，入睡的時間函數(shù)的觀察結(jié)果如下圖所示：

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

在圖上，每個數(shù)據(jù)節(jié)點被標(biāo)記為一個點，點的強度表明了在指定時間的觀測次數(shù)。我的表只記錄入睡的時間，所以為了擴大數(shù)據(jù)量，我在時間點兩邊的每一分鐘都加上了節(jié)點。如果手表記錄我在晚上10：05分睡著了，那么之前的每分鐘都被表示為0（清醒），之后的每分鐘都表示為1（入睡）。這將會把大約60個晚上的觀測擴展到11340個數(shù)據(jù)節(jié)點。

可以看到，我總是于晚上10:00之后入睡，我們想要創(chuàng)建一個獲取從醒到睡的概率轉(zhuǎn)換的模型。當(dāng)模型在一個精確的時間從清醒（0）到入睡（1）轉(zhuǎn)換的時候，我們可以給模型用一個簡單的階梯函數(shù)，因為我不會每晚都在同一時間入睡，所以我們需要一個函數(shù)來把轉(zhuǎn)換過程建模為漸變過程。給定數(shù)據(jù)的最佳選擇是一個在0和1之間平穩(wěn)轉(zhuǎn)換的邏輯函數(shù)。以下是入睡概率作為時間函數(shù)的邏輯方程。

β（beta）和α（alpha）是我們在MCMC過程中必須學(xué)習(xí)的模型參數(shù)。如下所示具有不同參數(shù)的邏輯函數(shù)：

邏輯函數(shù)適合這些數(shù)據(jù)，因為入睡轉(zhuǎn)換的概率是逐漸的，同時獲取了我的睡眠樣本的變化。我們希望能為函數(shù)增加一個時間變量t，并找出入睡的概率，它必須在0和1之間。為了創(chuàng)建這個模型，我們通過一個分類的技術(shù)作為MCMC，用數(shù)據(jù)來找到最佳的α和β參數(shù)。

馬爾可夫鏈蒙特卡羅（MCMC）

MCMC是從概率分布中抽樣以構(gòu)造最可能分布的一類方法。我們不能直接計算邏輯分布，所以我們?yōu)楹瘮?shù)的參數(shù)（α和β）生成數(shù)千個稱為樣本的數(shù)值，以創(chuàng)建分布的近似值。MCMC背后的思想是，隨著生成更多的樣本，近似值會越來越接近實際的分布。

MCMC方法有兩個部分，蒙特卡羅（Monte Carlo）是指使用重復(fù)隨機樣本來獲得數(shù)值結(jié)果的通用技術(shù)。蒙特卡羅可以被認(rèn)為是進(jìn)行很多次的實驗，每次改變模型中的變量并觀察反應(yīng)。通過選擇隨機值，我們可以研究參數(shù)空間的一大部分，可能的變量值的范圍。參數(shù)空間，如下圖所示：

很顯然，我們不能嘗試圖中的每個點，但是通過從較高概率（紅色）的區(qū)域隨機抽樣，我們可以創(chuàng)建最可能的模型。

馬爾可夫鏈（Markov Chain）

馬爾可夫鏈是當(dāng)前狀態(tài)被下一個狀態(tài)所依賴的過程。馬爾可夫鏈?zhǔn)遣挥涗洜顟B(tài)的，因為只有當(dāng)前的狀態(tài)才是重要的，而它如何到達(dá)那個狀態(tài)并不重要。如果這有點難以理解，那么考慮一個日?，F(xiàn)象，天氣。如果我們想預(yù)測明天的天氣，只能根據(jù)今天的天氣來進(jìn)行一個合理的預(yù)測。如果今天下雪了，我們查看一下下雪后第二天的天氣分布的歷史數(shù)據(jù)，以預(yù)測明天是什么天氣的概率。馬爾可夫鏈的概念是，我們不需要一定知道一個過程的整個歷史數(shù)據(jù)來預(yù)測下一個輸出，相似的現(xiàn)象在許多現(xiàn)實情況中都能很好地預(yù)測。

MCMC結(jié)合了馬爾可夫鏈和蒙特卡羅的思想，是一種基于當(dāng)前值對分布的參數(shù)重復(fù)抽取隨機值的方法。每個值的樣本都是隨機的，但是這些值的選擇受到當(dāng)前狀態(tài)和參數(shù)的假定先驗分布的限制。MCMC可以看作是逐漸收斂到實際分布的隨機游動。

為了抽取α和β的隨機值，我們需要假設(shè)這些值的先驗分布。由于我們沒有預(yù)先給參數(shù)假設(shè)，所以可以使用一個正態(tài)分布。正態(tài)分布或高斯分布由平均值定義，表明了數(shù)據(jù)的位置、方差和分布范圍。下圖是具有不同平均值和范圍的幾個正態(tài)分布:

我們使用的特定MCMC算法稱為Metropolis Hastings。為了將觀測數(shù)據(jù)與模型聯(lián)系起來，每次繪制一組隨機值時，該算法都會根據(jù)數(shù)據(jù)進(jìn)行評估。如果這些隨機值與數(shù)據(jù)不一致，則會被拒絕，并且模型保持在當(dāng)前狀態(tài)。反之，如果與數(shù)據(jù)一致，則將這些隨機值分配給參數(shù)并變成當(dāng)前狀態(tài)。這個過程持續(xù)一定數(shù)量的迭代，那么模型的精確度就會隨著迭代數(shù)量的增加而提高。

綜上所述，用MCMC解決問題的基本步驟如下:

（1）選擇α和β以及邏輯函數(shù)的參數(shù)初始值集合；

（2）根據(jù)當(dāng)前狀態(tài)給α和β隨機分配新的值；

（3）檢查新的隨機值是否符合觀測值。如果不符合，則拒絕這些隨機值，并返回到先前狀態(tài)，反之，要是符合，則接受，并更新為新的當(dāng)前狀態(tài)；

（4）如需迭代，則重復(fù)步驟2和3；

算法返回其生成的所有α和β的值。然后，我們可以使用這些值的平均值作為邏輯函數(shù)中最有可能的α和β最終值。MCMC不能返回“True”值，而是返回一個分布的近似值。由已有數(shù)據(jù)得到的睡眠概率，其最終模型是具有α和β平均值的邏輯函數(shù)。

Python的實現(xiàn)

為了在Python中實現(xiàn)MCMC，我們將使用貝葉斯推理庫PyMC3，它對大部分細(xì)節(jié)進(jìn)行了抽象，使我們能夠創(chuàng)建模型而不是空有理論。

下面的代碼創(chuàng)建了具有參數(shù)α和β、概率p和觀察值的完整模型。步驟變量引用特定的算法，并且變量sleep_trace保存了由模型生成的所有參數(shù)值。

為了更好地了解代碼的運行，我們可以看下所有的模型運行過程中所產(chǎn)生的α和β值。

上圖被稱為軌跡圖。我們可以看到，每個狀態(tài)都與先前的馬爾可夫鏈相關(guān)，但其值在蒙特卡羅抽樣中振蕩明顯。

在MCMC中，很常見的做法是丟棄多達(dá)90%的軌跡。算法不立即收斂到真實分布，而且初始值往往也不準(zhǔn)確。之后的參數(shù)值通常會更好一些，這意味著應(yīng)該用它們來構(gòu)建模型。我們使用了10000個樣本，并丟棄了前面的50%，但是在企業(yè)級應(yīng)用中可能會使用成百上千個或數(shù)百萬個樣本。

MCMC收斂到真實值，但評估收斂可能是很困難的。如果我們想要最精確的結(jié)果，這是一個重要的因素。PyMC3庫已經(jīng)創(chuàng)建了用于評估模型質(zhì)量的函數(shù)，包括軌跡圖和自相關(guān)圖。

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

睡眠模型

在最終建立和運行模型之后，現(xiàn)在應(yīng)該使用了。我們將最后5000個α和β樣本的平均值作為最有可能的參數(shù)值，這允許我們創(chuàng)建單個曲線圖來模擬指定時間之后的入睡概率:

該模型能很好地表示數(shù)據(jù)，還獲取了我入睡模式固有的變化，它不是給出一個結(jié)果，而是給出一個概率。例如，可以查詢該模型以找出我在指定時間入睡的概率，并找到概率超過50%的時間點:

晚上9點半入睡概率: 4.80%.

晚上9點半入睡概率: 27.44%.

晚上10點入睡概率: 73.91%.

在晚上10點14分，入睡概率提高到了50%以上。

可以看到我入睡的平均時間是晚上10點14分左右。

這些值是給定數(shù)據(jù)的最可能的估計值。然而，會有這些概率相關(guān)的不確定性，因為模型是近似的。為了表示這種不確定性，我們可以在一個給定的時間點使用所有α和β的樣本來進(jìn)行入睡概率預(yù)測，然后根據(jù)結(jié)果繪制柱狀：

上述結(jié)果給出了一個更好的MCMC模型能實際做到的指標(biāo)。這個方法并沒有找到一個正確答案，而是可能值的一個樣本。貝葉斯推理是有實際用處的，因為它以概率的形式表示預(yù)測的結(jié)果。我們可以說得到一個最有可能的答案，但是更準(zhǔn)確的說法應(yīng)該是任何預(yù)測都是一系列值的范圍。

睡醒模型

可以用我早上睡醒的相關(guān)數(shù)據(jù)來找到一個類似的模型。我通常在早上6點起床，下圖根據(jù)觀察結(jié)果顯示了從睡覺到睡醒的最終模型。

可以通過模型來發(fā)現(xiàn)我在某一指定時間入睡的概率和最有可能睡醒的時間。

早上5點半睡醒的概率: 14.10%.

早上6點睡醒的概率: 37.94%.

早上6點半睡醒的概率: 69.49%.

在早上6點11分睡醒的概率超過50%。

入睡持續(xù)時間

最后一個我想創(chuàng)建的是入睡時間模型。首先，我們需要找到一個函數(shù)來模擬數(shù)據(jù)的分布，但只能通過檢查數(shù)據(jù)來找到。

一個普通的分布即可實現(xiàn)，但它不會獲取右邊的離群點?？梢允褂脙蓚€相互獨立的正態(tài)分布來表示這兩種模式，但是，我會使用一個偏態(tài)分布。偏態(tài)分布有三個參數(shù)，平均值、方差、偏差，并且這三個參數(shù)必須從MCMC算法中進(jìn)行學(xué)習(xí)。下面的代碼創(chuàng)建模型并實現(xiàn)Metropolis Hastings抽樣：

現(xiàn)在，我們可以使用三個參數(shù)的平均值來構(gòu)造最有可能的分布。下圖表示數(shù)據(jù)的最終偏態(tài)分布:

上圖看起來很合乎實際情況，通過查詢模型可以發(fā)現(xiàn)我至少獲得一定的入睡持續(xù)時間，和最可能的入睡時間范圍的概率：

入睡至少持續(xù)6.5小時的概率= 99.16%.

入睡至少持續(xù)8小時的概率= 44.53%.

入睡至少持續(xù)9小時的概率= 10.94%.

入睡最可能持續(xù)的時間是7.67小時

結(jié)論

通過完成這個項目，告訴我們解決問題的重要性是最好解決實際生活中的問題。數(shù)據(jù)科學(xué)需要不斷地進(jìn)行知識積累，最有效的方法就是有了問題盡快開始。

文章原標(biāo)題《Markov Chain Monte Carlo in Python》