|
編譯：伯樂在線 - JLee 如有好文章投稿，請點擊 → 這里了解詳情
本系列：
第1節(jié)：估計模型參數(shù)
在這一節(jié)，我們將討論貝葉斯方法是如何思考數(shù)據的，我們怎樣通過 MCMC 技術估計模型參數(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ù)據？
當我開始學習如何運用貝葉斯方法的時候，我發(fā)現(xiàn)理解貝葉斯方法如何思考數(shù)據是很有用的。設想下述的場景：
一個好奇的男孩每天觀察從他家門前經過的汽車的數(shù)量。他每天努力地記錄汽車的總數(shù)。一個星期過去后，他的筆記本上記錄著下面的數(shù)字：12,33,20,29,20,30,18
從貝葉斯方法的角度看，這個數(shù)據是由隨機過程產生的。但是，既然數(shù)據被觀測，它便固定了并且不會改變。這個隨機過程有些模型參數(shù)被固定了。然而，貝葉斯方法用概率分布來表示這些模型參數(shù)的不確定性。
由于這個男孩調查的是計數(shù)（非負整數(shù)），一個通常的做法是用泊松分布對數(shù)據（如隨機過程）建模。泊松分布只有一個參數(shù) μ，它既是數(shù)據的平均數(shù)，也是方差。下面是三個不同 μ 值的泊松分布。
代碼：
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ù)據集。特別地，我對我的回復時間（response_time）感興趣。鑒于 response_time 是計數(shù)數(shù)據，我們可以用泊松分布對其建模，并估計參數(shù) μ 。我們將用頻率論方法和貝葉斯方法兩種方法來估計。
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')
頻率論方法估計μ
在進入貝葉斯方法之前，讓我們先看一下頻率論方法估計泊松分布參數(shù)。我們將使用優(yōu)化技術使似然函數(shù)最大。
下面的函數(shù)poisson_logprob()返回在給定泊松分布模型和參數(shù)值的條件下，觀測數(shù)據總體的可能性。用方法opt.minimize_scalar找到在觀測數(shù)據基礎上參數(shù)值 μ 的最可信值（最大似然）。該方法的機理是，這個優(yōu)化技術會自動迭代可能的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
所以，μ 的估計值是18.0413533867。優(yōu)化技術沒有對不確定度進行評估，它只返回一個點，效率很高。
下圖描述的是我們優(yōu)化的函數(shù)。對于每個μ值，圖線顯示給定數(shù)據和模型在μ處的似然度。優(yōu)化器以登山模式工作——從曲線上隨機一點開始，不停向上攀登直到達到最高點。
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)化方法估計泊松分布的參數(shù)值為 18 .我們知道，對任意一個泊松分布，參數(shù) μ 代表的是均值和方差。下圖描述了這個分布。
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])
上述泊松分布模型和 μ 的估計表明，觀測小于10 或大于 30 的可能性很小，絕大多數(shù)的概率分布在 10 和 30 之間。但是，這不能反映我們觀測到的介于 0 到 60 之間的數(shù)據。
貝葉斯方法估計 μ
如果你之前遇到過貝葉斯理論，下面的公式會看起來很熟悉。我從沒認同過這個框架直到我讀了John K. Kruschke的書“貝葉斯數(shù)據分析”，從他優(yōu)美簡潔的貝葉斯圖表模型視角認識這個公式。
代碼：
Image('graphics/Poisson-dag.png', width=320)
上述模式可以解讀如下（從下到上）：
對每一個對話（i）觀測計數(shù)數(shù)據（y） 數(shù)據由一個隨機過程產生，我們認為可以用泊松分布模擬 泊松分布只有一個參數(shù)，介于 0 到 60 之間
由于我們沒有任何依據對μ在這個范圍內進行預估，就用均勻分布對 μ 建模
神奇的方法MCMC
下面的動畫很好的演示了馬爾科夫鏈蒙特卡羅方法（MCMC）的過程。MCMC采樣器從先驗分布中選取參數(shù)值，計算從這些參數(shù)值決定的某個分布中得到觀測數(shù)據的可能性。
這個計算可以作為MCMC采樣器的導引。從參數(shù)的先驗分布中選取值，然后計算給定數(shù)據條件下這些參數(shù)值可能性——導引采樣器到更高概率的區(qū)域。
與上面討論的頻率論優(yōu)化技術在概念上有相似之處，MCMC采樣器向可能性最高的區(qū)域運動。但是，貝葉斯方法不關心尋找絕對最大值，而是遍歷和收集概率最高區(qū)域附近的樣本。所有收集到的樣本都可認為是一個可信的參數(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
上面的代碼通過遍歷 μ 的后驗分布的高概率區(qū)域，收集了 200,000 個 μ 的可信樣本。下圖（左）顯示這些值的分布，平均值與頻率論的估值（紅線）幾乎一樣。但是，我們同時也得到了不確定度，μ 值介于 17 到 19 之間都是可信的。我們稍后會看到這個不確定度很有價值。
_ = pm.traceplot(trace, varnames=['mu'], lines={'mu': freq_results['x']})
丟棄早期樣本（壞點）
你可能疑惑上面 MCMC 代碼中pm.find.MAP()的作用。MAP 代表最大后驗估計，它幫助 MCMC 采樣器尋找合適的采樣起始點。理想情況下，模型從高概率區(qū)域開始，但有時不是這樣。因此，樣本集中的早期樣本（壞點）經常被丟棄。
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'))
模型的收斂性
樣本軌跡
由于上面模型對 μ 的估計不能表明估計的效果很好,可以采用一些檢驗手段。首先，觀察樣本輸出，樣本的軌跡應該輕微上下浮動，就像一條毛蟲。如果你看到軌跡上下跳躍或滯留在某點，這就是模型不收斂的標志，通過 MCMC 采樣器得到的估計沒有意義。
自相關圖
也可采另一種測試，自相關測試（見下圖），這是評估MCMC采樣鏈中樣本前后相關關系的一種方法。相關性弱的樣本比相關性強的樣本能夠給參數(shù)估計提供更多的信息。
直觀上看，自相關圖應該快速衰減到0附近，然后在 0 附近上下波動。如果你的自相關圖沒有衰減，通常這是弱融合的標志，你需要重新審查你的模型選擇（如似然度）和抽樣方法（如Metropolis）
_ = pm.autocorrplot(trace[:2000], varnames=['mu'])
看完本文有收獲？請轉發(fā)分享給更多人 關注「Python開發(fā)者」，提升Python技能
|
