如何用python做模型_用 Python 进行贝叶斯模型建模（1）

第1节：估计模型参数

在这一节，我们将讨论贝叶斯方法是如何思考数据的，我们怎样通过 MCMC 技术估计模型参数。

fromIPython.displayimportImage

importmatplotlib.pyplotasplt

importnumpyasnp

importpandasaspd

importpymc3aspm

importscipy

importscipy.statsasstats

importscipy.optimizeasopt

importstatsmodels.apiassm

%matplotlibinline

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]

foriinnp.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

defpoisson_logprob(mu,sign=-1):

returnnp.sum(sign*stats.poisson.logpmf(y_obs,mu=mu))

freq_results=opt.minimize_scalar(poisson_logprob)

%timeprint("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)foriinx])

y_max=np.max([poisson_logprob(i,sign=1)foriinx])

fig=plt.figure(figsize=(6,4))

_=plt.plot(x,[poisson_logprob(i,sign=1)foriinx])

_=plt.fill_between(x,[poisson_logprob(i,sign=1)foriinx],

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

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

withpm.Model()asmodel:

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)

Appliedinterval-transform to muandadded transformed mu_interval tomodel.

[-----------------100%-----------------]200000of200000completein48.3sec

上面的代码通过遍历 μ 的后验分布的高概率区域，收集了 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'])

编译：伯乐在线 - JLee

Python开发整理发布，转载请联系作者获得授权。