基于伽马分布的变点检测.pdf

2023 年第 8 期73计算机应用信息技术与信息化基于伽马分布的变点检测王轩云1 黄介武1WANG Xuanyun HUANG Jiewu 摘要 针对伽马分布中变点个数及变点位置的检测问题，本文基于三类不同的惩罚项，采用不同方法对伽马分布中单变点及多变点进行检测，对变点个数及位置识别的准确性进行研究，从而比较不同方法的检测效果。数值模拟表明，针对单变点的情况，AMOC 方法的检测效果优于其他方法；而在多变点的情况下，所采用的 WBS 方法具有更好的检测性能，优于其他方法。最后，本文通过实证分析验证不同方法的有效性。关键词 变点检测；伽马分布；惩罚项；单变点；多变点doi：10.3969/j.issn.1672-9528.2023.08.0161.贵州民族大学数据科学与信息工程学院 贵州贵阳 550025 基金项目 贵州省科技计划基金项目（黔科合基础20171083 号）；贵州省基础研究计划（软科学）（黔科合支201920001）0 引言变点检测是一个普遍存在的问题，变点检测是确定时间序列的特征是否发生了变化的主要方法。更正式的说，假设存在一个有序的数据序列：y1,y2,yn，在某一时刻发生了变化使得发生变化之前和之后的序列的统计性质不同，则该点就是变点。关于变点检测方法，众多学者对其进行了研究，例如，Steinebach1利用似然比检验统计量，研究了一般线性回归模型中只有一个变点(AMOC)的检测问题；Xu 等人2基于剪枝精确线性时间(PELT)，提出了一种有效的在线优化和变化点检测算法；李等人3利用二元分割方法(BS)，对多元 Logistic 回归模型中存在的变点进行了研究；Gordon等人4利用野生二元分割法(WBS)，对非平稳时间序列中变点的数量及位置进行了研究。伽马分布在很多的领域都有着应用，比如产品的可靠性评估、生存分析、地震的预测等。例如，王等人5在层次贝叶斯框架下提出了的两阶段伽马过程模型对单调非光滑的退化数据建模，分析了一个抗辐射性能的可靠性案例。此外，Wang 等人6对存在变化点的降解情况，采用多相降解模型检测降解过程，找到了针对被监测产品的实时可靠性的评估方法。伽马分布在实际应用领域中具有重要意义，但其变点检测研究却较少。针对如何检测服从伽马分布的数据的变点问题，通过模拟实验来检验不同方法的检测效果，从而在实际应用中提供有价值的实施策略。本研究填补了伽马分布变点检测方面的研究空白，对于相关领域的进一步研究和实际应用具有重要意义。在实际应用中，如何选取最佳的变点检测方法仍是一个具有挑战性的问题。目前，针对检测变化点信号的方法，缺乏直接的比较来评估其相对性能，这使得方法之间的差异难以评价。因此，本文考虑了 atmost one change point（AMOC）7、pruned exact linear time（PELT）8、binary segmentation（BS）9、Wild binary segmentation（WBS）10这四种变点检测方法，通过探讨这四种方法识别的变点个数、变点位置的正确率，得出不同检测方法在伽马分布中的优良性，从而可以更好的选取检测方法。1 模型介绍1.1 伽马分布的变点模型伽马分布的概率密度函数为1(,)()yf yyeI=,记为y G(,)，其中，I 为示性函数，10()exp()(0)x xdx=。其中均值为()/E y=，方差为2()/Var y=。伽马分布的单变点模型为:111221(,),1,.,(1)(,),1,.,iGaibyGaibn=+（1）其中：需要估计的是：变点的个数和变点的位置 b1。伽马分布的两个变点模型为：（2）1 112212332(),1,.,(),1,.,(2)(),1,.,iGaibyGaibbGaibn =+=+2023 年第 8 期74计算机应用信息技术与信息化其中：需要估计的是：变点的个数，变点的位置 b1，b2。2 检测方法2.1 AMOC 方法只有一个变点（atmost one change point，AMOC）方法6是一种检测单个变点的方法，该方法通过检验适当的似然函数来检测变点。详细过程如下：假设存在两种关系111 (3)ya xb=+（3）222 (4)ya xb=+（4）式中：1和 2是独立正态分布的误差项，均值为零，标准差为 1、2。这两个真实的关系产生了一个观测值，需要估计的问题是系统从一个状态切换到另一个状态的点。假设第一个 t 观测值由(3)产生，最后一个 T-t 观测值由(4)产生。为了估计 t 和(3)，(4)中的参数，继续如下操作，第 i 点的 1和第 j 点的 2的密度为：2211111/2)exp(1/2)()iiya xb（和 2222221/2)exp(1/2)()jjya xb（因此，来自(3)的 t 观测样本和来自(4)的 T-t 观测样本的似然值为：211211111()exp()22ttjjiya xb=和222212211()exp()22TT tjjj tya xb=+则整个样本的似然值：21121111222212111()()exp()2221()2ttT tiiiTjjj tya xbya xb=+取似然函数的对数：12221122221112log2log()log11()()22tTiijjij tLTtTtya xbya xb=+=（5）将上式(5)中 1，b1，2，b2的偏导数设为零，求解得到最小二乘估计：211112211()tttiiiiiiittiiiitx yxytxx=；1111ttiiiiyxbatt=；同样将式(5)中的 1和 2的偏导数设为零，并代入 1，b1，2，b2的估计值，得到：211211)tiiiya xbt=（，222122)Tjjj tya xbTt=+=（把2212 ,代入(5)得到的最大似然估计：12()log 2log()log/2L tTtTtT=（6）它给出了给定值 T 的最大似然的对数，是 T 单独的函数。现在要找到(6)式中的 t 的值，计算(6)式中所有 t 可能值的似然函数的值，并选择对应于最大 t 值的 t 值作为最大似然估计值。AMOC 的执行程序如下：算法 1：AMOC 算法步骤 1：根据时间段11(,),.,(,)TTx yxy对观测结果进行排序，并将数据分为左手组和右手组，分别估计两组的回归线。步骤 2：然后将两组间的分界点向右移动一个时间单位，向左移动一个单位。计算每个新分界点下新的左右两组的回归线再次移动分界线，并以类似的方式迭代进行处理。2.2 PELT 方法(1)剪枝精确线性时间（pruned exact linear time,PELT）方法7是以最佳分区算法（optimal partitioning,OP）为基础，通过递归的方式来优化目标函数来检测变点。在处理的过程中通过加入一个剪枝步骤来降低计算成本。PELT 的执行程序如下:算法 2：PELT 算法Input:一组数据12,.,ny yy适合度测量 C()依赖于数据一个不依赖于变化点的数量或位置的惩罚常数 一个满足方程(1):(1):(1):()()()tsttTC yC yKC y+的常数 KInitialize:设 n 等于数据长度并设置1(0),(0),0FcpNULL R=Iterate:for *=1,n1.计算(1):()min()()RFFC y+=+2.设(1):()()argminRFC y+=3.设置()(),cpcp=4.设置1(1):()()()RRFC yKF+=+Output:在 cp(n)中记录的变化点其中iyR。(2)最佳分区算法是将成本函数的最优值与最后一个变化点之前数据的最优划分的成本加上从最后一个变化点到数据结束的段的成本联系起来。OP 执行程序如下：算法 3：OP 算法Input:一组数据12,.,ny yy适合度测量 C()依赖于数据一个不依赖于变化点的数量或位置的惩罚常数 一个满足方程(1):(1):(1):()()()tsttTC yC yKC y+的常数 KInitialize:设 n 等于数据长度并设置(0),(0)FcpNULL=Iterate:for*=1,n1.计算0(1):()min()()FFC y +=+2.设0(1):argmin()()FC y +=+3.设置()(),cpcp=Output:在 cp(n)中记录的变化点其中iyR。2023 年第 8 期75计算机应用信息技术与信息化2.3 BS 算法BS 方法执行程序如下：算法 4：BS 算法function (,)TBINSEG s eif 1es then将 b0添加到估计的变点集合中0(,)TBINSEG s b0(1,)TBINSEG be+else stopend ifend ifend function其中，T是一个阈值参数。二元分割(binary segmentation，BS)方法8的思想是：首先，在整个时间序列中搜索一个变点，如果检测到一个变点，就将数据在该点分为两段，然后在变点之前和之后的两个数据集上重复单一的更改点过程。如果在这两个数据集上检测到变点，则将它们分割，重复此过程，直到所有的变点都检测完。2.4 WBS 方法野生二元分割（wild binary segmentation，WBS）方法9是基于二元分割的思想，继承 BS 的主要优点改进二元分割的缺点，即 BS 对两个变点位置间隔较近的区间进行估计时，会出现少估计的问题。野生二元分割的思想是在第一阶段，不是使用整个数据样本1,.,TXX的全局 CUSUM 统计量，而是随机抽取一些子样本，即向量1,.,sseXXX+，其中 s 和 e为整数，1seT，并计算每个子样本的 CUSUM 统计量。然后，最大化每个 CUSUM，在整个 CUSUMs 集合中选择最大的最大化值，并将其作为特定阈值进行测试的第一个变点候选点。检测到第一个变点，则它的左右重复同样的过程。WBS 执行程序如下：算法 5：WBS 算法function(,)TWILDBINSEG s eif 1es then将 b0添加到估计的变点集合中0(,)TWILDBINSEG s b0(1,)TWILDBINSEG be+elsestopend ifend ifend function其中，T是一个阈值参数。3 随机模拟及实证分析3.1 随机模拟在伽马分布的基础上，分别讨论两种情况：（1）只有一个变点；利用 AMOC、PELT、BS、WBS 这四种方法进行检测。（2）有两个变点。由于 AMOC 方法是检测单变点的方法，所以对两个变点的情况这里利用 PELT、BS、WBS这三种方法进行检测。样本量选取为(1500,2000)N=。单变点的参数选取分别为 2 0.5 和 5 2；两个变点的参数的选取为 2 0.5 2 和 5 2 5。在本文中，对于这四种方法采用 AIC11、SIC12、mBIC13这三种信息准则作为惩罚项来进行检测。每种参数变化的情况生成 5 组数据来进行重复实验。通过模拟实验，对伽马分布中变点个数及变点位置检测的情况，比较 AMOC、PELT、BS、WBS这四种方法对伽马分布的检测效果，并给出方法的选择建议。3.1.1 伽马分布的单变点检测效果模拟伽马分布中样本量为1000的单变点图，如图1所示。图 1 模拟仿真数据的伽马分布单变点数据图由表 1 可以看出，不同变点检测方法在伽马分布变点检测问题中的性能表现存在一定差异。在存在一个变点的情况下，AMOC 方法能够准确地检测到所有的变点，显示出最佳的检测正确率；而 WBS 方法在此情况下的检测表现略逊于其他三种方法。2023 年第 8 期76计算机应用信息技术与信息化表 1 变点个数为一个时，不同方法在伽马分布中正确检测变点位置的个数方法样本量参数变化惩罚方式 正确的个数正确率AMOCN=1500 :2 0.5AICSICmBIC1,1,1,1,1 1,1,1,1,1 1,1,1,1,1100%100%100%:5 2AICSICmBIC1,1,1,1,1 1,1,1,1,1 1,1,1,1,1100%100%100%N=2000:2 0.5AICSICmBIC1,1,1,1,1 1,1,1,1,1 1,1,1,1,1100%100%100%:5 2AICSICmBIC1,1,1,1,1 1,1,1,1,1 1,1,1,1,1100%100%100%PELTN=1500 :2 0.5AICSICmBIC1,1,1,0,1 1,1,1,1,1 1,1,1,1,190%100%100%:5 2AICSICmBIC1,1,1,0,1 1,1,1,1,1 1,1,1,1,190%100%100%N=2000:2 0.5AICSICmBIC1,1,0,1,1 1,1,0,1,1 1,1,1,1,190%90%100%:5 2AICSICmBIC1,1,1,1,1 1,1,1,1,1 1,1,1,1,1100%100%100%Binary Seg-mentationN=1500 :2 0.5AICSICmBIC1,1,1,0,0 1,1,1,1,1 1,1,1,1,180%100%100%:5 2AICSICmBIC1,1,1,0,0 1,1,1,1,1 1,1,1,1,180%100%100%N=2000:2 0.5AICSICmBIC1,1,1,0,0 1,1,1,0,0 1,1,1,1,180%80%100%:5 2AICSICmBIC1,1,1,0,0 1,1,1,1,1 1,1,1,1,180%100%100%Wild Bina-ry Segmen-tationN=1500 :2 0.5AICSICmBIC1,1,0,0,11,1,1,0,11,1,1,1,180%90%100%:5 2AICSICmBIC1,1,1,0,01,1,1,0,11,1,1,0,180%90%90%N=2000:2 0.5AICSICmBIC0,1,0,1,10,1,1,1,11,1,1,1,180%90%100%:5 2AICSICmBIC1,0,1,1,11,0,1,1,11,1,1,1,190%90%100%由表 2 不同方法检测的伽马分布的变点总数可以看出，存在一个变化点时，AMOC 方法在三种惩罚下均能准确地检测到变点的总数，表现最为优异。此外，PELT 和 BS 方法在mBIC 惩罚下的表现也相对较好。相比之下，WBS 方法的检测效果不如其他方法，存在过估计的问题，特别是在 AIC 惩罚下，产生过多的过估计个数。表 2 变点个数为一个时，不同方法在检测伽马分布中的变点总数模拟条件 方法比较及标准惩罚方式AMOCPELTBinarySegmenta-tionWild BinarySegmentationN=1500,n=1:2 0.5AICSICmBIC1,1,1,1,11,1,1,1,11,1,1,1,13,1,4,4,11,1,2,3,11,1,1,1,15,5,5,3,53,2,2,1,11,1,1,1,150,49,50,49,493,1,4,4,13,2,2,3,1N=2000,n=1:2 0.5AICSICmBIC1,1,1,1,11,1,1,1,11,1,1,1,12,1,3,1,21,1,3,1,31,1,1,1,15,5,5,5,52,1,1,1,21,1,1,1,150,50,50,50,493,4,4,4,13,3,3,1,1N=1500,n=1:2 0.5AICSICmBIC1,1,1,1,11,1,1,1,11,1,1,1,13,2,2,2,41,1,1,2,11,1,1,1,15,5,5,5,32,1,1,1,11,1,1,1,150,50,49,49,83,1,4,4,13,1,2,2,1N=2000,n=1:2 0.5AICSICmBIC1,1,1,1,11,1,1,1,11,1,1,1,15,4,3,3,23,2,2,3,11,1,1,1,15,5,5,5,21,2,1,1,11,1,1,1,150,50,50,50,74,1,4,4,12,1,2,3,13.1.2 伽马分布的两个变点检测效果图 2 是模拟伽马分布中样本量为 1000 的两个变点图。图 2 模拟仿真数据的伽马分布两个变点数据图由表 3 可以看出，存在两个变点时，WBS 方法相比于其他两种方法具有稍优的变点识别正确率。研究发现，随着样本量的增大和参数的增加，变点识别的正确率明显提高，说明在伽马分布中，参数的选取对正确检测变点个数也有影响。此外，从惩罚类型的角度来看，在 mBIC 惩罚下正确识别变点的个数的正确率要优于 AIC 和 SIC 惩罚。表 3 变点个数为两个时，不同方法在伽马分布中正确检测变点位置的个数方法样本量参数变化惩罚方式 正确的个数正确率PELTN=1500 :2 0.5 2AICSICmBIC1,2,1,2,11,2,1,2,11,2,1,2,270%70%80%:5 2 5AICSICmBIC2,2,1,2,11,2,2,2,22,2,1,2,280%90%90%N=2000:2 0.5 2AICSICmBIC1,2,2,1,22,2,2,1,22,2,2,1,280%90%90%:5 2 5AICSICmBIC2,2,2,2,12,2,2,2,12,2,2,2,290%90%100%2023 年第 8 期77计算机应用信息技术与信息化表 3（续）方法样本量参数变化惩罚方式 正确的个数正确率Binary Segmen-tationN=1500 :2 0.5 2AICSICmBIC1,2,1,2,12,2,1,2,12,2,1,2,270%80%90%:5 2 5AICSICmBIC1,2,1,2,21,2,2,2,22,2,2,2,280%90%100%N=2000:2 0.5 2AICSICmBIC2,1,2,1,22,2,1,2,22,2,1,2,280%90%90%:5 2 5AICSICmBIC1,1,2,2,22,2,2,2,22,2,2,2,280%100%100%Wild Binary Segmen-tationN=1500:2 0.5 2AICSICmBIC2,2,1,2,12,2,1,2,12,2,2,2,280%80%100%:5 2 5AICSICmBIC1,1,2,2,21,2,2,2,22,2,2,2,280%90%100%N=2000:2 0.5 2AICSICmBIC2,2,2,1,12,2,2,2,12,2,2,2,280%90%100%:5 2 5AICSICmBIC2,2,1,1,22,2,2,1,22,2,2,2,280%90%100%表 4 变点个数为两个时，不同方法在伽马分布中检测的变点总数模拟条件 方法比较及标准惩罚方式AMOCPELTBinarySegmenta-tionWild BinarySegmentationN=1500,n=2 :20.52AICSICmBIC4,3,3,5,53,2,2,4,32,2,3,2,25,5,5,5,22,2,4,4,12,2,2,2,250,50,50,50,503,4,3,3,22,2,2,2,2N=2000,n=2 :20.52 AICSICmBIC3,3,3,4,32,2,2,3,32,2,2,3,25,5,5,5,52,2,2,3,22,2,2,3,250,49,47,50,503,4,3,3,22,2,2,3,2N=1500,n=2 :20.52AICSICmBIC3,3,3,3,32,2,3,4,32,2,2,2,25,5,5,5,53,4,3,3,22,2,2,2,250,50,50,50,493,4,3,3,22,2,2,2,2N=2000,n=2 :20.52AICSICmBIC3,7,4,4,42,4,3,2,42,3,2,2,25,5,5,5,43,3,3,4,32,2,3,2,250,50,50,50,493,4,3,3,22,2,2,2,2从表 4 中可以看出，在存在两个变点的情况下，不同变点检测方法的表现也存在差异。WBS 方法在检测总变点数方面表现最为优异。而 PELT、BS、WBS 方法在 AIC 和SIC 惩罚下存在过估计问题，特别是在 AIC 惩罚下，WBS方法的过估计问题比较严重。相比之下，在 mBIC 惩罚下，三种方法均有较好的检测效果，并且能够最准确地检测到总变点数。为了消除 SIC 产生的虚假点，把增加惩罚到，而不是。此更改是通过将惩罚类型更改为来实现。表 5 增加惩罚的 SIC 对伽马分布的检测效果变点个数参数变化样本量正确的个数 正确率 检测到的个数n=1:2 0.5N=1500N=20001,1,1,0,10,1,1,1,190%90%3,1,1,1,11,1,1,3,3:5 2N=1500N=20000,1,1,1,11,0,1,1,190%90%1,1,1,1,33,1,1,1,1n=2:2 0.5 2 N=1500N=20002,2,1,2,12,2,2,1,180%80%2,2,2,2,22,2,2,2,2:5 2 5N=1500N=20001,2,2,2,22,2,2,2,290%100%2,2,2,2,22,3,2,2,2由表 5 可以看出，对于增加惩罚的 SIC，不管是单变点和两个变点，正确识别的变点个数的正确率都是 80%以上；检测到的变点总数也是比较准确的，同时可以看出增加惩罚的 SIC 的检测效果是优于 AIC、SIC、mBIC 惩罚。3.2 实证分析根据 Maguire 等人13对英国 1875 至 1951 年期间发生地雷爆炸导致 10 人以上死亡的时间间隔这 109 个数据进行了研究，发现在 1890 年可能发生了数据变化。在 2004 年，Ramanayake 14对这组数据进行研究，发现在第 45 个数据发生了突变，同时验证数据服从伽马分布。针对这组数据，使用 AMOC、PELT、BS、WBS 这四种方法来进行验证，以验证方法的有效性。图 3 18751951 年英国地雷爆炸数据变化规律图由表 6 可以看出，在三种惩罚下 AMOC 方法能够准确识别数据中变点的位置和个数；PELT、BS 方法在 SIC、mBIC惩罚下同样能够对变点的位置和个数进行识别，但是在 AIC惩罚下会出现过估计的问题。WBS 方法在 mBIC 和增加惩罚的 SIC 能够较好地识别变点位置，但在 AIC 惩罚下也会出现过估计的情况。2023 年第 8 期78计算机应用信息技术与信息化表 6 AMOC、PELT、BS、WBS 方法对英国地雷爆炸数据变点位置的估计检测方法惩罚方式变点位置AMOCPELTBinary Segmenta-tionWild Binary Seg-mentationAICSICmBICAICSICmBICAICSICmBICAICSICmBIC增加惩罚的SIC46464622,24,53,73,80,102,104,107 464653,57,74,80,103,1094646103,74,46,75,53,108,78,57,77,80,54,59,61,46,67,69,87,98,49,1,38,73,90,79,91,24,56,20,40,26,94,19,97,86,85,76,33,83,51,106,4,28,14,50,70,63,5,42,105103,74,46,75,53,108,78,57,77,80,54103,46,75 103,46 4 结论本文采用AMOC、PELT、BS、WBS这四种变点检测方法，在三类不同的惩罚项下，对伽马分布中变点个数及位置进行了检测。模拟实验结果表明，在存在一个变点时，AMOC 方法的检测效果优于其他三种方法；而在存在两个变点时，采用 WBS 方法检测效果较好，且 mBIC 惩罚下的效果最佳。在 WBS 方法中，增加惩罚的 SIC 方法的检测效果要优于AIC、SIC、mBIC 惩罚。最后，将这四种检测方法应用于英国地雷爆炸灾害发生数据，检测结果表明，所检测的变点时间与当年实际发生变化的时间吻合，证明了这些方法的有效性。参考文献：1 STEINEBACH J.Change point and jump estimates in an AMOC renewal modelC/Asymptotic Statistics:Proceedings of the Fifth Prague Symposium.Physica-Verlag HD,1994:447-457.2 XU R,WU J,YUE X,et al.Online structural change-point detection of high-dimensional streaming data via dynamic sparse subspace learningJ.Technometrics,2023,65(1):19-32.3 李拂晓，陈占寿多元回归模型的结构变点估计 J.数学的实践与认识,2020,50(9):122-131.4 KORKAS K K,FRYZLEWICZ P.Multiple change-point detection for non-stationary time series using wild binary segmentationJ.Statistica sinica,2017,27(1):287-311.5 王平平,汤银才,程恭品.基于变点伽马过程的退化数据分析：贝叶斯角度（英文）J.应用概率统计,2022,38(5):674-692.6 WANG X,GUO P J B,CHENG Z.Real-time reliability evaluation for an individual product based on change-point gamma and wiener processJ.Quality&reliability engineering international,2014,30(4):513-525.7 KILLICK R,ECKLEY I.changepoint:An R package for changepoint analysisJ.Journal of statistical software,2014,58(3):1-19.8 KILLICK R,FEARNHEAD P,ECKLEY I A.Optimal detection of changepoints with a linear computational costJ.Journal of the american statistical association,2012,107:1590-1598.9 SEN A,SRIVASTAVA M S.On tests for detecting change in meanJ.The annals of statistics,1975,27(1):98-108.10FRYZLEWICZ P.Wild binary segmentation for multiple change-point detectionJ.The annals of statistics,2014,42(6):2243-2281.11AKAIKE H.A new look at the statistical model identificationJ.Automatic control,ieee transactions on,1974,19(6):16-723.12YAO Y C.Estimating the number of change-points via Schwarz criterionJ.Statistics&probability letters,1988,6(3):181-189.13ZHANG N R,SIEGMUND D O.A Modified Bayes Information Criterion with Applications to the Analysis of Comparative Genomic Hybridization DataJ.Biometrics,2007,63,22-32.14MAGUIRE B A,et al.The time intervals between industrial accidentsJ.Biometrika,1952,39:168-180.15RAMANAYAKE A.Tests for a Change Point in the Shape Parameter of Gamma Random VariablesJ.Communications in statistics,2005,33(4):821-833.【作者简介】王 轩 云(1998)，通 信 作 者（邮 箱：），女，贵州威宁人，硕士研究生，研究方向：统计模型与统计计算。黄介武(1977)，男，湖南沅江人，教授，博士，研究方向：统计模型与统计计算。（收稿日期：2023-02-28 修回日期：2023-04-26）