怎么用Python實(shí)現(xiàn)MCMC模型

這篇文章主要介紹“怎么用Python實(shí)現(xiàn)MCMC模型”，在日常操作中，相信很多人在怎么用Python實(shí)現(xiàn)MCMC模型問(wèn)題上存在疑惑，小編查閱了各式資料，整理出簡(jiǎn)單好用的操作方法，希望對(duì)大家解答”怎么用Python實(shí)現(xiàn)MCMC模型”的疑惑有所幫助！接下來(lái)，請(qǐng)跟著小編一起來(lái)學(xué)習(xí)吧！

介紹

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

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

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

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

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

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

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

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

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

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

馬爾可夫鏈（Markov Chain）

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

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

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

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

綜上所述，用MCMC解決問(wèn)題的基本步驟如下:

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

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

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

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

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

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

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

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

withpm.Model() as sleep_model:# Create the alpha and beta parameters# Assume a normal distributionalpha=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 functionp=pm.Deterministic('p', 1. / (1. +tt.exp(beta * time + alpha)))# Create the bernoulli parameter which uses observed data to inform the algorithmobserved=pm.Bernoulli('obs', p, observed=sleep_obs)# Using Metropolis Hastings Samplingstep=pm.Metropolis()# Draw the specified number of samplessleep_trace=pm.sample(N_SAMPLES, step=step);

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

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

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

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

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

軌跡圖(上) 和自相關(guān)圖(下)

睡眠模型

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

該模型能很好地表示數(shù)據(jù)，還獲取了我入睡模式固有的變化，它不是給出一個(gè)結(jié)果，而是給出一個(gè)概率。例如，可以查詢?cè)撃Ｐ鸵哉页鑫以谥付〞r(shí)間入睡的概率，并找到概率超過(guò)50%的時(shí)間點(diǎn):

晚上9點(diǎn)半入睡概率: 4.80%.

晚上9點(diǎn)半入睡概率: 27.44%.

晚上10點(diǎn)入睡概率: 73.91%.

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

可以看到我入睡的平均時(shí)間是晚上10點(diǎn)14分左右。

這些值是給定數(shù)據(jù)的最可能的估計(jì)值。然而，會(huì)有這些概率相關(guān)的不確定性，因?yàn)槟Ｐ褪墙频摹榱吮硎具@種不確定性，我們可以在一個(gè)給定的時(shí)間點(diǎn)使用所有α和β的樣本來(lái)進(jìn)行入睡概率預(yù)測(cè)，然后根據(jù)結(jié)果繪制柱狀：

上述結(jié)果給出了一個(gè)更好的MCMC模型能實(shí)際做到的指標(biāo)。這個(gè)方法并沒(méi)有找到一個(gè)正確答案，而是可能值的一個(gè)樣本。貝葉斯推理是有實(shí)際用處的，因?yàn)樗愿怕实男问奖硎绢A(yù)測(cè)的結(jié)果。我們可以說(shuō)得到一個(gè)最有可能的答案，但是更準(zhǔn)確的說(shuō)法應(yīng)該是任何預(yù)測(cè)都是一系列值的范圍。

睡醒模型

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

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

早上5點(diǎn)半睡醒的概率: 14.10%.

早上6點(diǎn)睡醒的概率: 37.94%.

早上6點(diǎn)半睡醒的概率: 69.49%.

在早上6點(diǎn)11分睡醒的概率超過(guò)50%。

入睡持續(xù)時(shí)間

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

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

withpm.Model() asduration_model:# Three parameters to samplealpha_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 variableduration_ =pm.SkewNormal('duration', alpha=alpha_skew, mu= mu_, sd=1/tau_, observed= duration)# Metropolis Hastings for samplingstep=pm.Metropolis()duration_trace=pm.sample(N_SAMPLES, step=step)

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

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

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

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

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

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

到此，關(guān)于“怎么用Python實(shí)現(xiàn)MCMC模型”的學(xué)習(xí)就結(jié)束了，希望能夠解決大家的疑惑。理論與實(shí)踐的搭配能更好的幫助大家學(xué)習(xí)，快去試試吧！若想繼續(xù)學(xué)習(xí)更多相關(guān)知識(shí)，請(qǐng)繼續(xù)關(guān)注億速云網(wǎng)站，小編會(huì)繼續(xù)努力為大家?guī)?lái)更多實(shí)用的文章！