Geometric分布与寿命分布复合族参数估计的MM算法

熊朝松

(云南师范大学 数学学院,云南 昆明 650500)

0 引言

寿命分布在可靠性分析中扮演着重要角色.在处理生活中记录的寿命数据时,最直接的方法就是用寿命分布进行拟合分析,一个好的拟合结果往往能起到重要的指导作用.然而随着科技的发展和时代的进步,现实生活中观测到的寿命数据越来越多,结构也越来越复杂.经典的寿命分布如Exponential分布和Weibull分布等已经无法完全适应多种多样的寿命数据,在大多数情况下并不能得到一个较好的拟合结果.因此,为了减小在拟合分析中产生的偏差,开发新的寿命分布成为可靠性统计领域的重点研究方向.学者们通过复合一个经典的寿命分布和不含零的离散分布,研究了大量的模拟寿命数据的多参数分布.Adamidis和Loukas[1]首先将Exponential分布和Geometric分布复合,得到了一个新的具有下降风险率的Exponential-Geometric(EG)分布,并详细讨论了包括生存函数、风险函数、似然推断等在内的统计性质.EG分布的提出,为开发新寿命分布提供了一个新的思路.受到Adamidis和Loukas的启发,Barreto-Souza等[2]和顾蓓青等[3]开发了Weibull-Geometric(WG)分布;姚惠等[4]开发了Pareto-Geometric(PG)分布;高艳红和周秀轻[5],Okorie等[6]开发了Rayleigh-Geometric(RG)分布;Alkarni[7]归纳总结了Geometric分布与常用寿命分布的复合形式,将这类分布统称为Geometric分布与寿命分布复合族,并进一步研究了这类分布族的相关统计性质，概率密度函数和分布函数形式,可靠性和失效率函数的显示公式等.在似然推断中,Alkarni[7]建议使用EM算法求解这类复合族分布的极大似然估计.但由于EM算法涉及引入潜变量问题,对于部分初学者来说不容易理解.因此开发一个简单易懂的参数估计方法是有必要的.

与EM算法类似,Minorization-maximization(MM)算法同样也是一种迭代优化算法,其工作原理是通过找到一个替代函数驱动目标函数最大化.值得注意的是,MM算法包含了EM算法,即每个EM算法都是MM算法,反之不然.且MM算法在算法操作上与EM算法有所区别,MM算法不需要像EM算法那样引入潜变量来寻找替代函数,而是直接利用目标函数的极小化函数作为替代函数.由于MM算法具有简单的结构原理和优秀的估计性质,因此被广泛应用于统计学的各个领域[8-10].

本文引入MM算法来计算Geometric分布与寿命分布复合族的极大似然估计,推导了这类分布族的MM算法参数估计式,随后将本文提出的方法应用于几个具体的分布中,最后通过随机模拟实验验证MM算法估计的有效性.

1 MM算法

假设l(θ)是待最大化的目标函数,θ∈Θ表示未知的参数向量,令l(θ)的最大值点为

MM算法的构造方式主要包括两个步骤,具体内容如下.

第一步:Minorization(M)步,主要利用极小化功能找到目标函数的替代函数Q(θ|θ(t)),或称为极小化函数,使其满足

Q(θ|θ(t))≤l(θ),∀θ,θ(t)∈Θ;Q(θ(t)|θ(t))=l(θ(t)),

第二步:Maximization(M)步,通过最大化Q(θ|θ(t))将θ(t)更新为θ(t+1),可以得到

l(θ(t+1))≥Q(θ(t+1)|θ(t))≥Q(θ(t)|θ(t))=l(θ(t)).

MM算法在每次迭代计算时都会增加目标函数值,并且这种上升性质为MM算法提供了显著的数值稳定性.

MM算法的关键步骤就在于第一步中找到合适的替代函数,在此本文介绍一种常用的极小化方法,该方法由Zhou和Lange[11]提出,Tian等[10]对其做出了详细介绍.令r(x)和s(x)是在相同支撑集上的两个密度函数,即

{x|r(x)>0}={x|s(x)>0}.

r(x)和s(x)之间的Kullback-Leibler(KL)距离被定义为

对应的离散形式为

在上述公式中，令n=2,(r1,r2)=(δ0,1-δ0)，(s1,s2)=(δ,1-δ),则有

δ0logδ0+(1-δ0)log(1-δ0)≥
δ0logδ+(1-δ0)log(1-δ),

其中δ,δ0∈(0,1),当且仅当δ=δ0时等号成立.将上述不等式重新排列,即有

(1)

2 Geometric分布与寿命分布复合族的MM算法估计

2.1 Geometric分布与寿命分布复合族

定义1给定一个随机变量Z服从一个零截断Geometric分布,其概率质量函数(pmf)满足

fZ(z|p)=(1-p)pz-1,0

假设T=(T1,…,TZ)是独立同分布,概率密度函数(pdf)为

fTi(x|α)=fT(x|β)|β=(β1,…,βm),m≥1,x,β∈R+

的随机变量,其中对于任意i=1,2,…,随机变量Z和Ti是相互独立的.令随机变量X=min(T1,…,TZ),则根据上述条件,可以得到X|Z=z的pdf为

同时可以获得X和Z的联合概率密度函数为

其中FT(x|β)表示T的累积分布函数(cdf).通过计算容易得到X的边际pdf和cdf分别为

(2)

和

(3)

如果一个随机变量X的pdf和cdf分别为公式(2)和(3)的形式,则称X服从Geometric分布与寿命分布复合族,记作X～GL(p,β).

2.2 极大似然估计的MM算法

l(p,β|Xobs)=l0(p,β|Xobs)+l1(p,β|Xobs),

其中

观察到,l1(p,β|Xobs)含有与-log(1-δ)类似的结构,则考虑通过极小化l1(p,β|Xobs)来获取l(p,β|Xobs)的替代函数.

已知0

0

将δ=p[1-FT(xi|β)]和δ0=p(t)[1-FT(xi|β(t))]代入公式(1),可以得到l1(p,β|Xobs)的极小化函数为

其中c是与参数(p,β)无关的常数项.结合l0(p,β|Xobs)与Q1(p,β|p(t),β(t))可以得到对数似然函数l(p,β|Xobs)的替代函数为

其中

(4)

3 应用

将本文提出的方法应用于四个常见的Geometric分布与寿命分布复合族分布,分别为Exponential-Geometric分布,Weibull-Geometric分布,Pareto-Geometric分布和Rayleigh-Geometric分布.为了方便推导出基于极大似然的MM算法估计公式,首先在此介绍Tian等[12]提出的LOG-BETA函数族和LOG-GAMMA函数族.

定义2假设一个函数在定义域[0,1]上的表达式为

g(θ)=c+alog(θ)+blog(1-θ),θ∈[0,1],

其中,c是与参数θ不相关的常数项,且有a,b≥0.则称此函数属于LOG-BETA函数族,记作g(θ)∈LB(θ).log(θ)和log(1-θ)分别为LOG-BETA函数族的两个基本组装函数.

注意到,当a,b>0时,g(θ)是严格凹函数,则可以得到

(5)

定义3假设一个函数在定义域正实数集R+上的表达式为

g(θ)=c+alog(θ)+b(-θ),θ∈R+,

其中,c是与参数θ不相关的常数项,且有a,b≥0.则称此函数属于LOG-GAMMA函数族,记作g(θ)∈LG(θ).log(θ)和-θ分别为LOG-GAMMA函数族的两个基本组装函数.

同样地,当a,b>0时,g(θ)是严格凹函数,则可以得到

(6)

基于上述的LOG-BETA函数族和LOG-GAMMA函数族的定义与性质,下面详细讨论MM算法在求解几种常见的复合寿命分布的极大似然估计时的应用.

3.1 Exponential-Geometric分布

Adamidis等[1]提出的Exponential-Geometric分布的概率密度函数为

f(x|p,β)=β(1-p)e-βx(1-pe-βx)-2,x>0，

其中00.则对应的对数似然函数为

利用本文2.2节的方法,l(p,β|Xobs)的替代函数为

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

(7)

观察到,MM算法将一个两参数的优化问题转变成了两个单参数优化问题之和,其中关于p的函数属于LOG-BETA函数族,关于β的函数属于LOG-GAMMA函数族.

分别令dQ(p|p(t),β(t))/dp=0和dQ(β|p(t),β(t))/dβ=0,同时利用公式(5)和(6),即可得到参数p和β的MM算法估计式.

3.2 Weibull-Geometric分布

Barreto-Souza等[2]将Weibull分布与Geometric分布复合,得到了三参数的Weibull-Geometric分布,其概率质量函数满足

f(x|p,β,α)=αβα(1-p)xα-1exp[-(βx)α]·{1-pexp[-(βx)α]}-2,x>0，

其中00,α>0.令上述的pmf中的参数β=1/β*,即可得到顾蓓青等[3]考虑的Weibull-Geometric分布.

根据Weibull-Geometric分布的概率密度函数,对应的对数似然函数可写为

利用本文2.2节的方法,l(p,β,α|Xobs)的替代函数为

Q(p,β,α|p(t),β(t),α(t))=Q(p|p(t),β(t),α(t))+Q(β,α|p(t),β(t),α(t))+c,

其中

(8)

容易看出,Q(p|p(t),β(t),α(t))∈LB(p),利用公式(5)可以得到参数p的估计式.参数(β,α)的MM算法迭代公式为

其中α(t+1)为非线性方程

的解.

3.3 Pareto-Geometric分布

姚惠等[4]研究了定义在x>1上的Pareto-Geometric分布.而在现实中观测到的寿命数据大多定义在x>0上,因此根据定义1,可以得到Pareto-Geometric分布的概率密度函数为

f(x|p,β)=β(1-p)(1+x)-(β+1)[1-p(1+x)-β]-2,x>0,

其中00.故对应的对数似然函数为

利用本文2.2节的方法,其替代函数可以写为

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

令dQ(p|p(t),β(t))/dp=0,dQ(β|p(t),β(t))/dβ=0,同时利用公式(5)和(6),即可得到Pareto-Geometric分布参数p和β的MM算法估计式.

3.4 Rayleigh-Geometric分布

高艳红等[5]提出的Rayleigh-Geometric分布的概率密度函数为

f(x|p,β)=2β(1-p)xe-βx2(1-pe-βx2)-2,x>0,

其中00.对应的对数似然函数为

l(p,β|Xobs)替代函数为

Q(p,β|p(t),β(t))=Q(p|p(t),β(t))+Q(β|p(t),β(t))+c,

其中

(10)

令dQ(p|p(t),β(t))/dp=0和dQ(β|p(t),β(t))/dβ=0,同时利用公式(5)和(6),即可得到Rayleigh-Geometric分布参数p和β的MM算法估计式.

4 随机模拟

本节基于第3节中提到的例子进行了数值实验,以此评估所提出的MM算法的实际性能.四组实验分别生成样本量为100和500的EG分布、WG分布、PG分布和RG分布独立随机样本,并在不同的参数值组合下,利用本文提出的MM算法计算参数的极大似然估计值.将每组试验重复进行1000次,计算出参数的极大似然估计平均值和估计的均方误差.

第一组实验中,分别设置EG分布的参数(p,β)真值为(0.1,0.1),(0.9,0.1),(0.1,1)和(0.9,1).利用3.1节定义的MM算法来计算(p,β)的极大似然估计,随机模拟结果如表1所列.

表1 EG分布参数的极大似然估计

第二组实验中,分别设置WG分布的参数(p,β,α)真值为(0.2,0.2,2)和(0.8,2,3).利用3.2节所定义的MM算法来计算(p,β,α)的极大似然估计,随机模拟结果如表2所列.

表2 WG分布参数的极大似然估计

第三组实验中,分别设置PG分布的参数(p,β)真值为(0.3,0.3),(0.7,0.3),(0.3,3)和(0.7,3).利用3.3节定义的MM算法来计算(p,β)的极大似然估计,随机模拟结果如表3所列.

表3 PG分布参数的极大似然估计

第四组实验中,分别设置RG分布的参数(p,β)真值为(0.4,0.4),(0.6,0.4),(0.4,4)和(0.6,4).利用3.4节定义的MM算法来计算(p,β)的极大似然估计,随机模拟结果如表4所列.

表4 RG分布参数的极大似然估计

从表1至表4的结果可以看出,针对于4个不同的分布,在不同的样本量和参数真值的组合下,各个参数的极大似然估计值均能收敛到参数真值,同时随着样本量的增加,估计的偏差和均方误差越来越小,表明参数估计的效果和稳定性越来越好.

5 结束语

本文首先给出了计算Geometric分布与寿命分布复合族的参数极大似然估计的MM算法,为这类寿命分布族的参数估计提供了一个新的方案.与现有的EM算法不同的是,MM算法是完全基于观测数据的,不再需要引入潜变量.此外,MM算法将多维优化问题转变成了低维优化问题之和,极大地减少了计算难度和成本.随后本文将提出的方法应用于EG分布、WG分布、PG分布和RG分布的极大似然估计中,得到了各个分布参数的MM算法估计式.随机模拟结果表明MM算法估计结果准确有效,并且随着样本量的增加,估计的效果也越来越好.本文提出的方法为Geometric分布与寿命分布复合族的参数估计提供了一个统一框架,能对后续开发的新寿命分布的参数估计有所帮助.