通過(guò)把馬爾科夫鏈蒙特卡羅（MCMC）應(yīng)用于一個(gè)具體問(wèn)題,，本文介紹了 Python 中 MCMC 的入門(mén)級(jí)應(yīng)用,。機(jī)器之心對(duì)本文進(jìn)行了編譯介紹。

過(guò)去幾月中,，我總是反復(fù)遇到同一個(gè)數(shù)據(jù)科學(xué)術(shù)語(yǔ)：馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo/MCMC）,。每當(dāng)我在實(shí)驗(yàn)室、博客,、文章中聽(tīng)到這個(gè)概念,，我常常點(diǎn)頭贊同,，覺(jué)得它很酷，但實(shí)際上并沒(méi)有一個(gè)清晰的認(rèn)知,。有幾次我嘗試著學(xué)習(xí) MCMC 和貝葉斯推理,，我每次從閱讀書(shū)籍開(kāi)始，結(jié)果卻很快放棄,。我感到很惱怒,，于是決定轉(zhuǎn)向一種學(xué)習(xí)任何新技能的最佳方法：將它應(yīng)用于一個(gè)具體問(wèn)題。

使用我的睡眠數(shù)據(jù)（我一直打算對(duì)此一探究竟）和一本實(shí)際應(yīng)用的書(shū)（Bayesian Methods for Hackers）,，我終于通過(guò)一個(gè)實(shí)際問(wèn)題學(xué)習(xí)了馬爾科夫鏈蒙特卡羅,。像往常一樣，比起閱讀抽象的概念,，將這些技術(shù)應(yīng)用到具體問(wèn)題中能讓學(xué)習(xí)變得更簡(jiǎn)單,、更愉快。本文介紹了 Python 中的馬爾科夫鏈蒙特卡羅的入門(mén)級(jí)應(yīng)用,，正是它教會(huì)了我使用這個(gè)強(qiáng)大的建模分析工具,。

我鼓勵(lì)大家參閱 GitHub，并將其用于自己的數(shù)據(jù)當(dāng)中,。本文將重點(diǎn)介紹它的應(yīng)用和結(jié)果,，所以會(huì)產(chǎn)生很多高層次的話題。如果你在閱讀后想了解更多,，可以閱讀文中提供的鏈接,。

簡(jiǎn)介

我的 Garmin Vivosmart 手表可以根據(jù)心率和運(yùn)動(dòng)情況追蹤我的睡眠和起床狀況。它并非 100％準(zhǔn)確,，不過(guò)真實(shí)數(shù)據(jù)從不完美,，我們?nèi)匀豢梢越柚_的模型從噪聲數(shù)據(jù)中提取有用的知識(shí)！

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

本項(xiàng)目的目標(biāo)是借助睡眠數(shù)據(jù)創(chuàng)建一個(gè)模型,，通過(guò)把睡眠看作時(shí)間函數(shù),，而確定睡眠的后驗(yàn)概率。由于時(shí)間是連續(xù)變量,，確定整個(gè)后驗(yàn)分布非常棘手,。因此我們轉(zhuǎn)而使用一些可實(shí)現(xiàn)近似分布的方法，比如馬爾可夫鏈蒙特卡羅（MCMC）,。

選擇一個(gè)概率分布

在開(kāi)始使用 MCMC 之前,，我們需要確定一個(gè)合適的函數(shù)來(lái)對(duì)睡眠的后驗(yàn)概率分布進(jìn)行建模。一個(gè)簡(jiǎn)單的方法是直觀檢查這些數(shù)據(jù),。對(duì)于我的睡眠的時(shí)間函數(shù)的觀察如下圖所示,。

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

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

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

其中,，β 和 α 是我們?cè)?MCMC 過(guò)程中必須學(xué)習(xí)的模型參數(shù),。具有不同參數(shù)的 logsitic 函數(shù)圖像如下所示。

logsitic 函數(shù)很適合本案例中的數(shù)據(jù),，因?yàn)槿胨目赡苄詴?huì)逐漸轉(zhuǎn)變,，此函數(shù)能捕捉睡眠模式之中的變化情況,。我們希望能夠在函數(shù)中插入時(shí)間 t,，獲得睡眠概率（其值在 0 和 1 之間）,。我們最終得到的不是在晚上 10:00 入睡與否的直接答案,，而是一個(gè)概率,。為了建立這個(gè)模型,，我們使用這些數(shù)據(jù),，通過(guò) MCMC 尋找最佳的 α 和 β 參數(shù),。

馬爾科夫鏈蒙特卡羅

馬爾可夫鏈蒙特卡羅指從概率分布中抽樣以構(gòu)建最大可能分布的一類方法,。我們不能直接構(gòu)建 logistic 分布，所以,，與之相反,，我們?yōu)楹瘮?shù)的參數(shù)（α 和 β）生成了上千個(gè)值——被稱為樣本——從而創(chuàng)造分布的近似值。MCMC 背后的思想是,，當(dāng)我們生成更多的樣本時(shí),，我們的近似值越來(lái)越接近實(shí)際的真實(shí)分布,。

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

顯然,，我們無(wú)法一一嘗試圖像中的每一個(gè)點(diǎn)。但是通過(guò)對(duì)較高概率區(qū)域（紅色區(qū)域）進(jìn)行隨機(jī)抽樣,，我們可以為問(wèn)題建立最可能的模型,。

馬爾科夫鏈

馬爾科夫鏈?zhǔn)且粋€(gè)隨機(jī)過(guò)程，其中次態(tài)僅依賴于當(dāng)前狀態(tài)（在此語(yǔ)境中,，一個(gè)狀態(tài)指的是參數(shù)的一次賦值）,。馬爾科夫鏈沒(méi)有記憶性，因?yàn)橹挥挟?dāng)前狀態(tài)對(duì)下一狀態(tài)起作用,，而與到達(dá)當(dāng)前狀態(tài)的方式無(wú)關(guān),。如果這種說(shuō)法還是有些難以理解，我們可以用日?，F(xiàn)象中的天氣來(lái)舉例,。如果我們想預(yù)測(cè)明日天氣，我們可以僅通過(guò)今日天氣來(lái)得到一個(gè)合理的估計(jì),。如果今天下雪了,，我們可以查看下雪次日天氣分布的歷史數(shù)據(jù)，估算明天天氣的概率,。馬爾科夫鏈的概念在于,，我們無(wú)需了解整個(gè)歷史過(guò)程就能預(yù)測(cè)下一狀態(tài)，這個(gè)近似在許多現(xiàn)實(shí)情況中就能很好地工作,。

綜合馬爾科夫鏈和蒙特卡羅的思想,，馬爾科夫鏈蒙特卡羅是一種基于當(dāng)前值重復(fù)繪制某一分布參數(shù)隨機(jī)值的方法。每個(gè)值的樣本都是隨機(jī)的,，但是值的選擇受限于當(dāng)前狀態(tài)和假定的參數(shù)先驗(yàn)分布,。MCMC 可以被認(rèn)為是一種隨機(jī)游走，在這個(gè)過(guò)程中逐漸收斂到真實(shí)分布,。

為了繪制 α 和 β 的隨機(jī)值，我們需要假設(shè)這些值的先驗(yàn)分布。由于我們對(duì)參數(shù)沒(méi)有任何提前的假設(shè),，我們可以使用正態(tài)分布,。正態(tài)分布也稱高斯分布，它由均值和方差定義,，分別顯示數(shù)據(jù)的位置以及擴(kuò)散情況,。下圖是具有不同均值和方差的幾種正態(tài)分布：

我們所使用的 MCMC 算法被稱為 Metropolis Hastings。為了將我們觀察的數(shù)據(jù)與模型聯(lián)系起來(lái),，每繪制一組隨機(jī)值,，算法會(huì)根據(jù)數(shù)據(jù)對(duì)其進(jìn)行評(píng)估。如果隨機(jī)值與數(shù)據(jù)不一致（這里稍微進(jìn)行了一些簡(jiǎn)化）,，這些值將被拒絕,，模型保持當(dāng)前狀態(tài)。反之,，如果隨機(jī)值與數(shù)據(jù)一致,，這些值將會(huì)分配給參數(shù)并成為當(dāng)前狀態(tài)。該過(guò)程將持續(xù)進(jìn)行指定的步驟數(shù)目,，模型的準(zhǔn)確率也隨著步驟數(shù)量的增加而改善,。

綜合而言，馬爾科夫鏈蒙特卡羅在我們的問(wèn)題當(dāng)中基本步驟如下：

為 logistic 函數(shù)選擇一組初始參數(shù) α 和 β,。

根據(jù)當(dāng)前狀態(tài),，把新的隨機(jī)值分配給 α 和 β。

檢查新的隨機(jī)值是否與觀察結(jié)果一致,。如果不一致,，拒絕這些隨機(jī)值并返回前一個(gè)狀態(tài)。如果一致,，則接受這些值,，將其作為新的當(dāng)前狀態(tài)。

對(duì)指定的迭代次數(shù)重復(fù)執(zhí)行步驟 2 和 3,。

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

Python 實(shí)現(xiàn)

上述細(xì)節(jié)在我腦海中徘徊已久，最后終于在 Python 中進(jìn)行了實(shí)現(xiàn)！親眼看到第一手的結(jié)果比讀取別人的描述有幫助得多,。要在 Python 中實(shí)現(xiàn) MCMC，我們需要使用 PyMC3 貝葉斯推理庫(kù),。它將大部分細(xì)節(jié)進(jìn)行了抽象,，從而讓我們能不迷失在理論中,，并建立我們的模型。

下面的代碼創(chuàng)建的模型帶有參數(shù) α 和 β,、概率 p 和觀察結(jié)果 observed。step 變量指的是特定的算法,，sleep_trace 則保存了模型生成的所有參數(shù)值,。

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

（請(qǐng)?jiān)?notebook 中查閱完整代碼）

為了了解運(yùn)行此代碼時(shí)發(fā)生的情況,，我們可以查看模型運(yùn)行過(guò)程中生成的 α 和 β 的所有值,。

它們被稱為軌跡圖。我們可以看到,，每個(gè)狀態(tài)都與之前的狀態(tài)有關(guān)（馬爾科夫鏈），但是這些值波動(dòng)顯著（蒙特卡羅采樣）,。

在 MCMC 中，通常高達(dá) 90% 的軌跡會(huì)被拋棄,。該算法無(wú)法立即收斂到真正的分布,，且初始值往往并不準(zhǔn)確,。后期的參數(shù)值通常更好，這意味著它們是適用于建立模型的參數(shù),。我們使用了 10000 個(gè)樣本并丟棄了前 50%，但是一個(gè)行業(yè)的應(yīng)用可能會(huì)使用數(shù)十萬(wàn)甚至上百萬(wàn)個(gè)樣本,。

給定足夠多的迭代次數(shù),，MCMC 將收斂于真實(shí)值,。但是，對(duì)收斂進(jìn)行評(píng)估可能比較困難,。對(duì)此我將不在本文討論（一個(gè)方法是測(cè)量軌跡的自相關(guān)），但是,，如果我們想要結(jié)果最準(zhǔn)確,，這是一個(gè)重要的考慮因素,。PyMC3 建立了評(píng)估模型好壞的函數(shù)，其中包括軌跡圖和自相關(guān)圖,。

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

軌跡圖（左）和自相關(guān)圖（右）

睡眠模型

最終建立并運(yùn)行模型之后，是時(shí)候使用結(jié)果了,。我們將最后 5000 個(gè) α 和 β 樣本的平均值作為參數(shù)最可能的值，這就讓我們能夠創(chuàng)建一條曲線,，建模睡眠后驗(yàn)概率：

該模型能很好地反映數(shù)據(jù)的結(jié)果。此外,，它捕捉了我睡眠模式當(dāng)中的固有變化,。該模型給出的不是一個(gè)簡(jiǎn)單的是非答案,，而是一個(gè)概率。例如,，我們可以通過(guò)該模型找到在給定時(shí)間我睡著的概率,，并能找到睡眠概率經(jīng)過(guò) 50% 的時(shí)間：

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 點(diǎn)上床睡覺(jué),，但這顯然不是大多數(shù)下的實(shí)際情況,。我們可以發(fā)現(xiàn)，我上床的平均時(shí)間是晚上 10:14 左右,。

在數(shù)據(jù)給定的情況下,，這些值是最有可能的估計(jì)值。然而,，因?yàn)槟Ｐ捅旧硎墙频?，所以存在與這些概率相關(guān)的不確定性。為了表示這種不確定性,，我們可以使用所有的 α 和 β 樣本（而不是它們的平均值）來(lái)預(yù)測(cè)某一給定時(shí)間的睡眠概率，然后據(jù)此繪制直方圖,。

這些結(jié)果更好地反映了 MCMC 模型真正做了什么,。MCMC 找到的不是一個(gè)簡(jiǎn)單的答案,，而是可能值的樣本。貝葉斯推理在現(xiàn)實(shí)世界中起到了重要作用,，是因?yàn)樗鼜母怕实慕嵌缺硎绢A(yù)測(cè)結(jié)果。我們可以說(shuō),，問(wèn)題會(huì)有一個(gè)可能性最大的答案,，但是更加準(zhǔn)確的回應(yīng)是任何預(yù)測(cè)都存在一系列的可能值,。

喚醒模型

我可以使用描述早晨醒來(lái)時(shí)間的數(shù)據(jù)建立一個(gè)類似的模型。我定了一個(gè)鬧鐘,，努力在早晨 6:00 起床。但是可以看到,，我并非每日都是如此,。下圖展現(xiàn)了我從入睡到醒來(lái)過(guò)渡過(guò)程的最終模型以及觀察數(shù)據(jù),。

通過(guò)查詢模型，我們可以找出在給定時(shí)間我睡著的概率以及最有可能醒來(lái)的時(shí)間,。

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.

看起來(lái)我得處理一下我的鬧鐘了！

睡眠時(shí)長(zhǎng)

出于好奇心和練習(xí)目的,，我最終想創(chuàng)造的是關(guān)于我睡眠時(shí)長(zhǎng)的模型,。首先，我們需要找到一個(gè)函數(shù)來(lái)模擬數(shù)據(jù)的分布,。我猜想結(jié)果應(yīng)該會(huì)是正態(tài)分布的形式，但是我們只有通過(guò)檢查數(shù)據(jù)才能得到最終結(jié)果,。

正態(tài)分布確實(shí)可行,，但這無(wú)法捕捉右側(cè)的偏離點(diǎn)（即我睡眠時(shí)間非常長(zhǎng)的情況）,。我們可以用兩個(gè)獨(dú)立的正態(tài)分布來(lái)表示兩個(gè)模型，但是,，我想使用偏正態(tài)分布,。偏正態(tài)分布有三個(gè)參數(shù)：均值,、方差、偏斜度 α,。以上三個(gè)參數(shù)都需要通過(guò) MCMC 來(lái)學(xué)習(xí),。下面的代碼建立了上述模型,，并進(jìn)行了 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)

現(xiàn)在,，我們可以使用這三個(gè)參數(shù)的均值來(lái)構(gòu)建最有可能的分布。以下是根據(jù)數(shù)據(jù)觀察得到的最終偏正態(tài)分布,。

看起來(lái)擬合得不錯(cuò)！我們可以通過(guò)查詢模型,，找到我至少可以獲得一定睡眠時(shí)長(zhǎng)的概率,，同時(shí)也能找到最可能的睡眠時(shí)長(zhǎng)：

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.

我對(duì)這個(gè)結(jié)果并不是完全滿意,，但是對(duì)一個(gè)研究生而言，這樣的結(jié)果已經(jīng)不錯(cuò)啦,。

小結(jié)

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