基于记录值的逆威布尔分布参数统计推断
Statistical Inference on Parameters of Inverse Weibull Distribution Based on Record Values
摘要: 本文基于下记录值样本,研究了两参数逆威布尔分布模型的点估计和区间估计问题。基于构造的服从卡方分布的枢轴量导出了参数的逆矩估计,同时得到了形状参数的传统等尾置信区间和最短置信区间以及形状参数和尺度参数的联合置信域,最后通过数值实例将得到的区间和已有文献基于构造的服从F分布的枢轴量得到的区间进行比较以说明本文方法的优良性。
Abstract: In this paper, point estimation and confidence intervals estimation of inverse Weibull distribution model with two parameters are studied based on the sample of lower record values. The inverse moment estimator of parameters is derived based on the structure of pivotal quantities subject to the chi-square distribution, and at the same time, the traditional equal tailed confidence interval and the shortest confidence interval of shape parameter, as well as the joint confidence field of shape and scale parameters, are obtained. Finally, a numerical example is given to compare the obtained intervals with the obtained intervals based on the structure of pivot quantities subject to F distribution in the existing literature, so as to illustrate the superiority of the method in this paper.
文章引用:段李永, 于秀君, 王涛. 基于记录值的逆威布尔分布参数统计推断[J]. 统计学与应用, 2022, 11(2): 346-353. https://doi.org/10.12677/SA.2022.112036

1. 引言

记录值作为一种刻画数据变化趋势的数据,自文献 [1] 提出以来已被广泛应用到天气预报、体育赛事、水文地震等诸多领域。如奥运会110米跨栏比赛的用时记录值;某些干旱地区年降雨量最小值变化的记录值等。因此,对记录值数据进行理论研究具有较大的实际应用价值。近些年,已有大量文献研究了记录值样本下参数的统计推断问题,如文献 [2] [3] 定义了记录值,并讨论了它们相关的统计性质;文献 [4] 基于构造的枢轴量研究了记录值样本下威布尔分布的参数估计以及应力强度可靠性统计推断问题;文献 [5] 研究了记录值样本下极值分布的参数点估计,区间估计以及联合置信域。

逆威布尔分布作为一种常见的寿命试验分布之一,最早由凯勒和卡玛特在文献 [6] 中对可靠性数据建模时提出,并且通过研究得出一些机械部件失效是由于其退化引起的结果。它现已在寿命试验、可靠性工程、电子设备等领域受到广泛关注和研究。文献 [7] 基于构造的服从F分布的枢轴量探讨了甘贝尔分布和逆威布尔分布参数的精确置信区间和置信域;文献 [8] 基于下记录值样本值,在形状参数已知、形状和尺度参数都未知的两种情况下研究了参数的不同估计方法。文献 [9] 基于记录值样本,讨论了逆威布尔分布的两个参数的最大似然估计和贝叶斯估计及预测问题。更多的关于记录值样本下参数的统计推断问题见文献 [10] [11] [12]。

本文将利用文献 [4] 的理论构造一个服从卡方分布的枢轴量,进一步探讨逆威布尔分布的参数估计问题,并把得到的形状参数精确置信区间和形状和尺度参数的联合置信域结果与文献 [7] 的相应的置信区间和联合置信域作比较,最后再把文献 [13] 中寻求最短置信区间方法用到该分布形状参数的置信区间的推断中,并把得到的最短置信区间与本文得到的传统等尾置信区间、文献 [7] 中的置信区间作比较。从比较结果看,本文得到的置信区间的长度确实比文献 [7] 中的要短。

设随机变量X服从逆威布尔分布,其概率密度函数如下

f ( x ; η , β ) = β η β x ( β + 1 ) e ( η x ) β , x > 0 (1)

相应的概率分布函数为

F ( x ) = e ( η x ) β , x > 0 (2)

其中 η , β > 0 ,分别为尺度和形状参数。

2. 最大似然估计

定义1 [1] 设 { X n , n 1 } 是一个独立同分布的随机变量序列,对于任意的 i > j ,如果有 x i < x j 成立,则称 x i 为一个下记录值。令 L 1 = 1 ; L n = min { i : X i < X L n 1 } , n 2 ,则称子序列 { Y n = X L n , n 1 } 为原始序列 { X n , n 1 } 的一个下记录值序列,而称 L n 为记录时间序列,类似可定义上记录值。

Y = { Y 1 , Y 2 , , Y n } 是一个来自逆威布尔分布(1)的n个下记录值序列, y = { y 1 , y 2 , , y n } 是一个样本观测值,则由文献 [3] 知样本似然函数为:

L ( η , β ) = f ( y n ) i = 1 n 1 f ( y i ) F ( y i ) = β n η n β e ( η y n ) β i = 1 n y i ( β + 1 ) (3)

而对数似然函数为:

ln L = n ln β + n β ln η ( η y n ) β ( β + 1 ) i = 1 n ln y i (4)

通过最大化(4)式

ln L β = n β + n ln η ( η y n ) β ln ( η y n ) i = 1 n ln y i = 0 (5)

ln L η = n β η β η β 1 y n β = 0 (6)

可以得到参数 η , β 的最大似然估计分别为:

β ^ = n i = 1 n 1 ln ( y i / y n )

η ^ = n 1 / β ^ y n

3. 参数的置信区间和联合置信域

Y 1 > Y 2 > > Y n 是来自逆威布尔分布(1)的n个下记录值,令 Z i = ln [ F ( Y i ; η , β ) ] = ( η Y i ) β ,则

Z 1 < Z 2 < < Z n 是来自标准指数分布的n个上记录值。再令 T 1 = Z 1 ; T 1 = Z i Z i 1 , i = 1 , 2 , , n 。则由文献 [4] 知 T 1 , T 2 , , T n 是服从标准指数分布相互独立的随机变量。

引理1 [4] 设 T 1 , T 2 , , T n 是相互独立同分布的标准指数随机变量。令 S i = T 1 + T 2 + + T i U i = ( S i S i + 1 ) i i = 1 , 2 , , n 1 U n = S n 则有

1) U 1 , U 2 , , U n 相互独立;

2) U 1 , U 2 , , U n 1 服从均匀分布 U ( 0 , 1 )

3) 2 U n 服从自由度为2n的卡方分布。

定理1 令 W ( β ) = 2 i = 1 n 1 ln ( U i ) = 2 β i = 1 n 1 ln ( Y i / Y n ) ,当 n > 1 时,有下面的结论成立。

1) W ( β ) 服从自由度为 2 n 2 的卡方分布。

2) 对任意的 t > 0 W ( β ) = t 关于当 β 有唯一解。

证明:1) 由 U 1 , U 2 , , U n 1 服从均匀分布 U ( 0 , 1 ) ,可知 2 ln ( U i ) ~ χ 2 ( 2 ) ,又 U 1 , U 2 , , U n 1 相互独立,所以 W ( β ) = 2 i = 1 n 1 ln ( U i ) = 2 β i = 1 n 1 ln ( Y i / Y n ) ~ χ 2 ( 2 n 2 )

2) 显然 W ( β ) 关于 β 是严格单调递增的,再由 lim β 0 W ( β ) = 0 lim β W ( β ) = 可知,对任意的 t > 0 W ( β ) = t 关于当 β 有唯一解。

定理2 设 Y 1 > Y 2 > > Y n 是来自逆威布尔分布(1)的n个下记录值,则参数 η , β 的逆矩估计量由下式确定

{ W ( β ) = 2 n 2 2 U n = 2 ( n Y n ) β = 2 n

定理3设 Y 1 > Y 2 > > Y n 是来自逆威布尔分布(1)的n个下记录值,则参数 β 1 α ( 0 < α < 1 ) 水平的置信区间为

( χ α / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) , χ 1 α / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) ) (7)

证明:从引理2可知枢轴量

W ( β ) = 2 β i = 1 n 1 ln ( Y i / Y n ) ~ χ 2 ( 2 n 2 )

从而有,

P ( χ α / 2 2 ( 2 n 2 ) < W ( β ) < χ 1 α / 2 2 ( 2 n 2 ) ) = 1 α

P ( χ α / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) < β < χ 1 α / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) ) = 1 α

定理4设 Y 1 > Y 2 > > Y n 是来自逆威布尔分布(1)的n个下记录值,则参数 η , β 1 α ( 0 < α < 1 ) 联合置信域为

{ χ ( 1 1 α ) / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) < β < χ ( 1 + 1 α ) / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) ( χ ( 1 1 α ) / 2 2 ( 2 n ) 2 ) 1 / β Y n < η < ( χ ( 1 + 1 α ) / 2 2 ( 2 n ) 2 ) 1 / β Y n

证明:因为 2 U n ~ χ 2 ( 2 n ) ,且由文献 [4] 知 W ( β ) 2 U n 是相互独立的。所以对于任意的 0 < α < 1 ,有

P ( χ ( 1 1 α ) / 2 2 ( 2 n ) < 2 U n < χ ( 1 + 1 α ) / 2 2 ( 2 n ) ) = 1 α

由定理2知

P ( χ ( 1 1 α ) / 2 2 ( 2 n 2 ) < W ( β ) < χ ( 1 + 1 α ) / 2 2 ( 2 n 2 ) ) = 1 α

于是有

P ( χ ( 1 1 α ) / 2 2 ( 2 n 2 ) < W ( β ) < χ ( 1 + 1 α ) / 2 2 ( 2 n 2 ) , χ ( 1 1 α ) / 2 2 ( 2 n ) < 2 U n < χ ( 1 + 1 α ) / 2 2 ( 2 n ) ) = 1 α

P ( χ ( 1 1 α ) / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) < β < χ ( 1 + 1 α ) / 2 2 ( 2 n 2 ) 2 i = 1 n 1 ln ( Y i / Y n ) , ( χ ( 1 1 α ) / 2 2 ( 2 n ) 2 ) 1 / β Y n < η < ( χ ( 1 + 1 α ) / 2 2 ( 2 n ) 2 ) 1 / β Y n ) = 1 α

4. 形状参数的最短置信区间

一般对参数作区间估计时,在固定的置信水平下,所寻求的置信区间应该越短越好,由于本文构造的枢轴量的分布服从单峰的卡方分布,它的密度函数关于其峰值是不对称的,因此,传统的等尾置信区间,即(7)式确定的置信区间不是最短的,下面在固定的置信水平 1 α 下寻求参数 β 的最短置信区间,所谓最短,就是要在满足 P ( b < η < c ) = 1 α 的所有置信区间 ( b , c ) 中,找到使得 c b 最小的区间 ( b , c ) 。为了方便,采用 1 2 W ( β ) 这个枢轴量来讨论参数 β 的最短置信区间。由于

P ( b < 1 2 W ( β ) < c ) = P ( W 1 ( 2 b ) < β < W 1 ( 2 c ) ) = 1 α ,所以本文要找的参数 β 的最短置信区间 ( W 1 ( 2 b ) , W 1 ( 2 c ) ) 要满足:

{ min L = W 1 ( 2 c ) W 1 ( 2 b ) s .t . P ( b < 1 2 W ( β ) < c ) = 1 α (8)

这里 W 1 ( t ) 表示 W ( β ) = t 关于 β 的解。

定理5 当 n > 2 时,上述条件极值(8)存在唯一驻点 ( b * , c * )

证明:由于 W ( β ) β 的线性单调函数,故使得枢轴量 1 2 W ( β ) 的置信区间达到最短的b和c,也使得参数 β 的置信区间最短。事实上,由 W ( β ) = 2 β i = 1 n 1 ln ( Y i / Y n ) 可知, L = W 1 ( 2 c ) W 1 ( 2 b ) = c b i = 1 n 1 ln ( Y i / Y n ) ,而 i = 1 n 1 ln ( Y i / Y n ) 只与样本有关,因此要使得L最小,只需 c b 最小,从而只需证明在 P ( b < 1 2 W ( β ) < c ) = 1 α 的约束下, L = c b 存在唯一驻点 ( b * , c * ) 使得本身最小即可。约束条件 P ( b < 1 2 W ( β ) < c ) = 1 α 等价于

{ 0 b g ( y ) d y = β 0 c g ( y ) d y = 1 α + β (9)

这里 g ( y ) = y n 2 e y Γ ( n 1 ) , y > 0 是形状参数为 n 1 ,尺度参数为1的伽马分布密度函数,且 0 < β < α

我们记相应的概率分布函数为 G ( y ) 。由(9)式可知,在固定置信水平 1 α 下,任意给定一个 β ,就有唯一的一个b和c与之对应,即b和c都可以看作是 β 的函数,从而 L = c b 也是 β 的函数,故有

d L d β = d c d β d b d β = 1 g ( c ) 1 g ( b )

d L d β = 0 ,可得 g ( b ) = g ( c ) ,即

b n 2 e b = c n 2 e c (10)

现在要求的唯一驻点 ( b , c ) 就是在 b n 2 e b = c n 2 e c 的约束下证明 P ( b < 1 2 W ( β ) < c ) = b c g ( y ) d y = G ( c ) G ( b ) = 1 α 的解唯一。

考虑函数

h ( x ) = x n 2 e x , x > 0

h ( x ) = 0 ,可得到 h ( x ) 的最大值点 x 0 = n 2 ;且当 0 < x x 0 时, h ( x ) 严格单调递增, x x 0 时, h ( x ) 严格单调递减,而当x趋于无穷大或者零时 h ( x ) 趋于零,为确保 b < c 和(10)成立,应有 b < n 2 c > n 2 ,这样对任意的c,从(10)就可以得到唯一的 b = u ( c )

可以知道的是函数 g ( y ) 的峰值也在 x = x 0 = n 2 处达到。由函数 h ( x ) 的单调性可知,要满足 b n 2 e b = c n 2 e c ,则当c增大时, b = u ( c ) 必然要减小。所以当 c > n 2 时, G ( c ) G ( u ( c ) ) 是c的严格增函数,而且当c趋于无穷大时, b = u ( c ) 从大于零的一侧趋于零;当c从大于 n 2 的一侧趋于 n 2 时, b = u ( c ) 从小于 n 2 的一侧趋于 n 2 ,因此有

lim c [ G ( c ) G ( u ( c ) ) ] = 1

lim c n 2 [ G ( c ) G ( u ( c ) ) ] = 0

由中值定理知,在 b n 2 e b = c n 2 e c 的约束下, G ( c ) G ( u ( c ) ) = 1 α ( 0 < α < 1 ) 的解是唯一的,即存在唯一的 c = c * 使得 G ( c * ) G ( u ( c * ) ) = 1 α ( 0 < α < 1 ) 成立,再取 b = u ( c * ) ,则 ( b * , c * ) 是条件极值(8)的唯一驻点,命题得证。下面将通过实例给出怎样求最短置信区间,并把得到的最短置信区间与传统等尾置信区间作比较。

5. 实例应用

本文采用文献 [7] 中5.3例子的样本,数据为一组来自宾夕法尼亚州哈里斯堡的萨斯奎赫纳河的最高洪水水位(每秒数百万立方英尺)且样本容量为 n = 20 ,而文献 [14] 已经证明了用逆威布尔分布来拟合该组数据的效果很好,数据见表1

Table 1. Maximum flood level data

表1. 最高洪水水位数据

可以得到该组样本的6个下记录值为: 0.654 , 0.613 , 0.315 , 0.297 , 0.269 , 0.265 ,于是有:

参数 η , β 的最大似然估计分别为: η ^ = 0.4879 , β ^ = 2.9357

参数 η , β 的逆矩估计分别为: η = 0.5512 , β = 2.4464

参数 β 的置信水平为95%的传统等尾置信区间为: ( 0.7943 , 5.0110 )

为了得到参数 η , β 的置信水平为95%的联合置信域,需要下面的分位数

χ 0.0127 2 ( 10 ) = 2.7182 ; χ 0.9873 2 ( 10 ) = 22.5120

χ 0.0127 2 ( 12 ) = 3.7659 ; χ 0.9873 2 ( 12 ) = 25.4812

由定理4可知参数 η , β 的置信水平为95%的联合置信域为:

{ 0.6650 < β < 5.5072 0.265 ( 3.7659 2 ) 1 / β < η < 0.265 ( 25.4812 2 ) 1 / β

文献[7]中参数 λ ( = η β ,也就是文献中 [7] 的参数 λ 相当于本文的 η β ), β 的置信水平为95%的联合置信域的面积为3.6918。而本文的参数 λ ( = η β ) , β 的置信水平为95%的联合置信域的面积为:

{ 0.6650 < β < 5.5072 3.7659 2 0.265 β < η < 25.4812 2 0.265 β

即联合置信域的面积为3.3751,可见比文献 [7] 中的参数 λ ( = η β ) , β 的联合置信域3.6918要小。

由定理4的证明过程可知,参数 β 的置信水平为95%的最短置信区间可以通过求解下面非线性规划问题:

{ min L = c b e b k = 0 n 2 b k k ! e c k = 0 n 2 c k k ! = 1 α b n 2 e b = c n 2 e c

得到唯一驻点 ( b * , c * ) = ( 1.2070 , 9.4302 ) ,从而参数 β 的置信水平为95%的最短置信区间就为 ( W 1 ( 2 b * ) , W 1 ( 2 c * ) ) = ( 0.5906 , 4.6140 )

文献 [7] 中得到的参数 β 的置信水平为95%的置信区间和本文得到的置信区间及其长度L比较见表2

Table 2. Length comparison of confidence intervals

表2. 置信区间的长度比较

6. 结论

本文基于下记录值样本,首先给出了逆威布尔分布参数 η , β 的最大似然估计,然后利用文献 [4] 的引理构造了服从卡方分布的枢轴量,并基于该枢轴量得到了这两个参数的逆矩估计,形状参数 β 的置信区间,以及参数 η , β 的联合置信域,接着基于文献 [13] 的结果获得了形状参数 β 的最短置信区间,最后基于文献 [7] 的数值实例得到了相关的具体结果,并把它们与文献 [7] 的结果作比较。从比较结果可以看出,在相同的置信水平下,本文基于构造的服从卡方分布的枢轴量得到的置信区间和联合置信域比文献 [7] 基于构造的服从F分布的枢轴量得到的结果要短、小,因此是一个可供相对更好的选择方式。另外,还可以看到在固定置信水平下,最短置信区间的长度确实比文献 [7] 中的结果和本文得到的传统等尾置信区间的长度都要短。

NOTES

*通讯作者。

参考文献

[1] Chandler, K.N. (1952) The Distribution and Frequency of Record Values. Journal of the Royal Statistical Society: Series B (Methodological), 14, 220-228.
https://doi.org/10.1111/j.2517-6161.1952.tb00115.x
[2] Ahsanullah, M. (1995) Record Statistics. Nova Science Publishers Commack, New York.
[3] Arnold, B.C., Balakrishnan, N. and Nagaraja, H.N. (1998) Records. Wiley, New York.
https://doi.org/10.1002/9781118150412
[4] Wang, B.X. and Ye, Z.-S. (2015) Inference on the Weibull Distribution Based on Record Values. Computational Statistics and Data Analysis, 83, 26-36.
https://doi.org/10.1016/j.csda.2014.09.005
[5] 罗嘉成, 陈勇明. 记录值样本下极值分布的统计推断[J]. 统计与决策, 2017(18): 10-12.
[6] Keller, A.Z., Kamath, A.R.R. and Perera, U.D. (1982) Reliability Analysis of CNC Machine Tools. Reliability Engineering, 3, 449-473.
https://doi.org/10.1016/0143-8174(82)90036-1
[7] Asgharzadeh, A., Abdi, M. and Nadarajah, S. (2016) Interval Estimation for Gumbel Distribution Using Climate Records. Bulletin of the Malaysian Mathematical Sciences Society, 39, 257-270.
https://doi.org/10.1007/s40840-015-0185-2
[8] Sultan. K.S. (2010) Record Values from the Inverse Weibull Lifetime Model: Different Methods of Estimation. Intelligent Information Management, 2, 631-636.
https://doi.org/10.4236/iim.2010.211072
[9] El-Din, M.M.M., Riad, F.H. and El-Sayed, M.A. (2014) Statistical Inference and Prediction for the Inverse Weibull Distribution Based On Record Data. Journal of Statistics Applications & Probability, 3, 171-177.
https://doi.org/10.12785/jsap/030207
[10] Lee, M.Y. (2018) On Characterizations of the Inverse Weibull Distribution Based on Record Values. Journal of Applied Mathematics & Informatics, 36, 429-433.
[11] Pak, A. and Dey, S. (2019) Statistical Inference for the Power Lindley Model Based on Record Values and Inter-Record Times. Journal of Computational and Applied Mathematics, 347, 156-172.
https://doi.org/10.1016/j.cam.2018.08.012
[12] Pak, A., Jafari, A.A. and Khoolenjani, N.B. (2020) On Reliability in a Multicomponent Stress-Strength Generalized Rayleigh Model Based on Record Values. Journal of Testing and Evaluation, 48, No. 6.
https://doi.org/10.1520/JTE20170732
[13] 蒋福坤, 刘正春. 指数分布参数的最短区间估计[J]. 数理统计与管理, 2004, 23(3): 43-45+10.
[14] Maswadah, M. (2003) Conditional Confidence Interval Estimation for the Inverse Weibull Distribution Based on Censored Generalized Order Statistics. Journal of Statistical Computation and Simulation, 73, 887-898.
https://doi.org/10.1080/0094965031000099140