具有单调失效率的寿命分布可以用来模拟可靠性、生存分析、保险精算等许多领域的寿命数据, 关于这些分布的统计推断引起了很多统计学者的兴趣.

随着寿命分布理论逐步发展, 人们提出了各种新型的寿命分布.设$\xi_1, \cdots, \xi_n, \cdots$为独立同分布随机变量序列, 它们的共同分布函数为已知的寿命分布函数$G(x)$; 另设$N$是与$\{\xi_n\}$相互独立的随机变量, 其分布律为$p_k=P(N=k), k=1, 2, \cdots$.令$X=\min\{\xi_{i}:1\leq i\leq n\}$, 称$X$为由$N$与$\{\xi_n\}$“混合”而成的新型寿命分布. Adamidis和Loukas ^[1]设$G(x)$为指数分布, $N$服从几何分布, 混合得到的新分布称之为Exponential-Geometric分布.随后Kus ^[2]假定$G(x)$服从指数分布, $N$服从泊松分布, 混合成新分布Exponential-Poisson分布. Tahmasbi和Rezaei ^[3]取$G(x)$为指数分布, $N$为logarithmic分布, 所得到的新分布称为Exponential-Logarithmic分布. Chahkandi和Ganjali ^[4]假定$G(x)$仍为指数分布, $N$为二项分布, 构成的新分布命名为Exponential-Binomial分布. Barreto-Souza等^[5]取$G(x)$为Weibull分布, $N$服从几何分布, 将混合得到的新分布称为Weibull-Geometric分布.姚惠等(2011)^[6]令$G(x)$服从Pareto分布, $N$服从Geometric分布, 生成新分布Pareto-Geometric分布.这些文献提出了一些具有单调失效率的寿命分布, 得到了这些分布的各种性质, 极大似然估计的存在唯一性, 并进行了数值模拟.

Pareto分布常常用来描述各种社会经济、物理以及生物现象, 并且在军事领域、天文领域也有应用. Soliman ^[7]提及了Pareto分布作为寿命分布模型的重要性. Johnson等^[8]详细研究了Logarithmic作为计数分布的重要性以及它的各种性质.本文将上述$G(x)$取为Pareto分布, $N$取为Logarithmic分布, 混合成一个新的寿命分布, 称为Pareto-Logarithmic分布, 研究了该分布的矩、熵、失效率函数、平均剩余寿命和参数极大似然估计, 应用EM算法求参数的极大似然估计, 并进行了蒙特卡罗模拟.

设$\xi_{1}, \xi_{2}, \cdots, \xi_{n}, \cdots$为独立同分布的随机变量序列, 其共同的分布是形状参数为$\beta>0$的Pareto分布, 即分布函数为$G(x;\beta)=[1-(x+1)^{-\beta}]I_{(0, +\infty)}(x), $又设N是与$\{\xi_{i}, i=1, 2, \cdots \}$相互独立的随机变量, 且服从参数为$p\in(0, 1)$的Logarithmic分布, 其分布律为

$P_{N}(n;p)=\frac{(1-p)^{n}}{-n\ln(p)}, n=1, 2, \cdots.$

令$X=\min\{\xi_{1}, \xi_{2}\cdots, \xi_{N}\}$, 则$X$的分布函数为

$\begin{equation} \label{eq:2.1} F_{X}(x;p, \beta)=\sum\limits_{n = 1}^{ + \infty } {{P_N}} (n;p)\{1-[1-G(x)]^{n}\}=\Big\{1-\frac{\ln [1-(1-p)(x+1)^{-\beta}]}{\ln p}\Big\}I_{(0, +\infty)}(x), \end{equation}$

(2.1)

其概率密度函数为

$ \begin{equation} \label{eq:2.2} f_{X}(x;p, \beta)= \frac{\beta}{-\ln p}\frac{(1-p)(x+1)^{-(\beta+1)}I_{(0, +\infty)}(x)}{1-(1-p)(x+1)^{-\beta}}. \end{equation}$

(2.2)

仿文献[1-6], 给出如下定义:

定义2.1  称(2.1) 或(2.2) 式给出的分布是参数为$p, \beta$的Pareto-Loagrithmic分布(简称PL分布), 记为$PL(p, \beta)$, 其中参数$p\in(0, 1), \beta\in(0, +\infty)$.

由(2.2) 式可得, $\forall p, \beta$, $f_{X}(x;p, \beta)$关于$x$都是单调递减的. 图 1给出$f_{X}(x;p, \beta)$在参数取某些特定值时的图像.

按照文[9]的定义, 称分布函数$F(x)$为厚尾分布, 如果它的尾函数$S(x)=1-F(x)$以幂函数形式衰减, 即满足$1-F(x)=x^{-\frac{1}{r}}l_{F}(x), r>0, $这里$l_{F}$是无穷远处的缓变函数, 即对所有的$\eta>0, \lim\limits_{x\rightarrow\infty}\frac{l_{F}(\eta x)}{l_{F}(x)}=1.$

由(2.1) 式, 易见$PL(p, \beta)$是厚尾分布.

由文[3]的定义, 密度函数为

$g_{X}(x;p, \beta)=\frac{\beta}{-\ln p}\frac{(1-p)e^{-\beta x}I_{(0, +\infty)}(x)}{1-(1-p)e^{-\beta x}}$

的分布称为Exponential-Logarithmic分布, 记为$EL(p, \beta)$分布, 其中参数$p\in(0, 1), \beta\in(0, +\infty)$.

如果$X\thicksim PL(p, \beta)$, 则由(2.2) 式易得$Y:=\ln(X+1)\thicksim EL(p, \beta).$

设$0<r<1$, 由(2.1) 式可得$PL(p, \beta)$的$r$分位数为$(\frac{1-p}{1-p^{1-r}})^{\frac{1}{\beta}}-1$.特别地, 中位数为$(1 +p^{\frac{1}{2}})^{\frac{1}{\beta}}-1.$

记$\Phi(z, s, v):=\sum\limits_{n=0}^{\infty}(v+n)^{-s}z^{n}$为高等超越函数(详见文[10]).

定理3.1 设$X\sim PL(p, \beta)$, 则有

(1) 设$k$为正整数, 当$\beta>k$时, $PL(p, \beta)$分布的$k$阶矩为

$E(X^{k})=-\frac{(1-p)}{\ln p}\sum\limits_{j = 0}^k {C_k^j} (-1)^{k-j}\Phi(1-p, 1, 1-\frac{j}{\beta}).$

特别地, 当$\beta>1$时,

$E(X)=-\frac{(1-p)}{\ln p} \Phi(1-p, 1, 1-\frac{1}{\beta})-1.$

(2) 当$\beta>2$时, $PL(p, \beta)$分布的方差为

${\rm Var}(X)=-\frac{(1-p)}{\ln p} \Phi(1-p, 1, 1-\frac{2}{\beta})-\frac{(1-p)^{2}}{(\ln p)^{2}} \Phi^{2}(1-p, 1, 1-\frac{1}{\beta}).$

证   (1) 当$\beta>k$时, 利用(2.2) 式和泰勒展开, 有

$\begin{aligned} E(X^{k}) &= \int_{0}^{\infty} x^{k}\frac{\beta}{-\ln p} \frac{(1-p)(x+1)^{-(\beta+1)}}{1-(1-p)(x+1)^{-\beta}} dx \\ &= \int_{1}^{\infty} (x-1)^{k}\frac{\beta}{-\ln p}(1-p)x^{-(\beta+1)} \sum\limits_{i=0}^{+\infty}(1-p)^{i}x^{-\beta i} dx \\ &= -\frac{(1-p)}{\ln p}\sum\limits_{j=0}^{k}C_{k}^{j}(-1)^{k-j}\Phi(1-p, 1, 1-\frac{j}{\beta}). \end{aligned}$

(2) 当$\beta>2$时,

$\begin{aligned} {\rm Var}(X) &= E(X^2)-E^{2}(X)\\ &= -\frac{(1-p)}{\ln p}\sum\limits_{j=0}^{2}C_{2}^{j}(-1)^{2-j}\Phi(1-p, 1, 1-\frac{j}{\beta})\\ &-(-\frac{(1-p)}{\ln p} \Phi(1-p, 1, 1-\frac{1}{\beta})-1)^2 \\ &= -\frac{(1-p)}{\ln p} \Phi(1-p, 1, 1-\frac{2}{\beta})-\frac{(1-p)^{2}}{(\ln p)^{2}} \Phi^{2}(1-p, 1, 1-\frac{1}{\beta}). \end{aligned}$

熵是度量随机变量不确定性的工具.熵有很多种, 其中Renyi熵^[11]和Shannon熵^[12]比较常用, 我们先回顾相关定义.

定义3.1  对于一个密度函数为$f(\cdot)$的非负随机变量$X$,

(1) 其Renyi熵为$I_{R}(r)=\displaystyle\frac{1}{1-r}\ln\int_{0}^{+\infty}f(x)^{r}dx$, 这里$r>0, r\neq1$;

(2) 其Shannon熵为$E[-\ln f(X)]$.

下面给出$PL(p, \beta)$分布的Renyi熵和Shannon熵,

$\begin{aligned} I_{R}(r) &= \frac{1}{1-r}\ln\int_{0}^{+\infty}f^{r}(x;p, \beta)dx= \frac{1}{1-r}\ln\int_{0}^{+\infty} \\ &\frac{\beta^{r}}{(-1)^{r}\ln^{r} p}\frac{(1-p)^{r}(x+1)^{-r(\beta+1)}}{(1-(1-p)(x+1)^{-\beta})^r}dx \\ &= \frac{1}{1-r}\ln\int_{0}^{+\infty} \frac{\beta^{r}(1-p)^{r}}{(-1)^{r}\ln^{r} p}(x+1)^{-r(\beta+1)}\\ &\big[1+\sum\limits_{j=1}^{+\infty}\frac{r\cdots(r+j-1)}{j!}(1-p)^{j}(x+1)^{-\beta j}\big]dx \\ &= \frac{1}{1-r}\ln \frac{\beta^{r}(1-p)^{r}}{(-1)^{r}\ln^{r} p}[\frac{1}{r(\beta+1)-1} \\ &+\sum\limits_{j=1}^{+\infty}\frac{r\cdots(r+j-1)}{j!}(1-p)^{j}\frac{1}{(j+r)\beta+r-1}], \end{aligned}$

这里要求$r>\frac{1}{\beta+1}, r\neq1$,

$\begin{aligned} E[-\ln f(X)] =& \ln(-\ln p)-\ln\beta-\ln(1-p)-\int_{0}^{+\infty}\ln(x+1)^{-(\beta+1)}f(x;p, \beta)dx\\ &+\int_{1}^{+\infty}\ln(1-(1-p)x^{-\beta})d[1-\frac{\ln(1-(1-p)x^{-\beta})}{\ln p}] \\ =& \ln(-\ln p)-\ln\beta-\ln(1-p)+\frac{\ln p}{2}+\frac{\beta+1}{-\ln p}\sum\limits_{j=0}^{+\infty}(1-p)^{j+1}\frac{1}{\beta(j+1)^{2}}. \end{aligned}$

下面研究PL分布的失效率函数和平均剩余寿命, 先回顾相关定义.对分布密度函数为$f(x)$, 分布函数为$F(x)$的非负随机变量$X$, 其生存函数为$S(x):=1-F(x)$, 失效率函数为$h(x):=\frac{f(x)}{S(x)}$, 平均剩余寿命为$m(x):=E(X-x|X\geq x).$

寿命分布类及其基本性质是可靠性理论的一个重要研究分支, 详细讨论可参看文献[13-14].近几年, 有关年龄性质的研究很活跃, 一些新的年龄性质被陆续引入, 年龄性质的可靠性特征也渐渐被人们发掘.我们先回顾相关年龄性质的定义.

定义3.2  设$X$是非负随机变量, 其分布函数、生存函数以及失效率函数分别为$F(t)$、$\overline{F}(t)$以及$h(t)$.称$X$具有

(1) DFR(decreasing failure rate)性质, 如果$h(t)$是单调递减函数;

(2) DFRA(decreasing failure rate average)性质, 如果$\displaystyle\frac{\int_{0}^{t}h(u) du}{t}$, $\forall t>0$单调递减;

(3) NWU(new worse than used)性质, 如果$\forall t, x>0$, 都有$\overline{F}(t+x)\geq \overline{F}(t)\overline{F}(x)$;

(4) NWUFR(new worse than used in failure rate)性质, 如果$\forall t>0$, 都有$h(t)\leq h(0)$;

(5) NWAFR(new worse than used average in failure rate)性质, 如果$\forall t>0$, 都有

$\displaystyle\frac{\int_{0}^{t}h(u) du}{t}\leq h(0).$

有关上述年龄性质的详细讨论见Nanda等^[15].

定理3.2  $PL(p, \beta)$分布具有如下特性:

(1) 失效率函数为

$h(x;p, \beta)=\frac{-\beta(1-p)(x+1)^{-(\beta+1)}I_{(0, +\infty)}(x)}{[1-(1-p)(x+1)^{-\beta}]\ln [1-(1-p)(x+1)^{-\beta}]}, $

对所有的$p, \beta$, $h(x;p, \beta)$关于$x$都单调递减, 即PL分布具有DFR性质, 且

$\mathop {\lim }\limits_{x \to {0^ + }} h(x;p, \beta)=\frac{-\beta(1-p)}{p\ln p}, \lim\limits_{x\rightarrow+\infty}h(x;p, \beta)=0.$

(2) 当$\beta>1$时, 平均剩余寿命为

$m(x;p, \beta)=\frac{-(x+1)^{-\beta+1}(1-p)\Phi((1-p)(x+1)^{-\beta}, 1, 1-\frac{1}{\beta})}{\ln(1-(1-p)(x+1)^{-\beta})}-(x+1), $

$m(x;p, \beta)$关于$x$单调递增, 且

$\lim\limits_{x\rightarrow0^{+}}m(x;p, \beta)=E(x), \lim\limits_{x\rightarrow+\infty}m(x;p, \beta)=+\infty.$

证   (1) 由(2.1) 式易见$PL(p, \beta)$分布的生存函数为

$S(x;p, \beta)=1-F(x;p, \beta)=\{\frac{\ln [1-(1-p)(x+1)^{-\beta}]}{\ln p}\}I_{(0, +\infty)}(x), $

从而可得其失效率函数为

$h(x;p, \beta)=\frac{-\beta(1-p)(x+1)^{-(\beta+1)}I_{(0, +\infty)}(x)}{[1-(1-p)(x+1)^{-\beta}]\ln [1-(1-p)(x+1)^{-\beta}]}, $

其中参数$p\in(0, 1), \beta\in(0, +\infty)$, 从而求得

$\begin{eqnarray*} &&\lim\limits_{x\rightarrow0^{+}}h(x;p, \beta)=\lim\limits_{x\rightarrow0^{+}}\frac{-\beta(1-p)(x+1)^{-(\beta+1)}}{[1-(1-p)(x+1)^{-\beta}]\ln [1-(1-p)(x+1)^{-\beta}]}=\frac{-\beta(1-p)}{p\ln p}, \\ &&\lim\limits_{x\rightarrow+\infty}h(x;p, \beta)=\lim\limits_{x\rightarrow+\infty}\frac{-\beta(1-p)(x+1)^{-(\beta+1)}}{[1-(1-p)(x+1)^{-\beta}]\ln [1-(1-p)(x+1)^{-\beta}]}=0. \end{eqnarray*}$

令

$\eta(x)=-\frac{f'(x)}{f(x)}=\frac{-(1-p)(x+1)^{-(\beta+1)}+(\beta+1)(x+1)^{-1}}{1-(1-p)(x+1)^{-\beta}}.$

记

$q(x)=-(1-p)(x+1)^{-(\beta+1)}+(\beta+1)(x+1)^{-1}, $

由于$q(x)$是单调递减的, 故$\eta(x)$单调递减, 根据文[16] Glaser定理可得$h(x;p, \beta)$是单调递减的, 即具有DFR性质.

(2) 由平均剩余寿命的定义以及(2.2) 式可得, 当$x>0, \beta>1$时, $PL(p, \beta)$分布的平均剩余寿命为

$\begin{aligned} m(x;p, \beta) &= \frac{E[(X-x)I_{(X>x)}]}{p(X>x)}\\ &=\frac{-(x+1)^{-\beta+1}(1-p)\Phi((1-p)(x+1)^{-\beta}, 1, 1-\frac{1}{\beta})}{\ln(1-(1-p)(x+1)^{-\beta})}-(x+1). \\ \end{aligned}$

由$h(x;p, \beta)$单调递减可知$m(x;p, \beta)$关于$x$单调递增, 且

$\begin{eqnarray*} \lim\limits_{x\rightarrow0^{+}}m(x;p, \beta)&=&-\frac{(1-p)}{\ln p} \Phi(1-p, 1, 1-\frac{1}{\beta})-1=E(x), \\ \lim\limits_{x\rightarrow+\infty}m(x;p, \beta)&=& \lim\limits_{x\rightarrow+\infty}\frac{-(1-p)(x+1)^{-\beta+1}\sum\limits_{j=0}^{+\infty}[\frac{1}{1+j-\frac{1}{\beta}}-\frac{1}{1+j}](1-p)^{j}(x+1)^{-\beta j}}{\ln(1-(1-p)(x+1)^{-\beta})}\\ &=& +\infty, \end{eqnarray*}$

这里因为

$\lim\limits_{x\rightarrow+\infty}\frac{-(1-p)(x+1)^{-\beta}}{\ln(1-(1-p)(x+1)^{-\beta})}=1.$

注  由于$PL(p, \beta)$分布具有DFR性质, 由文献[15]可以得到$PL(p, \beta)$分布具有DFRA、NWU、NWUFR、NWAFR性质.

图 2给出$h_{X}(x;p, \beta)$在参数取某些特定的值时的图像.

设$X_{1}, X_{2}, \cdots, X_{n}$为取自于$PL(p, \beta)$分布的简单随机样本, 样本观察值为$\{x_{i}, i=1, 2, \cdots, n\}$, 记$x_{obs}=(x_i, i=1, 2, \cdots, n)$, 则由(2.2) 式可得其似然函数为

$L(p, \beta;x_{obs})=\prod\limits_{i = 1}^n {\frac{\beta }{{ - \ln p}}} \frac{(1-p)(x_{i}+1)^{-(\beta+1)}I_{(0, +\infty)}(x_{i})}{1-(1-p)(x_{i}+1)^{-\beta}};$

对数似然函数为

$\begin{eqnarray*}l(p, \beta;x_{obs})&=&n\ln\beta-n\ln(-\ln p)+n\ln(1-p)\\ &&-(\beta+1)\sum\limits_{i=1}^{n}\ln(x_{i}+1)-\sum\limits_{i=1}^{n}\ln(1-(1-p)(x_{i}+1)^{-\beta});\end{eqnarray*}$

似然方程组为

$ \begin{eqnarray} \label{eq:4.1} &&\frac{\partial l}{\partial p}=-\frac{n}{p\ln p}-\frac{n}{1-p}-\sum\limits_{i=1}^{n}\frac{(x_{i}+1)^{-\beta}}{1-(1-p)(x_{i}+1)^{-\beta}}=0, \end{eqnarray}$

(4.1)

$ \begin{eqnarray} \label{eq:4.2} &&\frac{\partial l}{\partial \beta}=\frac{n}{\beta}-\sum\limits_{i=1}^{n}\ln(x_{i}+1)-\sum\limits_{i=1}^{n}\frac{(1-p)(x_{i}+1)^{-\beta}\ln(x_{i}+1)}{1-(1-p)(x_{i}+1)^{-\beta}}=0. \end{eqnarray}$

(4.2)

似然方程的解就是参数的极大似然估计, 下面讨论$PL(p, \beta)$分布两参数$p, \beta$在一个参数确定时, 另一个参数极大似然估计的存在唯一性, 得到了与文献[3]定理1和定理2类似的结论.

定理4.1  设$PL(p, \beta)$分布中参数$p$已知, 则关于$\beta$的方程(4.2) 在有唯一解$\widehat{\beta}$, 且解在$(pt^{-1}, t^{-1})$中, 这里$t=\frac{1}{n}\sum\limits_{i=1}^{n}\ln(x_{i}+1)$, 且$\widehat{\beta}$是参数$\beta$的相合估计.

证  记

$\begin{eqnarray*} &&g(\beta;p)=\frac{\partial l}{\partial \beta}=\frac{n}{\beta}-\sum\limits_{i=1}^{n}\ln(x_{i}+1)-\sum\limits_{i=1}^{n}\frac{(1-p)(x_{i}+1)^{-\beta}\ln(x_{i}+1)}{1-(1-p)(x_{i}+1)^{-\beta}}, \\ &&\omega(\beta)=\sum\limits_{i=1}^{n}\frac{(1-p)(x_{i}+1)^{-\beta}\ln(x_{i}+1)}{1-(1-p)(x_{i}+1)^{-\beta}}, \end{eqnarray*}$

因为$\omega(\beta)$关于$\beta$是单调递减的且$\omega(\beta)>0$, 故有

$g(\beta;p)<\frac{n}{\beta}-\sum\limits_{i=1}^{n}\ln(x_{i}+1), $

当$\beta>t^{-1}$时, 有$g(\beta;p)<0$, 另$\lim\limits_{\beta\rightarrow0}\omega(\beta)=\frac{1-p}{p}\sum\limits_{i=1}^{n}\ln(x_{i}+1)$, 故

$g(\beta;p)>\frac{n}{\beta}-\sum\limits_{i = 1}^n {\ln } ({x_i} + 1)-\frac{1-p}{p}\sum\limits_{i=1}^{n}\ln(x_{i}+1)=\frac{n}{\beta}-\frac{n}{p}t, $

故当$\beta<pt^{-1}$时, 有$g(\beta;p)>0$, 所以方程(4.2) 在$(pt^{-1}, t^{-1})$区域内至少有一解, 在其它区域无解.又因为

$g'(\beta;p)=-\frac{n}{\beta^{2}}+\mathop \sum \limits_{i = 1}^n \frac{(1-p)(x_{i}+1)^{-\beta}\ln^{2}(x_{i}+1)}{[1-(1-p)(x_{i}+1)^{-\beta}]^{2}}, $

如果$g'(\beta;p)$对所有的$\beta>0$都是单调的, 因为$\lim\limits_{\beta\rightarrow+\infty}g'(\beta;p)=0$, 故$g(\beta;p)$关于$\beta$单调, 解的唯一性得证.如果$g'(\beta;p)$并不是对所有的$\beta>0$都是单调的, 设$\beta_{0}$为$g'(\beta;p)=0$的根, 则

$\frac{n}{\beta_{0}}=\beta_{0}\mathop \sum \limits_{i = 1}^n \frac{(1-p)(x_{i}+1)^{-\beta_{0}}\ln^{2}(x_{i}+1)}{[1-(1-p)(x_{i}+1)^{-\beta_{0}}]^{2}}.$

得到

$g(\beta_{0};p)=\mathop \sum \limits_{i = 1}^n \ln(x_{i}+1)\frac{(1-p)\beta_{0}(x_{i}+1)^{-\beta_{0}}\ln(x_{i}+1)-1+(1-p)(x_{i}+1)^{-\beta_{0}}}{[1-(1-p)(x_{i}+1)^{-\beta_{0}}]^{2}}.$

令

$m(\beta_{0})=(1-p)\beta_{0}(x_{i}+1)^{-\beta_{0}}\ln(x_{i}+1)-1+(1-p)(x_{i}+1)^{-\beta_{0}}, $

因为$m(\beta_{0})$关于$p$单调递减, 故

$m(\beta_{0})<\beta_{0}(x_{i}+1)^{-\beta_{0}}\ln(x_{i}+1)+(x_{i}+1)^{-\beta_{0}}-1<0, $

从而对所有的$\beta_{0}$都有$g(\beta_{0};p)<0$, 即$g(\beta;p)$的稳定点都有$g(\beta;p)<0$, 又因为$\lim\limits_{\beta\rightarrow+\infty}g(\beta;p)=-nt<0$, 故(4.2) 式的解唯一.

容易验证文[17]中定理2.13的条件, 故极大似然估计$\widehat{\beta}$关于$\beta$是相合的.证毕.

注  当参数$p$已知时, 在实际工作中, 如果由历史数据知$\beta$不会很小, 可假定$\beta$的参数空间为$\Theta_{\beta}=(\varepsilon, \infty)$, 其中$\varepsilon>0$为很小的常数, 则经计算得

(1)

$\int_{0}^{+\infty}\frac{\partial f(x)}{\partial\beta}dx=\int_{0}^{+\infty}\frac{\partial^{2} f(x)}{\partial\beta^{2}}dx=0;$

(2)

$\begin{eqnarray*}&&\frac{\partial\ln f(x)}{\partial\beta}=\frac{1}{\beta}-\ln(x+1)-\frac{(1-p)(x+1)^{-\beta}\ln(x+1)}{1-(1-p)(x+1)^{-\beta}}, \\ &&E\big[\frac{1}{\beta}-\ln(X+1)-\frac{(1-p)(X+1)^{-\beta}\ln(X+1)}{1-(1-p)(X+1)^{-\beta}}\big]^{2}<+\infty;\end{eqnarray*}$

(3)

$\left|\frac{\partial^{3}\ln f(x)}{\partial\beta^{3}}\right|<\frac{2}{\varepsilon^{3}}+\frac{2[\ln(x+1)]^{3}}{p^{3}}=H(x), $

且$E[H(X)]<+\infty$.故根据文[17]定理2.14, 可知$\widehat{\beta}$是关于$\beta$的渐近正态估计.

定理4.2  设$PL(p, \beta)$分布中参数$\beta$已知, 则关于$p$的方程(4.1) 在

$\Lambda_{n}:=\{n<2\sum\limits_{i=1}^{n}(x_{i}+1)^{-\beta}\}$

上至少有一解, 且$P(\Lambda_{n})\rightarrow1(n\rightarrow\infty)$.

证  根据式(4.1) 可得

$\lim\limits_{p\rightarrow0}\frac{\partial l}{\partial p}=+\infty, \lim\limits_{p\rightarrow1}\frac{\partial l}{\partial p}=\frac{n}{2}-\sum\limits_{i=1}^{n}(x_{i}+1)^{-\beta}.$

故在$\Lambda_{n}$上方程(4.1) 在$p\in(0, 1)$内至少有一解, 存在性得证.根据(2.2) 式可求得

$E((X_{i}+1)^{-\beta})=\int_{0}^{\infty}(x+1)^{-\beta}f_{X}(x;p, \beta)dx=\frac{1}{1-p}+\frac{1}{\ln p}, $

显然$\frac{1}{1-p}+\frac{1}{\ln p}$关于$p$在$(0, 1)$内单调递减且$\lim\limits_{p\rightarrow1}\frac{1}{1-p}+\frac{1}{\ln p}=\frac{1}{2}, $故有$E((X_{i}+1)^{-\beta})>\frac{1}{2}$, 根据大数定律有$P(\Lambda_{n})\rightarrow1(n\rightarrow\infty)$.证毕.

由于极大似然方程组的表达式比较复杂, 不容易求解, 而EM算法是一种具有良好性质的稳定简单的迭代算法(详见文[18]).下面我们给出$PL(p, \beta)$分布的极大似然估计的EM算法.

设$\xi_{1}, \xi_{2}, \cdots, \xi_{s}\cdots$为取自于参数为$\beta$的Pareto分布的简单随机样本. Z是与$\xi_{s}(s=1, 2, \cdots)$独立的随机变量, 且服从参数为$p$的Logarithmic分布, 令$X=\min\{\xi_{1}, \xi_{2}\cdots, \xi_{Z}\}$, 则$X$和$Z$的联合概率密度函数为

$f(x, z;p, \beta)=\frac{(1-p)^{z}}{-\ln(p)}\beta(x+1)^{-(\beta z+1)}I_{(0, +\infty)}(x), $

这里$z=1, 2, \cdots$, 参数$p\in(0, 1), \beta\in(0, +\infty)$. $Z$关于$X$的条件概率密度为

$\begin{eqnarray*}P_{Z|X}(z|x;p, \beta)&=&P(Z=z|X=x)\\ &=&(1-p)^{z-1}(x+1)^{-\beta(z-1)}[1-(1-p)(x+1)^{-\beta}]I_{(0, +\infty)}(x), \end{eqnarray*}$

这里$z=1, 2, \cdots$, 参数$p\in(0, 1), \beta\in(0, +\infty)$.则$Z$关于$X$的条件期望为

$E(Z|X;p, \beta)=[1-(1-p)(x+1)^{-\beta}]^{-1}.$

当初值选取$(p^{(0)}, \beta^{(0)})$时, 不妨设设$ p^{(t)}, \beta^{(t)}$为第$t+1$次迭代开始时的估计值.其EM算法的迭代公式为

$\begin{eqnarray} \label{eq:3} &&\beta^{(t+1)}=n\big[\sum\limits_{i=1}^{n}\frac{\ln(x_{i}+1)}{1-(1-p^{(t)})(x_{i}+1)^{-\beta^{(t)}}}\big]^{-1};\end{eqnarray}$

(4.3)

$\begin{eqnarray} \label{eq:4} &&p^{(t+1)}=\frac{-n(1-p^{(t+1)})}{\ln(p^{(t+1)})\sum\limits_{i=1}^{n}[1-(1-p^{(t)})(x_{i}+1)^{-\beta^{(t)}}]^{-1}}. \end{eqnarray}$

(4.4)

因为$Mx\ln x+n(1-x)=0(M>n)$在$x\in(0, 1)$里只有一个解, 故可以采用二分法迭代求出$p^{(t+1)}$.

按照通常的大样本理论, 参数$\theta=(p, \beta)^{'}$的极大似然估计$(\widehat{p}, \widehat{\beta})^{'}$近似地服从二元正态分布$N(\theta, J(\theta)^{-1})$, 其中$J(\theta)=E(I;\theta)=\left( \begin{array} {lcr} J_{11}&J_{12} \\ J_{21}&J_{22} \end{array} \right) $, 而$I=I(\theta, x_{obs})=\left( \begin{array} {lcr} I_{11}&I_{12} \\ I_{21}&I_{22} \end{array} \right)$为观测信息矩阵, 由式(4.1) 和(4.2) 求得

$\begin{eqnarray*} &&I_{11}=-\frac{\partial^{2}l}{\partial^{2}p}=-\frac{n(\ln p+1)}{p^2\ln^2 p}+\frac{n}{(1-p)^2}-\sum\limits_{i=1}^{n}\frac{(x_{i}+1)^{-2\beta}}{[1-(1-p)(x_{i}+1)^{-\beta}]^{2}}, \\ &&I_{12}=I_{21}=-\frac{\partial^{2}l}{\partial p\partial\beta}=-\sum\limits_{i=1}^{n}\frac{\ln(x_{i}+1)(x_{i}+1)^{-\beta}}{[1-(1-p)(x_{i}+1)^{-\beta}]^{2}}, \\ &&I_{22}=-\frac{\partial^{2}l}{\partial^{2}\beta}=\frac{n}{\beta^2}-(1-p)\sum\limits_{i=1}^{n}\frac{\ln^2(x_{i}+1)(x_{i}+1)^{-\beta}}{[1-(1-p)(x_{i}+1)^{-\beta}]^{2}}.\end{eqnarray*}$

记

$g(a)=\int_{1}^{a}\frac{{\rm ln} x}{1-x}dx, $

应用如下积分公式

$\begin{eqnarray*} &&\int\frac{2u\ln u}{(1-u)^3}du=\frac{u^{2}\ln u}{(1-u)^2}-\frac{1}{1-u}-\ln(1-u).\\ &&\int\frac{u(\ln u)^2}{(1-u)^3}du=\frac{u^{2}\ln^{2}u}{2(1-u)^2}-\frac{u\ln u}{(1-u)}-\ln(1-u)[\ln u+1]+\int\frac{\ln(1-u)}{u}du.\end{eqnarray*}$

得到

$\begin{eqnarray*}&&E\big[\frac{(X+1)^{-2\beta}}{[1-(1-p)(X+1)^{-\beta}]^{2}}\big]=-\frac{1}{(1-p)^2\ln p}(\frac{3}{2}+\frac{1}{2p^2}-\frac{2}{p}-\ln p);\\ &&E\big[\frac{\ln(X+1)(X+1)^{-\beta}}{[1-(1-p)(X+1)^{-\beta}]^{2}}\big]=-\frac{1-p+p\ln p}{2\beta p(1-p)\ln p};\\ &&E\big[\frac{\ln^2(X+1)(X+1)^{-\beta}}{[1-(1-p)(X+1)^{-\beta}]^{2}}\big]=\frac{1}{\beta^{2}(1-p)}+\frac{g(p)}{\beta^{2}(1-p)\ln p}.\end{eqnarray*}$

从而求得

$\begin{eqnarray*}&&J_{11}=-\frac{n(\ln p+1)}{p^2\ln^2 p}+\frac{n}{(1-p)^2}+\frac{n}{(1-p)^2\ln p}(\frac{3}{2}+\frac{1}{2p^2}-\frac{2}{p}-\ln p);\\ &&J_{12}=J_{21}=n\frac{1-p+p\ln p}{2\beta p(1-p)\ln p};\\ &&J_{22}=-n\frac{g(p)}{\beta^{2}\ln p}.\end{eqnarray*}$

于是$J(p, \beta)^{-1}$在$\widehat{\theta}$的值就给出了参数$\theta$极大似然估计的渐近协方差阵, 记为$\widehat{J^{-1}}(\widehat{p}, \widehat{\beta})$.

以下给定不同初始值, 分别在不同的样本容量下进行1000次Monte-Carlo模拟(迭代容差为0.0001), 最后用EM算法求得$(p, \beta)$的极大似然估计.

假设1000次模拟所得到的的参数的估计值为$\widehat{p_{i}}, \widehat{\beta_{i}}$, $i=1, 2, \cdots, 1000$.

记样本均值为

${\rm MEAN}(p)=\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000} \widehat{p_{i}}, {\rm MEAN}(\beta)=\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000}\widehat{\beta_{i}};$

样本标准差为

${\rm SD}(p)=\sqrt{\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000} ({\rm MEAN}(p)-\widehat{p_{i}})^{2}}, {\rm SD}(\beta)=\sqrt{\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000}({\rm MEAN}(\beta)-\widehat{\beta_{i}})^{2}} ;$

样本均方误差为

${\rm MSE}(p)=\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000}(\widehat{p_{i}}-p)^{2}, {\rm MSE}(\beta)=\frac{1}{1000}\mathop \sum \limits_{i = 1}^{1000}(\widehat{\beta_{i}}-\beta)^{2}.$

多元正态分布$N_{2}(0, \widehat{J^{-1}}(\widehat{p}, \widehat{\beta}))$可以用来构造参数$p, \beta$的渐近置信区间.其中$p$的水平为$1-\eta$的渐近置信区间为$\widehat{p}\pm z_{\frac{\eta}{2}}[\widehat{{\rm var}}(\widehat{p})]^{\frac{1}{2}}$, $\beta$的$1-\eta$渐近置信区间为$\widehat{\beta}\pm z_{\frac{\eta}{2}}[\widehat{{\rm var}}(\widehat{\beta})]^{\frac{1}{2}}$.这里$\widehat{{\rm var}}(\cdot)$为$\widehat{J^{-1}}(\widehat{p}, \widehat{\beta})$中对应参数对角元素的值, 而$z_{\frac{\eta}{2}}$表示标准正态分布的$1-\frac{\eta}{2}$分位点.在这1000次模拟中, $\alpha_{1}$为$p$的$95\%$置信区间的实际覆盖率, $\alpha_{2}$为$\beta$的$95\%$置信区间的实际覆盖率.

从表 1可以看出, 用EM算法求得的极大似然估计可以很好的反应参数的真实值, 选定不同的初值对估计的影响不大, 并且随着样本容量$n$的增大, 估计的偏差越来越小, 估计的精度就越来越高.当$p$固定时, $\beta$越小, 极大似然估计对$\beta$的估计效果越好, 对$p$估计变化效果不明显.当$\beta$固定时, $p$越小, 极大似然估计对$p$的估计效果越好, 对$\beta$估计变化效果不明显.总的来说, 当$p, \beta$的值越小, MLE估计的效果越好.