用 Python 进行贝叶斯模型建模（1）

506 查看

本系列：

第1节：估计模型参数

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

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值直到找到可能性最大的值。

贝叶斯方法如何思考数据？

当我开始学习如何运用贝叶斯方法的时候，我发现理解贝叶斯方法如何思考数据是很有用的。设想下述的场景：

一个好奇的男孩每天观察从他家门前经过的汽车的数量。他每天努力地记录汽车的总数。一个星期过去后，他的笔记本上记录着下面的数字：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]])

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

频率论方法估计μ

在进入贝叶斯方法之前，让我们先看一下频率论方法估计泊松分布参数。我们将使用优化技术使似然函数最大。

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

on-num" data-line="crayon-5812b1408f85c085651845-3">3

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。优化技术没有对不确定度进行评估，它只返回一个点，效率很高。

下图描述的是我们优化的函数。对于每个μ值，图线显示给定数据和模型在μ处的似然度。优化器以登山模式工作——从曲线上随机一点开始，不停向上攀登直到达到最高点。