一种基于无信息先验下泊松分布变点的贝叶斯估计

许 婷，吴有富，张英雪

（1.贵州民族大学数据科学与信息工程学院，贵州 贵阳 550025；2.贵州交通职业技术学院 贵州 贵阳 550025）

变点问题的贝叶斯估计注重的是先验分布的选取，不同先验分布的选取会对估计精度有不同的影响。对于泊松分布参数变点的研究，张梦琇[1]选取先验分布为共轭先验，并通过二分法实现了对泊松分布多变点模型的参数估计。范元静[2]在选择先验分布为共轭先验伽马分布的基础上，运用RJMCMC算法对泊松分布的变点参数进行数量及位置的估计。何朝兵[3]选取共轭先验分布作为泊松分布的先验分布，并利用贝叶斯方法对IIRCT下单变点泊松分布的参数进行估计。钟颖和田茂[4]再选取泊松分布的先验分布为共轭先验分布，并使用贝叶斯方法对单变点泊松分布进行参数估计。胡兴[5]使用Fisher信息阵确定了泊松分布的先验分布，并通过MCMC方法对推导的后验分布进行模拟研究。对于泊松分布变点问题的研究，人们几乎都是选择泊松分布的共轭先验分布伽马分布作为其先验分布对其变点问题进行研究，在确定无信息先验分布的方法下对于泊松分布先验分布的选取也仅仅是利用Fisher信息阵确定了其先验分布，而对于其他无信息先验下的泊松分布变点问题的估计目前尚未报道。文章基于此，在Zhang等[6]讨论的确定无信息先验分布方法的基础上，基于“变换不变性”原则推导出Poisson分布参数的无信息先验分布，在参数的全条件后验分布下对Poisson分布尺度参数变点模型的参数估计问题进行研究。

1 Poisson分布单变点模型

2 极大似然估计

3 贝叶斯估计

贝叶斯估计作为常用的对参数进行统计推断的方法，在使用贝叶斯方法对变点问题进行研究时，是将模型中把变点涵盖在内的所有参数都视为随机变量。通过引入先验分布，先验分布是反映在抽样前对参数的认识，对参数设定不同的先验分布，根据设定的先验分布和已知样本服从的分布计算出每个参数的后验分布，后验分布是反映在抽样后对参数的认识，根据每个参数全条件下的后验分布对参数进行估计。由于不同的先验分布会对参数估计结果有不同的影响，所以在选取先验分布的时候很重要。文章选用的是无信息先验分布，在此基础上对泊松分布的变点进行贝叶斯估计。

3.1 不变性原则下泊松分布的无信息先验分布

3.2 无信息先验下泊松分布单变点参数的贝叶斯估计

首先确各参数的先验分布：

4 数值仿真

考虑下面的变点模型：

（1）极大似然估计法

根据推导出的Poisson分布的对数似然函数式（4）和各参数的极大似然估计式（8）对参数，k进行估计，使用R软件进行100次模拟的结果如表1所示：

表1 参数k12的极大似然估计

表1 参数k12的极大似然估计

参数 真值 估计值 相对误差1 10 11.26 0.126 20 21.51 0.076 k 40 42.73 0.068 2

（2）贝叶斯估计

表2 参数k12的贝叶斯估计

表2 参数k12的贝叶斯估计

参数 真值Gibbs样本1 10 10.01 0.0010.000139.672 9.99510.346 20000 2 20 19.410.0295 0.0002919.03 19.41 19.79 20000 k 40 39.980.0005 0.00057 40 39.98 40.00 20000估计值相对误差MC误差2.5%分位数中位数97.5%分位数

图1 1迭代轨迹图

图2 2迭代轨迹图

图3 k迭代轨迹图

图4 1后验分布的核密度估计图

图5 2后验分布的核密度估计图

图6 k后验分布的核密度估计图

虽然从图1到图3中可以看出各参数产生的马尔科夫链已经是平稳状态了，但是为了判定生成的Gibbs样本是不是收敛的，在用MCMC方法进行模拟的时候，会同时产生多个迭代链，若生成的多个迭代链逐渐稳定并且趋于重合，那么可以说明Gibbs抽样是收敛的。由于文章主要是是为了判定变点的位置，因此在模拟的过程中，固定了变点k的位置，用不同的的值分别进行20000次迭代，其两组初始值分别为，产生的迭代链轨迹图如图 7—图 9所示。其中图 7是下的迭代链，从图7中可以看出，其抽样基本上都在变点附近波动，且波动的幅度较小，说明估计的效果较好；图8是下的迭代链，同样可以看出，其抽样基本在变点附近波动，抽样的波动程度同样较小，估计的效果也较好；图9是两组初始值下同时产生的迭代链，从图中可以看出，产生的两个马尔科夫链都是稳定且趋于重合的，表明Gibbs样本是收敛的。

图7 k第一组初始值下的迭代链

图8 k第二组初始值下的迭代链

图9 k两组初始值下的迭代链

表3 不同样本量下参数k12的贝叶斯估计

表3 不同样本量下参数k12的贝叶斯估计

总样本量n 100 200 500 1000 2000 5000 1样本量 40 110 300 400 1100 3000 2样本量 60 90 200 600 900 2000 1估计值 10.01 10.30 10.16 9.94 10.05 2估计值 19.41 20.25 19.94 20.09 19.81 19.85 k估计值 39.98 111.35300.07 399.861099.64 2999.05

文章还给出了总样本量 n=100，200，500，1000，2000，5000时的三个参数的估计值，如表2所示，其中样本量的值就是变点的位置。从表中可以看出，在给定的的情况下，不管n是小样本还是大样本，每个参数的估计值都与真实值非常接近。且不管样本量增大到多少，其估计精度都还是挺高的。因此从模拟的结果来看，参数的估计值与真实值都是很接近的，表明MCMC的模拟效果是比较好的。

5 实证分析

我们将我们的方法应用于Jarrett[8]给出的英国煤灾事故数据集，他修正并扩展了Maguire[9]等人给出的数据集。

表4 英国年度煤矿灾害数据，1851-1962年

Carlin等人[10]用共轭先验得出的矿难发生变化的点是在41，也就是在1890年底至1892年初之间矿难次数开始减少，用我们的方法得到的结果与其是差不多的，Raftery和Akman[11]用模糊先验也给出了相似的结果。

从图10中可以看出，发生突变的位置是40附近，也就是位置41的后验概率最大，其次是位置40和位置39，对应下来也就是1890年附近，也即是意味着这种变化最有可能发生在1889年末至1892年初之间。且从图11中可以看到，发生变化之前平均每年要发生3次左右的矿难，而在变化点之后每年发生的矿难次数平均是1次左右。

图10 变点位置后验分布密度曲线

图11 变点前年均矿难次数密度曲线

图12 变点后年均矿难次数密度曲线

用MCMC抽样方法对矿难灾害的估计结果如下：

表5 参数估计

从表5中可以看出，MC误差较小，说明MCMC方法的收敛性较好。将用文章的方法得出的估计值与Carlin等人研究的结果相比，相对误差都小于0.1，说明估计效果较好。从估计结果来看，英国煤矿灾害变点位置之后发生灾难的次数很大程度降低的原因可能是与19世纪英国颁布的一系列有关采矿安全等法案有关，实施之后对采矿安全有了一定程度上的影响。

6 结论

文章主要在已有的推导无信息先验分布的方法下，基于不变性理论下推导出Poisson分布的无信息先验分布，利用贝叶斯方法实现对参数的估计，并和极大似然方法做比较。

根据表1和表2中的估计值可以看出，两种方法都能估计出变点位置参数 k的位置和尺度参数的值。但是通过比较两种方法的相对误差可以知道极大似然估计的精确度没有贝叶斯估计的高，极大似然方法得出的估计值与给定的真实值相差较大，而贝叶斯方法得出的估计值与给定的真实值相差较小。因此可以知道用贝叶斯方法估计变点位置参数k的效果更好一些。根据表3中不同样本量下的估计值，k可以知道，总样本量n增加时，各个参数的估计值依然与真实值很接近，说明贝叶斯方法在小样本和大样本下都可以用来估计参数值，且估计效果比极大似然方法好。

并用文章使用的先验应用于英国煤灾事故数据集，通过实证分析的结果可以看出，与前人的研究结果是基本一致的。综上所述，用贝叶斯方法对变点位置参数k和尺度参数的估计效果都是非常不错的，因此对于推导出来的Poisson分布的尺度参数的无信息先验为是可行且精度较高的，使用此先验分布对Poisson分布的变点问题进行研究是有效的。