本系列：

第1節：估計模型參數

在這一節，我們將討論貝葉斯方法是如何思考數據的，我們怎樣通過 MCMC 技術估計模型參數。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

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

貝葉斯方法如何思考數據？

當我開始學習如何運用貝葉斯方法的時候，我發現理解貝葉斯方法如何思考數據是很有用的。設想下述的場景：

一個好奇的男孩每天觀察從他家門前經過的汽車的數量。他每天努力地記錄汽車的總數。一個星期過去后，他的筆記本上記錄著下面的數字：12,33,20,29,20,30,18

從貝葉斯方法的角度看，這個數據是由隨機過程產生的。但是，既然數據被觀測，它便固定了并且不會改變。這個隨機過程有些模型參數被固定了。然而，貝葉斯方法用概率分布來表示這些模型參數的不確定性。

由于這個男孩調查的是計數(非負整數)，一個通常的做法是用泊松分布對數據(如隨機過程)建模。泊松分布只有一個參數?μ，它既是數據的平均數，也是方差。下面是三個不同?μ 值的泊松分布。

代碼：

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

在上一節中，我們引入我的 hangout 聊天數據集。特別地，我對我的回復時間(response_time)感興趣。鑒于 response_time 是計數數據，我們可以用泊松分布對其建模，并估計參數 μ 。我們將用頻率論方法和貝葉斯方法兩種方法來估計。

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

頻率論方法估計μ

在進入貝葉斯方法之前，讓我們先看一下頻率論方法估計泊松分布參數。我們將使用優化技術使似然函數最大。

下面的函數poisson_logprob()返回在給定泊松分布模型和參數值的條件下，觀測數據總體的可能性。用方法opt.minimize_scalar找到在觀測數據基礎上參數值 μ 的最可信值(最大似然)。該方法的機理是，這個優化技術會自動迭代可能的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。優化技術沒有對不確定度進行評估，它只返回一個點，效率很高。

下圖描述的是我們優化的函數。對于每個μ值，圖線顯示給定數據和模型在μ處的似然度。優化器以登山模式工作——從曲線上隨機一點開始，不停向上攀登直到達到最高點。

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)

上述優化方法估計泊松分布的參數值為 ?18 .我們知道，對任意一個泊松分布，參數 μ 代表的是均值和方差。下圖描述了這個分布。

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 的可能性很小，絕大多數的概率分布在 10 和 30 之間。但是，這不能反映我們觀測到的介于 0 到 60 之間的數據。

貝葉斯方法估計?μ

如果你之前遇到過貝葉斯理論，下面的公式會看起來很熟悉。我從沒認同過這個框架直到我讀了John K. Kruschke的書“貝葉斯數據分析”，從他優美簡潔的貝葉斯圖表模型視角認識這個公式。

代碼：

Image('graphics/Poisson-dag.png', width=320)

上述模式可以解讀如下(從下到上)：

對每一個對話(i)觀測計數數據(y)

數據由一個隨機過程產生，我們認為可以用泊松分布模擬

泊松分布只有一個參數，介于 0 到 60 之間

由于我們沒有任何依據對μ在這個范圍內進行預估，就用均勻分布對 μ 建模

神奇的方法MCMC

下面的動畫很好的演示了馬爾科夫鏈蒙特卡羅方法(MCMC)的過程。MCMC采樣器從先驗分布中選取參數值，計算從這些參數值決定的某個分布中得到觀測數據的可能性。

這個計算可以作為MCMC采樣器的導引。從參數的先驗分布中選取值，然后計算給定數據條件下這些參數值可能性——導引采樣器到更高概率的區域。

與上面討論的頻率論優化技術在概念上有相似之處，MCMC采樣器向可能性最高的區域運動。但是，貝葉斯方法不關心尋找絕對最大值，而是遍歷和收集概率最高區域附近的樣本。所有收集到的樣本都可認為是一個可信的參數。

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

上面的代碼通過遍歷 μ 的后驗分布的高概率區域，收集了 200,000 個 μ 的可信樣本。下圖(左)顯示這些值的分布，平均值與頻率論的估值(紅線)幾乎一樣。但是，我們同時也得到了不確定度，μ 值介于 17 到 19 之間都是可信的。我們稍后會看到這個不確定度很有價值。

_ = pm.traceplot(trace, varnames=['mu'], lines={'mu': freq_results['x']})

丟棄早期樣本(壞點)

你可能疑惑上面 MCMC 代碼中pm.find.MAP()的作用。MAP 代表最大后驗估計，它幫助 MCMC 采樣器尋找合適的采樣起始點。理想情況下，模型從高概率區域開始，但有時不是這樣。因此，樣本集中的早期樣本(壞點)經常被丟棄。

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采樣鏈中樣本前后相關關系的一種方法。相關性弱的樣本比相關性強的樣本能夠給參數估計提供更多的信息。

直觀上看，自相關圖應該快速衰減到0附近，然后在 0 附近上下波動。如果你的自相關圖沒有衰減，通常這是弱融合的標志，你需要重新審查你的模型選擇(如似然度)和抽樣方法(如Metropolis)

_ = pm.autocorrplot(trace[:2000], varnames=['mu'])

打賞支持我翻譯更多好文章，謝謝！

打賞譯者

打賞支持我翻譯更多好文章，謝謝！

任選一種支付方式

總結

以上是生活随笔為你收集整理的如何用python进行建模_用 Python 进行贝叶斯模型建模（1）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。