用Python實現馬爾可夫鏈蒙特卡羅

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

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

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

本項目的完整代碼和相關數據在GitHub上能夠找到。本文重點討論了應用程序和結果，涵蓋了不少有深度的內容。

介紹

實際生活中的數據永遠不是完美的，但咱們仍然能夠經過正確的模型從噪音數據中提取有價值的信息。

典型睡眠數據

本項目的目的是使用睡眠數據建立一個模型，該模型指定了睡眠的後驗機率做爲一個時間函數。因爲時間是一個連續變量，由於指定整個後驗分佈是比較困難的，因此咱們轉向接近於分佈的方法-MCMC。

機率分佈的選擇

在使用MCMC以前，咱們須要肯定適當的函數來給睡眠的後驗機率分佈建模。最簡單方法是經過可視化的方式進行數據檢查，入睡的時間函數的觀察結果以下圖所示：

 

睡眠數據

在圖上，每一個數據節點被標記爲一個點，點的強度代表了在指定時間的觀測次數。個人表只記錄入睡的時間，因此爲了擴大數據量，我在時間點兩邊的每一分鐘都加上了節點。若是手錶記錄我在晚上10：05分睡着了，那麼以前的每分鐘都被表示爲0（清醒），以後的每分鐘都表示爲1（入睡）。這將會把大約60個晚上的觀測擴展到11340個數據節點。

能夠看到，我老是於晚上10:00以後入睡，咱們想要建立一個獲取從醒到睡的機率轉換的模型。當模型在一個精確的時間從清醒（0）到入睡（1）轉換的時候，咱們能夠給模型用一個簡單的階梯函數，由於我不會每晚都在同一時間入睡，因此咱們須要一個函數來把轉換過程建模爲漸變過程。給定數據的最佳選擇是一個在0和1之間平穩轉換的邏輯函數。如下是入睡機率做爲時間函數的邏輯方程。

 

β（beta）和α（alpha）是咱們在MCMC過程當中必須學習的模型參數。以下所示具備不一樣參數的邏輯函數：

 

邏輯函數適合這些數據，由於入睡轉換的機率是逐漸的，同時獲取了個人睡眠樣本的變化。咱們但願能爲函數增長一個時間變量t，並找出入睡的機率，它必須在0和1之間。爲了建立這個模型，咱們經過一個分類的技術做爲MCMC，用數據來找到最佳的α和β參數。

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

 

MCMC是從機率分佈中抽樣以構造最可能分佈的一類方法。咱們不能直接計算邏輯分佈，因此咱們爲函數的參數（α和β）生成數千個稱爲樣本的數值，以建立分佈的近似值。MCMC背後的思想是，隨着生成更多的樣本，近似值會愈來愈接近實際的分佈。

MCMC方法有兩個部分，蒙特卡羅（Monte Carlo）是指使用重複隨機樣原本得到數值結果的通用技術。蒙特卡羅能夠被認爲是進行不少次的實驗，每次改變模型中的變量並觀察反應。經過選擇隨機值，咱們能夠研究參數空間的一大部分，可能的變量值的範圍。參數空間，以下圖所示：

 

很顯然，咱們不能嘗試圖中的每一個點，可是經過從較高几率（紅色）的區域隨機抽樣，咱們能夠建立最可能的模型。

馬爾可夫鏈（Markov Chain）

 

馬爾可夫鏈是當前狀態被下一個狀態所依賴的過程。馬爾可夫鏈是不記錄狀態的，由於只有當前的狀態纔是重要的，而它如何到達那個狀態並不重要。若是這有點難以理解，那麼考慮一個平常現象，天氣。若是咱們想預測明天的天氣，只能根據今天的天氣來進行一個合理的預測。若是今天下雪了，咱們查看一下下雪後次日的天氣分佈的歷史數據，以預測明天是什麼天氣的機率。馬爾可夫鏈的概念是，咱們不須要必定知道一個過程的整個歷史數據來預測下一個輸出，類似的現象在許多現實狀況中都能很好地預測。

MCMC結合了馬爾可夫鏈和蒙特卡羅的思想，是一種基於當前值對分佈的參數重複抽取隨機值的方法。每一個值的樣本都是隨機的，可是這些值的選擇受到當前狀態和參數的假定先驗分佈的限制。MCMC能夠看做是逐漸收斂到實際分佈的隨機遊動。

爲了抽取α和β的隨機值，咱們須要假設這些值的先驗分佈。因爲咱們沒有預先給參數假設，因此可使用一個正態分佈。正態分佈或高斯分佈由平均值定義，代表了數據的位置、方差和分佈範圍。下圖是具備不一樣平均值和範圍的幾個正態分佈:

 

咱們使用的特定MCMC算法稱爲Metropolis Hastings。爲了將觀測數據與模型聯繫起來，每次繪製一組隨機值時，該算法都會根據數據進行評估。若是這些隨機值與數據不一致，則會被拒絕，而且模型保持在當前狀態。反之，若是與數據一致，則將這些隨機值分配給參數並變成當前狀態。這個過程持續必定數量的迭代，那麼模型的精確度就會隨着迭代數量的增長而提升。

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

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

（2）根據當前狀態給α和β隨機分配新的值；

（3）檢查新的隨機值是否符合觀測值。若是不符合，則拒絕這些隨機值，並返回到先前狀態，反之，要是符合，則接受，並更新爲新的當前狀態；

（4）如需迭代，則重複步驟2和3；

算法返回其生成的全部α和β的值。而後，咱們可使用這些值的平均值做爲邏輯函數中最有可能的α和β最終值。MCMC不能返回「True」值，而是返回一個分佈的近似值。由已有數據獲得的睡眠機率，其最終模型是具備α和β平均值的邏輯函數。

Python的實現

爲了在Python中實現MCMC，咱們將使用貝葉斯推理庫PyMC3，它對大部分細節進行了抽象，使咱們可以建立模型而不是空有理論。

下面的代碼建立了具備參數α和β、機率p和觀察值的完整模型。步驟變量引用特定的算法，而且變量sleep_trace保存了由模型生成的全部參數值。

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

爲了更好地瞭解代碼的運行，咱們能夠看下全部的模型運行過程當中所產生的α和β值。

 

上圖被稱爲軌跡圖。咱們能夠看到，每一個狀態都與先前的馬爾可夫鏈相關，但其值在蒙特卡羅抽樣中振盪明顯。

在MCMC中，很常見的作法是丟棄多達90%的軌跡。算法不當即收斂到真實分佈，並且初始值每每也不許確。以後的參數值一般會更好一些，這意味着應該用它們來構建模型。咱們使用了10000個樣本，並丟棄了前面的50%，可是在企業級應用中可能會使用成百上千個或數百萬個樣本。

MCMC收斂到真實值，但評估收斂多是很困難的。若是咱們想要最精確的結果，這是一個重要的因素。PyMC3庫已經建立了用於評估模型質量的函數，包括軌跡圖和自相關圖。

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

 

軌跡圖(上) 和自相關圖(下)

睡眠模型

在最終創建和運行模型以後，如今應該使用了。咱們將最後5000個α和β樣本的平均值做爲最有可能的參數值，這容許咱們建立單個曲線圖來模擬指定時間以後的入睡機率:

 

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

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

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

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

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

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

這些值是給定數據的最可能的估計值。然而，會有這些機率相關的不肯定性，由於模型是近似的。爲了表示這種不肯定性，咱們能夠在一個給定的時間點使用全部α和β的樣原本進行入睡機率預測，而後根據結果繪製柱狀：

 

 

上述結果給出了一個更好的MCMC模型能實際作到的指標。這個方法並無找到一個正確答案，而是可能值的一個樣本。貝葉斯推理是有實際用處的，由於它以機率的形式表示預測的結果。咱們能夠說獲得一個最有可能的答案，可是更準確的說法應該是任何預測都是一系列值的範圍。

睡醒模型

能夠用我早上睡醒的相關數據來找到一個相似的模型。我一般在早上6點起牀，下圖根據觀察結果顯示了從睡覺到睡醒的最終模型。

 

能夠經過模型來發現我在某一指定時間入睡的機率和最有可能睡醒的時間。

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

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

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

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

入睡持續時間

最後一個我想建立的是入睡時間模型。首先，咱們須要找到一個函數來模擬數據的分佈，但只能經過檢查數據來找到。

 

一個普通的分佈便可實現，但它不會獲取右邊的離羣點。可使用兩個相互獨立的正態分佈來表示這兩種模式，可是，我會使用一個偏態分佈。偏態分佈有三個參數，平均值、方差、誤差，而且這三個參數必須從MCMC算法中進行學習。下面的代碼建立模型並實現Metropolis Hastings抽樣：

withpm.Model() asduration_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)

如今，咱們可使用三個參數的平均值來構造最有可能的分佈。下圖表示數據的最終偏態分佈:

 

上圖看起來很合乎實際狀況，經過查詢模型能夠發現我至少得到必定的入睡持續時間，和最可能的入睡時間範圍的機率：

入睡至少持續6.5小時的機率= 99.16%.

入睡至少持續8小時的機率= 44.53%.

入睡至少持續9小時的機率= 10.94%.

入睡最可能持續的時間是7.67小時

結論

經過完成這個項目，告訴咱們解決問題的重要性是最好解決實際生活中的問題。數據科學須要不斷地進行知識積累，最有效的方法就是有了問題儘快開始。

原文連接