0 引言

寿命分布的研究在科学技术领域中占有重要地位，随着寿命分布理论的逐步发展，人们陆续提出了多种产生寿命分布的机制，为寿命数据的研究提供了便利。1979年，马逢时和刘德辅提出复合极值分布概念^[1]，将一个连续分布与一个离散分布混合得到一个新的连续型分布，这种产生分布的方法在极值事件的研究中受到重视，其基本思想是设ξ₁, ⋯, ξ_z, ⋯为独立随机变量序列且有相同的分布函数G (x)，又设Z是与{ ξ_z, z =1, 2, ⋯ }相互独立的离散随机变量且其分布律为P_z=P（Z=z），记X=min{ξ₁, ξ₂, ⋯, ξ_z}，则称X服从由G (x)和Pz复合而成的分布，它表示串联系统的寿命分布。由于故对于极大值的研究与极小值相似。根据这种机制，人们提出了很多新型的寿命分布。Kuş^[2]（2007年）将指数分布与泊松分布复合，提出两参数的Exponential-Possion分布；随后，Tahmasbi和Rezaei^[3]（2008年）将指数分布与Logarithmic分布复合，提出两参数的Exponential-Logarithmic分布。在指数分布与离散分布复合取得进展后，人们开始转向更多连续分布与离散分布复合的探索。Lu等^[4]（2012年）将Weibull分布和Possion分布复合，提出三参数的Weibull-Possion分布，并进行了实际应用；Zakerzadeh和Mahmoudi^[5]（2012年）提出Lindley-Geometric分布，并表示Lindley分布为该分布的一种特殊形式；Gui等^[6]（2014年）将Lindley分布和Possion分布复合，提出Lindley-Possion分布；吕晓星等^[7]（2015年）将Pareto分布和Possion分布复合，提出Pareto-Logarithmic分布；Liyanage等^[8]（2015年）提出Lindley-Logarithmic分布，并与Exponential-Possion分布、Exponential-Logarithmic分布等进行了比较，表明了该分布的良好适用性。这些文献对所提出的分布进行了研究，得到其各阶矩、熵、参数的极大似然估计（maximum likelihood estimation，MLE）及参数估计的EM算法等。本文在研究新分布的性质中，借鉴了前人的思路，同时由于极大似然估计形式的复杂性，也采用了便于求解的EM算法。

1958年Lindley^[9]提出一种以权重和将Exp(θ)分布与Gamma (2，θ)分布混合得到的分布，称为参数为θ的Lindley分布，其密度函数为

Ghitany等^[10]（2008年）对Lindley分布进行了详细地研究，发现多数情况下Lindley分布具有比指数分布更好的性质，在生存分析与可靠性的研究中意义深远。Ghitany等^[11]（2013年）对Lindley分布进行推广得到两参数的Power Lindley分布，其密度函数为

(1)

和Possion分布、几何分布一样，Logarithmic分布是一种离散计数分布，该分布来源于-ln(1 - p)的麦克劳林展开即-ln(1 - p)= p + p² /2 + p³ /3 + ⋯。Johnson等^[12]详细研究了Logarithmic分布作为计数分布的重要性以及其各种性质。该分布在植物以及生物学领域研究中有广泛应用。

本文将Power Lindley（PL）分布与Logarithmic分布进行复合得到一个危险率形式多样的新的三参数寿命分布：Power Lindley-Logarithmic（PLL）分布，研究其相关性质，给出了参数的极大似然估计，并验证了估计的相合性和渐近正态性，给出了求解极大似然估计的EM算法，并进行了Monte Carlo数值模拟。

1 Power Lindley-Logarithmic分布的定义

设ξ₁，⋯，ξ_z，⋯为独立同分布的非负随机变量序列，其共同分布是参数为(θ，α)的Power Lindley分布，记为PL(θ，α)，其密度函数如（1）式所述，分布函数为

(2)

又设Z是与{ ξ_z，z = 1，2，⋯}相互独立的随机变量，服从参数为p ∈(0，1)的Logarithmic分布，其分布律为

(3)

令X = min { ξ₁，ξ₂⋯，ξ_Z}，则易见X的分布函数为

(4)

称此分布为以(p，θ，α)为参数的Power Lindley-Logarithmic分布，记为PLL(p，θ，α)，易见其密度函数为

(5)

特别地，当α = 1时，PLL(p，θ，α)就是Liyanage^[8]等提出的参数为p，θ的Lindley-Logarithmic分布，简记为LL(p，θ)分布。

2 Power Lindley-Logarithmic分布的性质

由（5）式，容易得出：当α > 1且θ ≥ 2α - 1或α = 1且p ≥ θ²时，PLL(p，θ，α)为单峰分布；当1且或或α = 1且p ≤ θ²时，PLL(p，θ，α)分布的函数单调递减。图 1给出f (x；p，θ，α)在参数取某些特定值时PLL(p，θ，α)分布的密度函数。

定理1   设0 < u < 1. PLL(p，θ，α)分布的u-分位数为

其中W_-1(·)表示Lambert W函数的负数部分，即方程的负数解^[13]。

证   令，其中0 < u < 1，得。又由文献[11]知Power Lindley分布的ν-分位数为

故可得。

定理2   设r为正整数。PLL(p，θ，α)的r阶原点矩存在，且

证   由（5）式易看出，PLL(p，θ，α)分布的密度函数是与PL(θ，α)分布的分布函数G (x；θ，α)和密度函(x；θ，α)有关的函数，即

又文献[11]证明了PL(θ，α)分布存在有限的各阶矩，故

所以PLL(p，θ，α)的r阶原点矩存在。又因为，故可利用公式(1- ，对进行展开，易得

下面给出PLL(p，θ，α)分布的危险率函数，先回顾相关定义。非负随机变量X的分布函数为F (x)，密度函数为f (x)，其生存函数为S (x)= 1 - F (x)，危险率函数为。

由此得PLL(p，θ，α)分布的生存函数为

(6)

定理3   若随机变量X~PLL(p，θ，α)分布，则其危险率函数为

(7)

且

(8)

(9)

证由（6）式和危险率函数的定义即得（7）式。（8）式显然成立，利用洛必达法则易得（9）式。证毕。

图 2和图 3给出不同参数下PLL(p，θ，α)分布的危险率函数。可见PLL(p，θ，α)分布的危险率函数可以有单调递增型、浴盆型、先单调增再浴盆型、单调递减型四种类型。

定理4   设X₁，X₂，⋯，X_n是来自PLL(p，θ，α)分布的简单随机样本，其分布函数F (·)= F (·；p，θ，α)如（4）式所述，记X₍₁₎ = min { X₁，X₂，⋯，X_n }，则有

其中。

证   记PLL(p，θ，α)分布的密度函数为f (·)= f (·；p，θ，α)，由（5）式及洛必达法则有

再利用文献[14]定理2. 1. 5及定理2. 4. 4即知定理得证。

3 参数的极大似然估计

设X₁，X₂，⋯，X_n是来自PLL(p，θ，α)分布的简单随机样本，样本观测值为x₁，x₂，⋯，x_n，记x_obs={ x₁，x₂，⋯，x_n }，则由（5）式知对数似然函数为

为方便，记ϕ=(p，θ，α)，其极大似然估计为。下面讨论极大似然估计的存在性与相合性。

PLL(ϕ)分布的参数空间记为Φ = { ϕ =(p，θ，α)：p ∈ (0，1)，θ ∈ (0，∞)，α ∈ (0，∞) }，令

则似然方程组可表示为

(10)

对任意的ϕ₀∈ Φ及充分小的ϵ > 0使得球B(ϕ₀，ϵ)={ ϕ：|ϕ - ϕ₀| ≤ ϵ }⊂ Φ，记

定理5    设X₁，X₂，⋯，X_n是来自PLL(ϕ)分布的简单随机样本，则存在ϕ的一个强相合估计，满足

证   由（5）式易得

根据文献[15]定理4.7，为证MLE的存在性与相合性，只需证明对∀ϕ₀ ∈ Φ和充分小的ϵ > 0，有E _ϕ0 [H (X，ϕ₀，ϵ)] < ∞和E_ϕ0 [|lnf (X；ϕ₀)|] < ∞。

事实上，h (x，ϕ)作为Φ上的函数在B (ϕ₀，ϵ)连续，故存在ϕ₁= ϕ₁(x)=(p₁(x), θ₁(x), α₁(x))∈ B (ϕ₀，ϵ), 使得|h (x，ϕ₁)| = H (x，ϕ₀，ϵ)，为方便，分别将p₁ (x)，θ₁ (x)，α₁ (x)记为p₁，θ₁，α₁，所以

显然E_ϕ0[1- G (X；θ₁，α₁)] < ∞。由定理2知E_ϕ0[(1 + θ₁) X ^α₁+ 1] < ∞和, 故只需证明。

因为，且关于x为[0, 1]上的有界函数，故

又因为和关于x分别为[0, 1]和[0, ∞]上的有界函数，且，所以有

从而有E_ϕ0[ H (X，ϕ₀，ϵ)] < ∞。

又因为

证毕。

定理5给出了在n充分大时参数ϕ的极大似然估计的存在性，但并不能保证对任意固定的n，参数ϕ的极大似然估计存在。下面考虑在已知两个参数的情况下，另一个未知参数的极大似然估计的存在性和渐近正态性，将得到类似于文献[3]定理1的结论。由于方法和结论是类似的，我们只考虑当参数θ和α已知时，参数p的极大似然估计。此时关于p的似然方程为

(11)

定理6   设PLL(p，θ，α)分布的参数θ，α已知，p未知。令

则

1）；

2) ∀n≥ 1, 似然方程（11）在上至少有一解(在上补充定义), 称为p的极大似然估计；

3）记f(x；p)=f(x；p，θ，α), 为Fisher信息量，则p的极大似然估计是参数p的渐近正态估计，其渐近分布为N (p，[nI (p)]^-1)。

证   1）由（5）式可知

因为关于p在（0，1）内单调递减，故，由大数定律，以概率1，当n充分大时有，即成立，也就是，因此。

2）由似然方程（11）计算可知

故在Λ_n上方程（10）关于p在（0，1）内至少有一解。

3）仿定理5的证明，将文献[15]定理4. 7用到目前情形，可知是相合的。记参数p的真值为p₀，为证明的渐近正态性，验证文献[16]的定理2. 14的4个条件均被满足：

a) 当x > 0时, 都存在；

b) 对任意充分小的正数ϵ, 当时，计算可知

易证E (Q (X)) < ∞，从而文献[16]的定理2. 14的条件（2）得到验证；

c）令，对∀p > 0，总存在充分小的正数ε₀，满足p - ε₀ > 0。由中值定理，[p - ϵ₀，p + ϵ₀]，使得

因为

易知∫₀^∞ K (x)dx < ∞，所以由控制收敛定理，对∀p有

同理有

d）Fisher信息量为

又Power Lindley-Logarithmic分布为非退化分布，故I (p) > 0。从而文献[16]的定理2. 14的所有条件均满足，于是是p的渐近正态估计。

证毕。

4 数值模拟

由于PLL(p，θ，α)分布的似然方程组表达式复杂不易求出显式解，因此用EM算法（算法介绍详见文献[17]）给出其参数的极大似然估计。

设ξ₁，ξ₂，⋯，ξ_z，⋯是来自PL(θ，α)分布的简单随机样本。Z是与ξ_z(z = 1，2，⋯)相互独立的随机变量，且服从参数为p的Logarithmic分布。令X = min { ξ₁，ξ₂，⋯，ξ_Z}，易得Z关于X的条件期望为

记x_obs=(x₁，x₂，⋯，x_n)为来自PLL(ϕ)分布的观测数据，z =(z₁，z₂，⋯，z_n)为来自数据Z但不能观测的潜在数据，设ϕ^(t)=(p^(t)，θ^(t)，α^(t))为第t + 1次迭代的初始值，则由EM算法可得

(12)

(13)

(14)

其中，且z_i^(t) > 1。因为z_i^(t)已知且大于 1，方程Mplnp + n (1 - p)= 0，当M > n时在(0，1)上有惟一解，故用二分法求解（12），得到p。方程（13）和（14）关于θ和α均在(0，∞)存在解，文献[18]指出了求解非线性方程组的牛顿迭代法具有二阶收敛速度，因此对于（13）和（14）式组成的非线性方程组用牛顿迭代法求解，得到θ和α。

记n为样本容量，则通过以下步骤得到参数ϕ =(p，θ，α)的极大似然估计：

步骤1   产生n个服从U (0，1)分布的随机数，记为u =(u₁，u₂，⋯，u_n)；

步骤2   将u₁，u₂，⋯，u_n分别代入到PLL(ϕ)分布的分位数函数（定理1），计算得到x_i = F^-1 (u_i；ϕ)，i=1，2，⋯，n，则x_obs =(x₁，x₂，⋯，x_n)为来自PLL(ϕ)分布的容量为n的样本；

步骤3   给定参数的初始值ϕ⁽⁰⁾=(p⁽⁰⁾，θ⁽⁰⁾，α⁽⁰⁾)，设定迭代精度为0. 000 1，用EM算法得到参数ϕ的极大似然估计；

步骤4   重复步骤1~3 1 000次，计算得到的1 000个估计值的均值作为参数ϕ的极大似然估计的最终结果。

假设模拟得到的极大似然估计为。记这1000个估计值的均值为

样本均方误差为

样本标准误差为

表 1给出了在样本容量分别为50，100，500时模拟的结果。

从表 1可以看到，用EM算法得到的参数的极大似然估计很好地反映了参数的真值，选取不同的初始值对似然估计值的影响不大，且从整体上来看，随着样本量的增加，估计值越来越高，均方误差和标准误差也越来越小。

下面考虑在两个参数已知情况下，另一参数的极大似然估计及95%渐近置信区间的覆盖率。由于方法类似，我们只考虑当θ和α已知，p未知时的情形。在求解极大似然估计时，同样使用EM算法。

沿用前面的符号，记，得到EM算法求解p^{( t + 1)} 的更新公式

(15)

同方程（12），这里仍用二分法求p^{( t + 1)}。

根据定理6可知固定参数θ和α时，参数p的渐近方差为，因为利用Fisher信息量I (p)计算期望时比较复杂，这里可以用观测信息量代替，记p的渐近方差为var(p)，则有

记参数p的极大似然估计为，其95%渐近置信区间为。

记n为样本容量，则通过以下步骤得到在参数θ，α已知时，未知参数p的极大似然估计：

步骤1   产生n个服从U (0, 1)分布的随机数，记为u =(u₁, u₂, ⋯, u_n)；

步骤2   将u₁，u₂，⋯，u_n分别代入到PLL(ϕ)分布的分位数函数（定理1），计算得到x_i = F^-1 (u_i；ϕ)，i=1，2，⋯，n，则x_obs =(x₁，x₂，⋯，x_n)为来自PLL(ϕ)分布的容量为n的样本；

步骤3   给定参数的初始值p⁽⁰⁾，设定迭代精度为0. 000 1，用EM算法得到参数p的极大似然估计，并判断真值p是否落入95%渐近置信区间；

步骤4   重复步骤1~3 1 000次，计算得到的1 000个估计值的均值作为参数p的极大似然估计的最终结果。

假设模拟得到的极大似然估计为，i= 1，2，⋯，1 000，相对应的渐近置信区间为CI_i= ，记这1 000个估计值的均值为

样本均方误差为

样本标准误差为

参数p的95%置信区间平均长度为

95%置信区间的覆盖率为

通过模拟得到的θ和α已知时，参数p的极大似然估计及渐近置信区间的覆盖率如表 2所示。

从表 2可以看出，对于相同的固定值θ和α，参数p越小，极大似然估计对p的估计效果越好；选取不同的初始值对p的估计影响不大；随着样本容量n的增大，p的95%渐近置信区间越来越精确；p的95%渐近置信区间的实际覆盖率在0. 95左右，说明当θ和α固定时，p的极大似然估计具有良好的渐近正态性。