本系列:
在这一节,我们将讨论贝叶斯方法是如何思考数据的,我们怎样通过 MCMC 技术估计模型参数。
123456789101112131415161718 | from IPython.display import Image import matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport pymc3 as pmimport scipyimport scipy.stats as statsimport scipy.optimize as optimport 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
从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。
由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。
代码:
1234567891011121314 | fig = plt.figure(figsize=(11,3))ax = fig.add_subplot(111)x_lim = 60mu = [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
是计数数据,我们可以用泊松分布对其建模,并估计参数 μ 。我们将用频率论方法和贝叶斯方法两种方法来估计。
123456 | 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
从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。
由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。
代码:
1234567891011121314 | fig = plt.figure(figsize=(11,3))ax = fig.add_subplot(111)x_lim = 60mu = [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
是计数数据,我们可以用泊松分布对其建模,并估计参数 μ 。我们将用频率论方法和贝叶斯方法两种方法来估计。
123456 | 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值直到找到可能性最大的值。
1234567 | 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\’]) |
12on-num\” data-line=\”crayon-5812b1408f85c085651845-3\”>3 | The estimated value of mu is: 18.0413533867CPU times: user 63 µs, sys: 1e+03 ns, total: 64 µsWall time: 67.9 µs |
所以,μ 的估计值是18.0413533867。优化技术没有对不确定度进行评估,它只返回一个点,效率很高。
下图描述的是我们优化的函数。对于每个μ值,图线显示给定数据和模型在μ处的似然度。优化器以登山模式工作——从曲线上随机一点开始,不停向上攀登直到达到最高点。