广义Pareto 分布参数的Bootstrap置信区间及应用

张艳芳, 赵宜宾, 任晴晴

(1.防灾科技学院，河北 三河 065201; 2.河北省地震动力学重点实验室，河北 三河 065201)

0 引言

极值理论是数学在工程、环境及风险管理问题中应用非常广泛的统计方法。由阈值模型导出的广义Pareto 分布(Generalized Pareto Distribution, GPD)能充分利用数据中包含的极值信息而被广泛应用于各行业数据的分析[1–3]。在广义Pareto 分布的应用中，参数的估计和分位数的计算非常重要。很多学者从不同角度改进了传统的广义Pareto 模型以便得到更好的参数估计。胡尧等[4]建立了GPD 变点模型，能够较好地捕获传统GPD 模型所不能捕捉的内在规律，弥补了传统GPD 模型的不足。Nagatsuka 和Balakrishnan[5]建立了一种任何形状参数都适用的广义Pareto 分布方法。Sauer 等[6]基于广义Pareto 分布的逐步II 型截尾样本的极大似然估计研究了系统的可靠性。Zhao 等[7]应用广义Pareto 分布的阈值峰值框架提出了一种新的方法，其利用部分L-矩导出了GPD 参数计算的有效估计量。还有学者将其它统计方法与广义Pareto 模型相结合分析问题[8–10]。

在用GPD 模型拟合数据时，阈值u的选取一直是一个棘手的问题[11]，u值太大只有很少几个超出量，估计量的方差就较大；u值太小，超出量的分布与广义Pareto 分布相差较大，估计量成为有偏估计。因此在实际数据处理中，一方面为了选取足够大的阈值确保超出量的分布近似为GPD 分布，而这时通常超过阈值的数据太少，导致方差太大，很难得到有意义的推断。另一方面，Pareto 分布中参数的区间估计是使用Fisher 信息矩阵得到参数的渐近正态协方差矩阵，从而计算出参数的近似置信区间[12]。基于以上原因，本文在确定阈值和分布参数的极大似然估计后，应用Bootstrap 方法求出参数的Bootstrap 置信区间估计。文中对巴颜喀拉地震带震级数据进行分析获得了相应的参数估计，并且应用Bootstrap 置信区间方法比应用Fisher 信息矩阵方法得到的参数的近似置信区间更短。本文分为四个部分，第1 部分介绍广义Pareto 分布的基本概念及参数估计方法；第2 部分介绍广义Pareto 分布的Bootstrap 置信区间估计构造方法；第3 部分将上述方法应用于分析巴颜喀拉地震带震级数据分析中；第4 部分是模型讨论与展望。

1 广义Pareto 分布

1.1 广义Pareto 分布的基本概念

Pickands[13]于1975 年首次提出了广义Pareto 分布，证明了若选定一个较大的阈值，则取自总体的样本中超过阈值的样本渐近地服从GPD。该理论一经提出便被广泛应用于风险评估、寿命分析、地震预测等领域。本文研究的是双参数广义Pareto 分布，以下按照文献[13]给出相关的定义。

若随机变量X的分布函数为

其中称ξ ∈R为形状参数，σ> 0 为比例参数，且x ≥0, 1+ξx/σ> 0，则称X服从广义Pareto 分布。

定理1 假设X1,X2,···,Xn独立且服从广义Pareto 分布，则对于固定的阈值u，超阈值Y=Xi-u，近似服从双参数的广义Pareto 分布GP(y,ξ,˜σ)，且两者具有相同的形状参数ξ。

1.2 广义Pareto 分布的参数估计

参数估计有多种方法，本文采用极大似然估计，假设Y1,Y2,···,Yn服从广义Pareto 分布，不失一般性，设ξ ̸=0，由概率密度函数

上述似然函数取得最大值的参数没有解析解，通过数值解法可以得到ξ和σ的极大似然估计的数值解ˆξ和ˆσ。

令信息矩阵为IE(ξ,σ)，当n →∞时，由观测信息矩阵的逆矩阵可得到ξ和σ的极大似然估计ˆξ和ˆσ的渐近正态协方差矩阵[12]，进而求得各参数估计的近似置信区间

1.3 广义Pareto 分布的分位数

在实际应用中，极值分析非常重要的指标重现水平是由分位数xp确定的，以下给出分位数的估计方法。如果X ∼GP(x,ξ,σ)，则它的尾部分布为

上述的公式也可以转化为地震危险性的评价指标。

指标1 重现期：对于给定的时间序列X1,X2,···,Xn，及阈值u，序列中变量相邻两次超阈值平均时间间隔T(u)称为重现期：

2 广义Pareto 分布参数的Bootstrap 置信区间估计

2.1 广义Pareto 分布参数极大似然估计的Bootstrap 分位数置信区间估计

Bootstrap 方法是Efron 在20 世纪70 年代后期建立的。本文采用参数Bootstrap 方法求广义Pareto 分布两个参数的区间估计。具体步骤如下：

1) 对于给定的来自广义Pareto 分布GP(x,ξ,σ)的一组样本X1,X2,···,Xn，x1,x2,···,xn是一组已知的样本值；

2) 由1.2 节(4)式和1.3 节(9)式，应用数值解法可以求得参数ξ和σ的极大似然估计值ˆξ和ˆσ；

此方法将传统的极大似然估计与Bootstrap 置信区间相结合得到参数的Bootstrap 置信区间。本文采用的是参数的Bootstrap 分位数区间估计，这类型的区间估计具有一定的优良性质。

2.2 广义Pareto 分布参数极大似然估计Bootstrap 分位数置信区间的性质

2.2.1 Bootstrap 分位数置信区间的性质

假设总体X →F(x,θ)，总体的样本X1,X2,···,Xn。未知参数θ的极大似然估计为ˆθ，针对B个Bootstrap 样本计算极大似然估计，则参数θ的Bootstrap 分位数置信区间定义为

以上结论表明，在引理1 成立的情况下，Bootstrap 分位数置信区间与精确的置信区间是一致的。如果引理1 中(14)式对较大的n近似成立，则Bootstrap 分位数置信区间的也是渐近有效的。

2.2.2 数值模拟研究Bootstrap 分位数置信区间的表现

为研究置信区间的优劣，应用数值模拟计算置信区间的覆盖率和区间长度，其中覆盖率记为

根据文献[16]可知，一个较好的置信区间应该满足：覆盖率接近预先设定的置信水平；置信区间长度短。

数值模拟中研究以下情形：ξ=1.5,σ=0.5,α=0.05，样本量n=30,50,100,B=100,200,500。表1 给出了这些模拟中置信区间覆盖率和相应的区间长度。

表1 置信区间覆盖率和相应的区间长度

由表1 可以看出，在选取两组不同的参数情况下，Bootstrap 分位数置信区间具有以下特点：在样本量n固定时，随着抽样次数的不断增大所得的置信区间覆盖率基本变化趋势是不断增大，区间长度减小；另一方面，随着样本量不断增加，置信区间覆盖率基本也是增加，且置信区间平均长度减小。模拟结果表明，Bootstrap 分位数置信区间在样本较小的情形下(如n= 30 时)，抽样的次数B 较大的时候，置信区间的覆盖率比较接近给定的置信水平0.95，随着样本容量的增加表现相对稳定，与极大似然估计结合给定参数的置信区间是可行的。

以下将第1 节和第2 节的分布拟合和参数估计应用于巴颜喀拉地震带震级数据分析，进行地震危险性分析。

3 巴颜喀拉地震数据分析

本文使用的数据来源于中国地震信息网提供的测震目录，从1920∼2019 年的巴颜喀拉断裂带的245 384 个震级数据，地震震级采用面波震级Ms。4 级以上震例从1950 年开始丰富，2 级以上震例从1965 年记录完整，1965 年之前较小的震级没有被记录或没有被观测到，为了尽量保留历史数据提供的信息，可以用阈值法建立一个近似的极值模型。

3.1 阈值选取

设定一个阈值，把超过这个阈值的事件作为极端事件。这里根据平均剩余寿命图选择阈值，图1 是震级的平均剩余寿命图及相应的95%置信区间。由图1 可以看到，从u=0 直到u≈6 图形近似呈直线，超过u=6，图形急剧下降，这表明应该取u=6。

图1 震级平均剩余寿命图

图2 和图3 是巴颜喀拉地震带震级超出量拟合的广义Pareto 分布，在不同阈值选取下参数的极大似然估计值变化图。图中形状参数和修正的尺度参数在区间[4,6.5]基本保持不变，确定可以选取阈值u=6。

图2 阈值–修正的尺度参数

图3 阈值–形状参数图

为了达到定理1 的要求，不降低阈值，这里选取阈值为u= 6。图4 为确定阈值后的数据样本散点图，也是下文中拟合模型的数据。当u= 6 时，共有79 个数据超过阈值，超过阈值的数据占总量的79/245 384 = 3.219×10-4，相对整个数据比例较少，应用Bootstrap 方法构建置信区间能充分应用样本信息。

图4 超过阈值6 的震级散点图

3.2 广义Pareto 分布模型的拟合与检验

依据1.2 节的广义Pareto 分布参数估计理论，通过数值计算可求得参数的极大似然估计值ˆξ=-0.356, ˆσ= 0.884，相应的对数似然函数值为41.063。按(4)式可计算出各参数95%的置信区间，如表2 所示。参数估计ˆξ< 0，从表2 可以看到，ξ的整个置信区间，ξ总为负表示分布支撑上端点是有限的，也就是对应震级有上限，数据拟合广义Pareto 分布是合理的。

表2 广义Pareto 分布参数的极大似然估计

采用第3 节参数的Bootstrap 置信区间估计方法，取B= 5 000,α= 0.05，计算得到ξ和σ的Bootstrap 置信区间估计，如表3 所示。

表3 广义Pareto 分布参数的置信区间比较

由表3 可以看出，在相同的置信度0.05 下，采用Bootstrap 方法计算的置信区间相比由渐近正态协方差矩阵所计算的区间长度更短。Bootstrap 置信区间计算方法可以充分利用给定的样本信息，参数ξ和σ的置信区间为参数估计的不确定提供了定量验证的依据。

模型检验结果如图5 所示。

图5 巴颜喀拉地震带震级拟合检验图

由图5 可以看出，分布的理论密度曲线与震级直方图趋势基本吻合，广义Pareto 分布在震级越高的段位出现的偏离越小，这是由于广义Pareto 分布理论主要是针对尾部数据建模，模型的适用性在于强震段。P-P 图和Q-Q 图的散点总体在一条直线附近小幅波动，说明理论分位数值与观测值分位数几乎是一致的。重现水平曲线渐近的趋于某个有限值，由(7)式可以计算出震级的最大值为Mmax= 8.48，与该地震带历史最大震级8.1 级相吻合。根据(12)式重现期为100 年的重现水平为7.45，即7.5 级左右的超强震约100 年一遇。由(10)式计算出xp的置信区间为[7.23,7.66]，百年重现水平的区间估计上限为7.66，说明巴颜喀拉地震带百年里有可能发生强震。重现水平的置信区间说明数据与模型的偏离不大，模型拟合较合理。

4 模型讨论与展望

本文将广义Pareto 分布模型应用于巴颜喀拉地震带强震震级的分析中，从分布拟合检验结果中可以看出，广义Pareto 分布对于超阈值数据拟合具有明显的优势。合理的分布模型更需要精确的参数估计支撑模型的确定性，而传统极值分析中超阈值统计数据量少使得分析结果不确定性增加。文中针对这种情况，应用Bootstrap 抽样方法，通过重复抽样，生成大量子样本以得到更精确的参数的区间估计。Bootstrap 抽样方法计算简便，结果表明在巴颜喀拉地震带震级拟合中，置信度相同的情况下，相比传统的应用近似协方差矩阵计算所得置信区间长度更短。同时，Bootstrap 方法的应用降低了由于样本变异引起的模型的不确定性。文中的分析过程也可应用于其它行业数据的尾部分析中。