|
編譯：伯樂(lè)在線(xiàn) - JLee 如有好文章投稿，請(qǐng)點(diǎn)擊 → 這里了解詳情
本系列：
第1節(jié)：估計(jì)模型參數(shù)
在這一節(jié)，我們將討論貝葉斯方法是如何思考數(shù)據(jù)的，我們?cè)鯓油ㄟ^(guò) MCMC 技術(shù)估計(jì)模型參數(shù)。
from IPython.display import Image import matplotlib.pyplot as plt import numpy as np import pandas as pd import pymc3 as pm import scipy import scipy.stats as stats import scipy.optimize as opt import statsmodels.api as sm %matplotlib inline plt.style.use('bmh') colors = ['#348ABD', '#A60628', '#7A68A6', '#467821', '#D55E00', '#CC79A7', '#56B4E9', '#009E73', '#F0E442', '#0072B2'] messages = pd.read_csv('data/hangout_chat_data.csv')
貝葉斯方法如何思考數(shù)據(jù)？
當(dāng)我開(kāi)始學(xué)習(xí)如何運(yùn)用貝葉斯方法的時(shí)候，我發(fā)現(xiàn)理解貝葉斯方法如何思考數(shù)據(jù)是很有用的。設(shè)想下述的場(chǎng)景：
一個(gè)好奇的男孩每天觀(guān)察從他家門(mén)前經(jīng)過(guò)的汽車(chē)的數(shù)量。他每天努力地記錄汽車(chē)的總數(shù)。一個(gè)星期過(guò)去后，他的筆記本上記錄著下面的數(shù)字：12,33,20,29,20,30,18
從貝葉斯方法的角度看，這個(gè)數(shù)據(jù)是由隨機(jī)過(guò)程產(chǎn)生的。但是，既然數(shù)據(jù)被觀(guān)測(cè)，它便固定了并且不會(huì)改變。這個(gè)隨機(jī)過(guò)程有些模型參數(shù)被固定了。然而，貝葉斯方法用概率分布來(lái)表示這些模型參數(shù)的不確定性。
由于這個(gè)男孩調(diào)查的是計(jì)數(shù)（非負(fù)整數(shù)），一個(gè)通常的做法是用泊松分布對(duì)數(shù)據(jù)（如隨機(jī)過(guò)程）建模。泊松分布只有一個(gè)參數(shù) μ，它既是數(shù)據(jù)的平均數(shù)，也是方差。下面是三個(gè)不同 μ 值的泊松分布。
代碼：
fig = plt.figure(figsize=(11,3)) ax = fig.add_subplot(111) x_lim = 60 mu = [5, 20, 40] for i in np.arange(x_lim): plt.bar(i, stats.poisson.pmf(mu[0], i), color=colors[3]) plt.bar(i, stats.poisson.pmf(mu[1], i), color=colors[4]) plt.bar(i, stats.poisson.pmf(mu[2], i), color=colors[5]) _ = ax.set_xlim(0, x_lim) _ = ax.set_ylim(0, 0.2) _ = ax.set_ylabel('Probability mass') _ = ax.set_title('Poisson distribution') _ = plt.legend(['$mu$ = %s' % mu[0], '$mu$ = %s' % mu[1], '$mu$ = %s' % mu[2]])
在上一節(jié)中，我們引入我的 hangout 聊天數(shù)據(jù)集。特別地，我對(duì)我的回復(fù)時(shí)間（response_time）感興趣。鑒于 response_time 是計(jì)數(shù)數(shù)據(jù)，我們可以用泊松分布對(duì)其建模，并估計(jì)參數(shù) μ 。我們將用頻率論方法和貝葉斯方法兩種方法來(lái)估計(jì)。
fig = plt.figure(figsize=(11,3)) _ = plt.title('Frequency of messages by response time') _ = plt.xlabel('Response time (seconds)') _ = plt.ylabel('Number of messages') _ = plt.hist(messages['time_delay_seconds'].values, range=[0, 60], bins=60, histtype='stepfilled')
頻率論方法估計(jì)μ
在進(jìn)入貝葉斯方法之前，讓我們先看一下頻率論方法估計(jì)泊松分布參數(shù)。我們將使用優(yōu)化技術(shù)使似然函數(shù)最大。
下面的函數(shù)poisson_logprob()返回在給定泊松分布模型和參數(shù)值的條件下，觀(guān)測(cè)數(shù)據(jù)總體的可能性。用方法opt.minimize_scalar找到在觀(guān)測(cè)數(shù)據(jù)基礎(chǔ)上參數(shù)值 μ 的最可信值（最大似然）。該方法的機(jī)理是，這個(gè)優(yōu)化技術(shù)會(huì)自動(dòng)迭代可能的mu值直到找到可能性最大的值。
y_obs = messages['time_delay_seconds'].values def poisson_logprob(mu, sign=-1): return np.sum(sign*stats.poisson.logpmf(y_obs, mu=mu)) freq_results = opt.minimize_scalar(poisson_logprob) %time print('The estimated value of mu is: %s' % freq_results['x'])
The estimated value of mu is: 18.0413533867 CPU times: user 63 μs, sys: 1e+03 ns, total: 64 μs Wall time: 67.9 μs
所以，μ 的估計(jì)值是18.0413533867。優(yōu)化技術(shù)沒(méi)有對(duì)不確定度進(jìn)行評(píng)估，它只返回一個(gè)點(diǎn)，效率很高。
下圖描述的是我們優(yōu)化的函數(shù)。對(duì)于每個(gè)μ值，圖線(xiàn)顯示給定數(shù)據(jù)和模型在μ處的似然度。優(yōu)化器以登山模式工作——從曲線(xiàn)上隨機(jī)一點(diǎn)開(kāi)始，不停向上攀登直到達(dá)到最高點(diǎn)。
x = np.linspace(1, 60) y_min = np.min([poisson_logprob(i, sign=1) for i in x]) y_max = np.max([poisson_logprob(i, sign=1) for i in x]) fig = plt.figure(figsize=(6,4)) _ = plt.plot(x, [poisson_logprob(i, sign=1) for i in x]) _ = plt.fill_between(x, [poisson_logprob(i, sign=1) for i in x], y_min, color=colors[0], alpha=0.3) _ = plt.title('Optimization of $mu$') _ = plt.xlabel('$mu$') _ = plt.ylabel('Log probability of $mu$ given data') _ = plt.vlines(freq_results['x'], y_max, y_min, colors='red', linestyles='dashed') _ = plt.scatter(freq_results['x'], y_max, s=110, c='red', zorder=3) _ = plt.ylim(ymin=y_min, ymax=0) _ = plt.xlim(xmin=1, xmax=60)
上述優(yōu)化方法估計(jì)泊松分布的參數(shù)值為 18 .我們知道，對(duì)任意一個(gè)泊松分布，參數(shù) μ 代表的是均值和方差。下圖描述了這個(gè)分布。
fig = plt.figure(figsize=(11,3)) ax = fig.add_subplot(111) x_lim = 60 mu = np.int(freq_results['x']) for i in np.arange(x_lim): plt.bar(i, stats.poisson.pmf(mu, i), color=colors[3]) _ = ax.set_xlim(0, x_lim) _ = ax.set_ylim(0, 0.1) _ = ax.set_xlabel('Response time in seconds') _ = ax.set_ylabel('Probability mass') _ = ax.set_title('Estimated Poisson distribution for Hangout chat response time') _ = plt.legend(['$lambda$ = %s' % mu])
上述泊松分布模型和 μ 的估計(jì)表明，觀(guān)測(cè)小于10 或大于 30 的可能性很小，絕大多數(shù)的概率分布在 10 和 30 之間。但是，這不能反映我們觀(guān)測(cè)到的介于 0 到 60 之間的數(shù)據(jù)。
貝葉斯方法估計(jì) μ
如果你之前遇到過(guò)貝葉斯理論，下面的公式會(huì)看起來(lái)很熟悉。我從沒(méi)認(rèn)同過(guò)這個(gè)框架直到我讀了John K. Kruschke的書(shū)“貝葉斯數(shù)據(jù)分析”，從他優(yōu)美簡(jiǎn)潔的貝葉斯圖表模型視角認(rèn)識(shí)這個(gè)公式。
代碼：
Image('graphics/Poisson-dag.png', width=320)
上述模式可以解讀如下（從下到上）：
對(duì)每一個(gè)對(duì)話(huà)（i）觀(guān)測(cè)計(jì)數(shù)數(shù)據(jù)（y） 數(shù)據(jù)由一個(gè)隨機(jī)過(guò)程產(chǎn)生，我們認(rèn)為可以用泊松分布模擬 泊松分布只有一個(gè)參數(shù)，介于 0 到 60 之間
由于我們沒(méi)有任何依據(jù)對(duì)μ在這個(gè)范圍內(nèi)進(jìn)行預(yù)估，就用均勻分布對(duì) μ 建模
神奇的方法MCMC
下面的動(dòng)畫(huà)很好的演示了馬爾科夫鏈蒙特卡羅方法（MCMC）的過(guò)程。MCMC采樣器從先驗(yàn)分布中選取參數(shù)值，計(jì)算從這些參數(shù)值決定的某個(gè)分布中得到觀(guān)測(cè)數(shù)據(jù)的可能性。
這個(gè)計(jì)算可以作為MCMC采樣器的導(dǎo)引。從參數(shù)的先驗(yàn)分布中選取值，然后計(jì)算給定數(shù)據(jù)條件下這些參數(shù)值可能性——導(dǎo)引采樣器到更高概率的區(qū)域。
與上面討論的頻率論優(yōu)化技術(shù)在概念上有相似之處，MCMC采樣器向可能性最高的區(qū)域運(yùn)動(dòng)。但是，貝葉斯方法不關(guān)心尋找絕對(duì)最大值，而是遍歷和收集概率最高區(qū)域附近的樣本。所有收集到的樣本都可認(rèn)為是一個(gè)可信的參數(shù)。
Image(url='graphics/mcmc-animate.gif')
with pm.Model() as model: mu = pm.Uniform('mu', lower=0, upper=60) likelihood = pm.Poisson('likelihood', mu=mu, observed=messages['time_delay_seconds'].values) start = pm.find_MAP() step = pm.Metropolis() trace = pm.sample(200000, step, start=start, progressbar=True)
Applied interval-transform to mu and added transformed mu_interval to model. [-----------------100%-----------------] 200000 of 200000 complete in 48.3 sec
上面的代碼通過(guò)遍歷 μ 的后驗(yàn)分布的高概率區(qū)域，收集了 200,000 個(gè) μ 的可信樣本。下圖（左）顯示這些值的分布，平均值與頻率論的估值（紅線(xiàn)）幾乎一樣。但是，我們同時(shí)也得到了不確定度，μ 值介于 17 到 19 之間都是可信的。我們稍后會(huì)看到這個(gè)不確定度很有價(jià)值。
_ = pm.traceplot(trace, varnames=['mu'], lines={'mu': freq_results['x']})
丟棄早期樣本（壞點(diǎn)）
你可能疑惑上面 MCMC 代碼中pm.find.MAP()的作用。MAP 代表最大后驗(yàn)估計(jì)，它幫助 MCMC 采樣器尋找合適的采樣起始點(diǎn)。理想情況下，模型從高概率區(qū)域開(kāi)始，但有時(shí)不是這樣。因此，樣本集中的早期樣本（壞點(diǎn)）經(jīng)常被丟棄。
fig = plt.figure(figsize=(11,3)) plt.subplot(121) _ = plt.title('Burnin trace') _ = plt.ylim(ymin=16.5, ymax=19.5) _ = plt.plot(trace.get_values('mu')[:1000]) fig = plt.subplot(122) _ = plt.title('Full trace') _ = plt.ylim(ymin=16.5, ymax=19.5) _ = plt.plot(trace.get_values('mu'))
模型的收斂性
樣本軌跡
由于上面模型對(duì) μ 的估計(jì)不能表明估計(jì)的效果很好,可以采用一些檢驗(yàn)手段。首先，觀(guān)察樣本輸出，樣本的軌跡應(yīng)該輕微上下浮動(dòng)，就像一條毛蟲(chóng)。如果你看到軌跡上下跳躍或滯留在某點(diǎn)，這就是模型不收斂的標(biāo)志，通過(guò) MCMC 采樣器得到的估計(jì)沒(méi)有意義。
自相關(guān)圖
也可采另一種測(cè)試，自相關(guān)測(cè)試（見(jiàn)下圖），這是評(píng)估MCMC采樣鏈中樣本前后相關(guān)關(guān)系的一種方法。相關(guān)性弱的樣本比相關(guān)性強(qiáng)的樣本能夠給參數(shù)估計(jì)提供更多的信息。
直觀(guān)上看，自相關(guān)圖應(yīng)該快速衰減到0附近，然后在 0 附近上下波動(dòng)。如果你的自相關(guān)圖沒(méi)有衰減，通常這是弱融合的標(biāo)志，你需要重新審查你的模型選擇（如似然度）和抽樣方法（如Metropolis）
_ = pm.autocorrplot(trace[:2000], varnames=['mu'])
看完本文有收獲？請(qǐng)轉(zhuǎn)發(fā)分享給更多人 關(guān)注「Python開(kāi)發(fā)者」，提升Python技能
|
