全文链接：https://tecdat.cn/?p=33416
在这里，我们将帮助客户将 PyMC3 用于两个贝叶斯推理案例研究：抛硬币和保险索赔发生（点击文末“阅读原文”获取完整代码数据）。
方法：
回想一下，我们最初的贝叶斯推理方法是：
- 设置先前的假设，并根据启发式、历史或样本数据建立我们数据的“已知已知”。
- 形式化问题空间和先前假设的数学模型。
- 正式化先前的分布。
- 应用贝叶定理从观察到的样本数据中推导出后验参数值。
- 重复步骤 1-4，以获取更多数据样本。
使用 PyMC3，我们现在可以简化和压缩这些步骤。
首先，我们设定先验信念和先验β-二项分布。
prior_beta = prior_beta.pdf(theta) / prior_beta.pdf(theta).sum() # sample integral pmf plt.legend();
其次，我们定义并检查我们的样本观察数据
print(f'Observed P(tails) = {tails/trials}')
第三，我们定义并运行我们的数学模型
请注意，PyMC3 提供了一种干净有效的语法来描述先验分布和观测数据，我们可以从中包括或单独启动模型抽样。另请注意，PyMC3 允许我们定义先验、引入样本观察数据并启动后验模拟。
# NUTS, the No U-Turn Sampler (Hamiltonian) step = pm.NUTS()
结果
或者通过更多的采样和更多的链。然后，跟踪摘要返回有用的模型性能摘要统计信息：
- mc_error通过将迹线分解为批次，计算每个批次的平均值，然后计算这些平均值的标准偏差来估计模拟误差。
- hpd_* 给出最高的后密度区间。2.5 和 97.5 标签有点误导。有很多 95% 的可信区间，具体取决于左右尾巴的相对权重。95% HPD 区间是这 95% 区间中最窄的。
- Rhat有时被称为潜在的规模缩减因子，它为我们提供了一个因子，如果我们的MCMC链更长，则可以减少方差。它是根据链与每个链内的方差来计算的。接近 1 的值很好。
summary
我们使用迹线手动绘制和比较先验分布和后验分布。确认这些与手动获得的相似，后验分布均值为 P（Tails|观测数据）= 0.35。
点击标题查阅往期内容
贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析免疫球蛋白、前列腺癌数据
01
02
03
04
但是，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]。
- 伽马和泊松属于同一分布家族。
- 伽马的峰值接近于零。
- 伽马尾巴走向无穷大。
伽马泊松先验为：
其中 a 是伽马形状，b 是伽马速率参数。伽马密度函数为：
其中 a>0 是形状参数，b>0 是速率参数，以及
和
注意在 scipy 中，伽马分布使用形状 a 和尺度参数化，其中速率 b 等于尺度的倒数（速率 = 1/尺度）。
prior = lambda x: stats.gamma.pdf(x, a=a, scale=rate,loc=0) priors = prior(x) # plot 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: # Define the prior of the parameter 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-泊松分布的保险索赔发生。
