用 Python 进行贝叶斯模型建模(1)
admin
2023-07-31 00:46:24
0

本系列:

  • 《第 0 节:导论》
  • 《第 1 节:估计模型参数》
  • 《第 2 节:模型检验》
  • 《第 3 节:分层模型》
  • 《第 4 节:贝叶斯回归》
  • 待续

第1节:估计模型参数

在这一节,我们将讨论贝叶斯方法是如何思考数据的,我们怎样通过 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

从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。

由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。

 p(x \\ | \\ \\mu) = \\frac{e^{-\\mu}\\mu^{x}} {x!} \\mbox{ for }x = 0, 1, 2, \\cdots

\\lambda = E(x) = Var(\\mu)

代码:

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

从贝叶斯方法的角度看,这个数据是由随机过程产生的。但是,既然数据被观测,它便固定了并且不会改变。这个随机过程有些模型参数被固定了。然而,贝叶斯方法用概率分布来表示这些模型参数的不确定性。

由于这个男孩调查的是计数(非负整数),一个通常的做法是用泊松分布对数据(如随机过程)建模。泊松分布只有一个参数 μ,它既是数据的平均数,也是方差。下面是三个不同 μ 值的泊松分布。

 p(x \\ | \\ \\mu) = \\frac{e^{-\\mu}\\mu^{x}} {x!} \\mbox{ for }x = 0, 1, 2, \\cdots

\\lambda = E(x) = Var(\\mu)

代码:

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

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

相关内容

热门资讯

Mobi、epub格式电子书如... 在wps里全局设置里有一个文件关联,打开,勾选电子书文件选项就可以了。
500 行 Python 代码... 语法分析器描述了一个句子的语法结构,用来帮助其他的应用进行推理。自然语言引入了很多意外的歧义,以我们...
定时清理删除C:\Progra... C:\Program Files (x86)下面很多scoped_dir开头的文件夹 写个批处理 定...
scoped_dir32_70... 一台虚拟机C盘总是莫名奇妙的空间用完,导致很多软件没法再运行。经过仔细检查发现是C:\Program...
65536是2的几次方 计算2... 65536是2的16次方:65536=2⁶ 65536是256的2次方:65536=256 6553...
小程序支付时提示:appid和... [Q]小程序支付时提示:appid和mch_id不匹配 [A]小程序和微信支付没有进行关联,访问“小...
pycparser 是一个用... `pycparser` 是一个用 Python 编写的 C 语言解析器。它可以用来解析 C 代码并构...
微信小程序使用slider实现... 众所周知哈,微信小程序里面的音频播放是没有进度条的,但最近有个项目呢,客户要求音频要有进度条控制,所...
Apache Doris 2.... 亲爱的社区小伙伴们,我们很高兴地向大家宣布,Apache Doris 2.0.0 版本已于...
python清除字符串里非数字... 本文实例讲述了python清除字符串里非数字字符的方法。分享给大家供大家参考。具体如下: impor...