在这里,我们将帮助客户将 PyMC3 用于两个贝叶斯推理案例研究:抛硬币和保险索赔发生。
回想一下,我们最初的贝叶斯推理方法是:
- 设置先前的假设,并根据启发式、历史或样本数据建立我们数据的“已知已知”。
- 形式化问题空间和先前假设的数学模型。
- 正式化先前的分布。
- 应用贝叶定理从观察到的样本数据中推导出后验参数值。
- 重复步骤 1-4,以获取更多数据样本。
可下载资源
使用 PyMC3,我们现在可以简化和压缩这些步骤。
首先,我们设定先验信念和先验β-二项分布。
prior_beta = prior_beta.pdf(theta) / prior_beta.pdf(theta).sum() # 样本积分 [pmf]()
plt.legend();
其次,我们定义并检查我们的样本观察数据
视频
R语言bnlearn包:贝叶斯网络的构造及参数学习的原理和实例
视频
贝叶斯推断线性回归与R语言预测工人工资数据
视频
R语言中RStan贝叶斯层次模型分析示例
print(f'Observed P(tails) = {tails/trials}')
第三,我们定义并运行我们的数学模型
请注意,PyMC3 提供了一种干净有效的语法来描述先验分布和观测数据,我们可以从中包括或单独启动模型抽样。另请注意,PyMC3 允许我们定义先验、引入样本观察数据并启动后验模拟。
# [NUTS](),采样器(汉密尔顿式)
step = pm.NUTS()
结果
或者通过更多的采样和更多的链。然后,跟踪摘要返回有用的模型性能摘要统计信息:
- mc_error通过将迹线分解为批次,计算每个批次的平均值,然后计算这些平均值的标准偏差来估计模拟误差。
- hpd_* 给出最高的后密度区间。2.5 和 97.5 标签有点误导。有很多 95% 的可信区间,具体取决于左右尾巴的相对权重。95% HPD 区间是这 95% 区间中最窄的。
- Rhat有时被称为潜在的规模缩减因子,它为我们提供了一个因子,如果我们的MCMC链更长,则可以减少方差。它是根据链与每个链内的方差来计算的。接近 1 的值很好。
summary
我们使用迹线手动绘制和比较先验分布和后验分布。确认这些与手动获得的相似,后验分布均值为 P(Tails|观测数据)= 0.35。
但是,PyMC3还提供了创建迹线图,后验分布图。
pm.traceplot(trace)
pm.plot_posterior(trace,ref_val=0.5);
我们有它。PyMC3 和其他类似软件包提供了一组简单的函数来组装和运行概率模拟,例如贝叶斯推理。
个案研究:
使用贝叶斯推理评估保险索赔发生率
保险索赔通常被建模为由于泊松分布式过程而发生。
泊松分布由下式给出:
其中 lambda λ 是事件的“速率”,由事件总数 (k) 除以数据中的单位数 (n) 给出 (λ = k/n)。在泊松分布中,泊松分布的期望值 E(Y)、均值 E(X) 和方差 Var(Y) 相同;
随时关注您喜欢的主题
例如,E(Y) = E(X) = Var(X) = λ。
请注意,如果方差大于均值,则称数据过于分散。这在具有大量零的保险索赔数据中很常见,并且最好由负二项式和零膨胀模型(如 ZIP 和 ZINB)处理。
一、建立先验分布
在这里,我们生成一些观测数据,这些数据遵循泊松分布,速率为 lambda,λ = 2。
n = 1000
lam_ = 2
axs.set_title('Histogram: Simulated Poisson $y$')
axs.set_xlabel('Poisson lambda=λ')
axs.set_ylabel('P(λ)')
axs.legend();
我们可以使用β泊松,或任何类似于观察到的λ数据形状的分布,但是伽马泊松最适合:
泊松可以取任何正数到无穷大(0,∞),而β或均匀是[0-100]。- 伽马和泊松属于同一分布家族。
- 伽马的峰值接近于零。
- 伽马尾巴走向无穷大。
伽马泊松先验为:
注意在 scipy 中,伽马分布使用形状 a 和尺度参数化,其中速率 b 等于尺度的倒数(速率 = 1/尺度)。
其中 a 是伽马形状,b 是伽马速率参数。伽马密度函数为:
其中 a>0 是形状参数,b>0 是速率参数,以及
和
prior = lambda x: stats.gamma.pdf(x, a=a, scale=rate,loc=0)
priors = prior(x)
# 画图
axs.plot(x, priors, 'r-',label='Gamma')
二、似然函数与后验
伽马函数通常被称为广义阶乘,因为:
sp.gamma(n+1) == math.factorial(n)
True
则似然函数为:
然后作为
后向分布再次为伽马
def posterior(lam,y):
shape = a + y.sum()
如图所示,后验平均值(蓝色)以我们在开始时设置的真实 lambda 速率为中心。后验平均值为:
即后验平均值是先验平均值和观测样本平均值的加权平均值
posterior mean: {(a+y.sum()) / (b+y.size)}
sample mean:{y.mean()}""")
现在让我们在 PyMC3 中重现上述步骤。
print(a,b,lam_,y.shape)
with model:
# 定义参数 lambda 的先验值。
prior_lam = pm.Gamma('prior-gamma-lambda', alpha=a, beta=b)
迹线图显示每个模拟的结果。
低于平均值、分位数、可信区间 (HPD) 94% 和任意参考值(橙色垂直)。
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore")
您可能已经注意到,在这个例子中,我们已经根据观察到的数据定义了我们的先验分布,并对该数据应用贝叶斯推理来推导出后验分布,确认 lambda 为 2。
结论:
在这篇文章中,PyMC3 被应用于对两个示例进行贝叶斯推理:使用 β-二项分布的抛硬币偏差,以及使用 gamma-泊松分布的保险索赔发生。
可下载资源
关于作者
Kaizong Ye是拓端研究室(TRL)的研究员。在此对他对本文所作的贡献表示诚挚感谢,他在上海财经大学完成了统计学专业的硕士学位,专注人工智能领域。擅长Python.Matlab仿真、视觉处理、神经网络、数据分析。
本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。
非常感谢您阅读本文,如需帮助请联系我们!