Kumaraswamy Marshall-Olkin Rayleigh分布的统计分析

Statistical Analysis of Kumaraswamy Marshall-Olkin Rayleigh Distributions

摘要: 文章通过Kumaraswamy Marshall-Olkin扩展，将Rayleigh分布扩展为一个新的具有四参数的分布Kumaraswamy Marshall-Olkin Rayleigh (KwMO-R)分布，然后计算了该分布的矩、危险率函数、Renyi熵、次序统计量的极限分布以及参数的极大似然估计，并验证了估计的相合性，最后利用Monte-Carlo算法对参数的估计进行了数值模拟。

Abstract: In this paper, through the Kumaraswamy Marshall-Olkin extension, the Rayleigh distribution is extended to a new distribution with four parameters. Then, the moment, risk function, Renyi entropy, limit distribution of order statistics and maximum likelihood estimation of the parameters are calculated, and the consistency of the estimation is verified. Finally, the parameter estimation is simulated by Monte-Carlo algorithm.

1. 引言

可靠性统计研究的核心是基于寿命分布的分析。然而，由于各种实际因素和产品特性的差异，寿命数据展现出显著的多样性。这种多样性体现在数据的类型、分布形态、变化趋势等多个方面，使得单一的寿命分布模型往往难以全面、准确地描述实际数据的特征。在过去的一段时间里，许多国内外学者致力于研究通过在原有基础连续分布基础上增加尺度函数，得到一种新的扩展分布族。增加的参数已被证明在探索偏度和尾部是很有用的，也提高了新分布的拟合优度。Eugene et al. (2002) and Jones (2004)提出了The beta-G扩展分布族，得到了Beta Weibull-geometric分布，可应用于生物领域的数据建模分析[1] [2]；Marshall和Olkin首次利用分布生成技术，提出了一种通过增加形状参数得到扩展分布族的方法：M-O扩展，并成功应用于指数分布和威布尔分布中[3]；Cordeiroandde Castro (2011)提出了The Kumaraswamy-G (Kw-G)扩展分布族，并应用Gumbel分布[4]；Ramos (2014)提出了The Kumaraswamy Poisson-G (Kw-G)扩展分布，对应的概率密度函数和相应的风险率函数更灵活[5]；刘焱哲等通过Kumaraswamy Marshall-Olkin扩展方法引进并研究了一个新的五参数寿命分布[6]；常帅等利用Kumaraswamy分布结构推广到了倒帕累托分布[7]；这些研究都是通过引入额外的参数，定义了新的多参数分布，同时研究了它们的各阶矩、熵、参数的估计，并给出了新分布参数估计值的算法。

本文将KwMO方法对Rayleigh分布进行扩展，得到KwMO-R分布，对研究在无线电通信工程、工程测量等领域有广泛的意义。在已有的Rayleigh分布研究基础上，本文进一步讨论了新分布的分位数、矩、次序统计量的性质、Renyi熵、并验证了估计的相合性，最后进行了数值模拟。

2. 分布的定义

瑞利分布的分布函数和概率密度函数分别为：

$G\left(x;\sigma \right)=1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{I}_{\left(0,\infty \right)}\left(x\right)$ (2.1)

$g_R\left(x;\sigma \right)=\frac{x}{{\sigma }^2}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{I}_{\left(0,\infty \right)}\left(x\right)$ (2.2)

其中$\sigma >0$ 为尺度参数。

引理：设变量X的分布函数为$G\left(x;\xi \right)$ ，密度函数为$g\left(x;\xi \right)$ ，则KwMO方法扩展得到的新分布族为：

$F\left(x;a,b,p,\xi \right)=1-{\left\{1-{\left[\frac{G\left(x\right)}{1-p\overline{G}\left(x\right)}\right]}^a\right\}}^b$

对应的密度函数为$f\left(x;a,b,p,\xi \right)=\frac{ab\left(1-p\right)g\left(x\right)G{\left(x\right)}^{a-1}}{{\left[1-p\overline{G}\left(x\right)\right]}^{a+1}}{\left\{1-{\left[\frac{G\left(x\right)}{1-p\overline{G}\left(x\right)}\right]}^a\right\}}^{b-1}$

根据引理把(2.1)、(2.2)分别代入上述两式，得到：

$F_{\text{KwMO-R}}\left(x;\theta \right)=1-{\left\{1-{\left[\frac{1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}\right]}^a\right\}}^b$ (2.3)

其密度函数为 $f_{\text{KwMO-R}}\left(x;\theta \right)=\frac{abx\left(1-p\right){\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\left[1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right]}^{a-1}}{{\sigma }^2{\left(1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{a+1}}{\left\{1-{\left[\frac{1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}\right]}^a\right\}}^{b-1}$ (2.4)

定义2.1 称(2.3)或(2.4)式给出的分布为KwMO-R分布，记为$X~\text{KwMO-R}\left(x;a,b,p,\sigma \right)$ ，$a>0,b>0,0<p<1$ 为形状参数，$\sigma >0$ 为尺度参数。

(2.4)式也可以写成：$f_{\text{KwMO-R}}\left(x;\theta \right)=\frac{ab\cdot \frac{1}{2{\sigma }^2}\cdot 2\cdot x^{2-1}\left(1-p\right){\text{e}}^{-\frac{1}{2{\sigma }^2}{x}^2}{\left[1-{\text{e}}^{-\frac{1}{2{\sigma }^2}{x}^2}\right]}^{a-1}}{{\left(1-p{\text{e}}^{-\frac{1}{2{\sigma }^2}{x}^2}\right)}^{a+1}}{\left\{1-{\left[\frac{1-{\text{e}}^{-\frac{1}{2{\sigma }^2}{x}^2}}{1-p{\text{e}}^{-\frac{1}{2{\sigma }^2}{x}^2}}\right]}^a\right\}}^{b-1}$ ，即$\text{KwMO-W}\left(a,b,p,\frac{1}{2{\sigma }^2},2\right)$ 。

(2.4)式中$a=1,b=1$ 时，该分布为$RG\left(p,\sigma \right)$ 。

$f_{\text{KwMO-R}}\left(x;\theta \right)$ 关于x递减，且当$x\to \infty$ 时，$f\left(x\right)\to 0$ ；$x\to 0$ 时，$f\left(x\right)\to 0$ 。

3. 分布的性质

定理3.1

1) KwMO-R分布的$q\left(0<q<1\right)$ 分位数为$Q\left(q\right)=\sqrt{2}\sigma {\left\{\log \left(\frac{1-p{\left[1-{\left(1-u\right)}^{\frac{1}{b}}\right]}^{\frac{1}{a}}}{1-{\left[1-{\left(1-u\right)}^{\frac{1}{b}}\right]}^{\frac{1}{a}}}\right)\right\}}^{\frac{1}{2}}$ ，即中位数$Q_{0.5}=\sqrt{2}\sigma {\left\{\log \left(\frac{1-p{\left[1-{0.5}^{\frac{1}{b}}\right]}^{\frac{1}{a}}}{1-{\left[1-{0.5}^{\frac{1}{b}}\right]}^{\frac{1}{a}}}\right)\right\}}^{\frac{1}{2}}$ 。

2) KwMO-R分布的r阶中心距$E\left(X^r\right)={\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i}}{s}_{k+1}\left(\begin{array}{c}k\\ i\end{array}\right)\frac{\left(k+1\right){\left(\sqrt{2}\sigma \right)}^r}{{\left(i+1\right)}^{\frac{r}{2}+1}}\Gamma \left(\frac{r}{2}+1\right)$ ，特别地KwMO-R的数学期望为$E\left(X\right)=\frac{\sqrt{2\pi }\sigma }{2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i{\left(i+1\right)}^{-\frac{3}{2}}\left(k+1\right)\left(\begin{array}{c}k\\ i\end{array}\right)s_{k+1}}}$ ，其中$s_{k+1}$ 见下面所述。

证明：1) 由KwMO-R分布的定义及分位数的定义可得，略。

2) 由$F\left(x;a,b,p,\xi \right)=1-{\left\{1-{\left[\frac{G\left(x\right)}{1-p\overline{G}\left(x\right)}\right]}^a\right\}}^b$

由广义二项式展开定理可得：

$\begin{array}{c}F\left(x;a,b,p,\xi \right)=1-{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i\left(\begin{array}{c}b\\ i\end{array}\right)}{\left[\frac{G\left(x\right)}{1-p\overline{G}\left(x\right)}\right]}^{ai}\\ =1-{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i\left(\begin{array}{c}b\\ i\end{array}\right)}\frac{{\displaystyle \sum _{k=0}^{\infty }{u}_{k}G{\left(x\right)}^k}}{{\displaystyle \sum _{k=0}^{\infty }{v}_{k}G{\left(x\right)}^k}}\\ =1-{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i\left(\begin{array}{c}b\\ i\end{array}\right)}{\displaystyle \sum _{k=0}^{\infty }{l}_k}G{\left(x\right)}^k\end{array}$

其中，$u_k={\displaystyle \sum _{j=k}^{\infty }{\left(-1\right)}^{j+k}}\left(\begin{array}{c}ai\\ j\end{array}\right)\left(\begin{array}{c}j\\ k\end{array}\right)$ ，$v_k={\displaystyle \sum _{j=k}^{\infty }{\left(-1\right)}^{j+k}}{p}^j\left(\begin{array}{c}ai\\ j\end{array}\right)\left(\begin{array}{c}j\\ k\end{array}\right)$ ，$l_k=\frac{1}{v_0}\left(u_k-\frac{1}{v_0}{\displaystyle \sum _{s=1}^k{u}_s}{v}_{k-s}\right),\left(k\ge 1\right)$ ，$l_0=u_0/v_0$

当$k\ge 1$ 时，$F\left(x;a,b,p,\xi \right)={\displaystyle \sum _{k=0}^{\infty }{s}_{k}G{\left(x\right)}^k}={\displaystyle \sum _{k=0}^{\infty }{s}_k{H}_k\left(x\right)}$ (3.1)

这里$s_k=-{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i}\left(\begin{array}{c}b\\ i\end{array}\right)l_k,s_0=1-{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i}\left(\begin{array}{c}b\\ i\end{array}\right)l_0$

$G{\left(x\right)}^k$ 是参数为k的广义指数分布的分布函数，即$H_k\left(x\right)=G{\left(x\right)}^k={\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^k{I}_{\left(0,\infty \right)}\left(x\right)$ 。

类似地，KwMO-R分布的密度函数：

$f\left(x;a,b,p,\xi \right)={\displaystyle \sum _{k=0}^{\infty }{s}_{k+1}{h}_{k+1}\left(x\right)}={\displaystyle \sum _{k=0}^{\infty }{s}_{k+1}g{\left(x\right)}^{k+1}},\left(k\ge 1\right)$ (3.2)

$g{\left(x\right)}^k$ 是参数为k的广义指数分布的密度函数，即$g{\left(x\right)}^k=\frac{kx}{{\sigma }^2}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{k-1}{I}_{\left(0,\infty \right)}\left(x\right)$

故KwMO-R分布的r阶中心距为：

$\begin{array}{c}E\left(X^r\right)={\displaystyle {\int }_0^{\infty }{x}^r}{\displaystyle \sum _{k=0}^{\infty }{s}_{k+1}}g{\left(x\right)}^{k+1}\text{d}x\\ ={\displaystyle \sum _{k=0}^{\infty }{s}_{k+1}}{\displaystyle {\int }_0^{\infty }{x}^r\frac{k+1}{{\sigma }^2}x{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i}}\left(\begin{array}{c}k\\ i\end{array}\right){\left({\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^i\text{d}x\\ ={\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i\left(\begin{array}{c}k\\ i\end{array}\right)s_{k+1}}}\frac{k+1}{{\sigma }^2}{\displaystyle {\int }_0^{\infty }{x}^{r+1}}{\text{e}}^{-\frac{x^2\left(i+1\right)}{2{\sigma }^2}}\text{d}x\\ ={\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i}}{s}_{k+1}\left(\begin{array}{c}k\\ i\end{array}\right)\frac{\left(k+1\right){\left(\sqrt{2}\sigma \right)}^r}{{\left(i+1\right)}^{\frac{r}{2}+1}}\Gamma \left(\frac{r}{2}+1\right)\end{array}$ (3.3)

当$r=1$ 时，KwMO-R的数学期望为：

$E\left(X\right)=\frac{\sqrt{2\pi }\sigma }{2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{i=0}^{\infty }{\left(-1\right)}^i{\left(i+1\right)}^{-\frac{3}{2}}\left(k+1\right)\left(\begin{array}{c}k\\ i\end{array}\right)s_{k+1}}}$ (3.4)

根据生存函数、危险率函数的定义，得到$\text{KwMO-R}\left(x;a,b,p,\sigma \right)$ 分布的生存函数为：

$S_{\text{KwMO-R}}\left(x;a,b,p,\sigma \right)=F_{\text{KwMO-R}}\left(x;a,b,p,\sigma \right)={\left\{1-{\left[\frac{1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}\right]}^a\right\}}^b$ (3.5)

定理3.2 若随机变量$X~\text{KwMO-R}\left(x;a,b,p,\sigma \right)$ 分布，则其危险率函数为$h\left(x;a,b,p,\sigma \right)=\frac{ab\left(1-p\right)x{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{a-1}}{{\sigma }^2\left(1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)\left[{\left(1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^a-{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^a\right]}$ $\left(x>0\right)$ ，并且有$\underset{x\to 0}{\lim }h\left(x;a,b,p,\sigma \right)=0$ ，$\underset{x\to +\infty }{\lim }h\left(x;a,b,p,\sigma \right)=\frac{bx}{{\sigma }^2}$ ，同时$h\left(x;a,b,p,\sigma \right)$ 所有参数关于x都是倒浴盆曲线，

证明：根据分布函数和危险率函数的定义可得：

$h\left(x;a,b,p,\sigma \right)=\frac{ab\left(1-p\right)x{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{a-1}}{{\sigma }^2\left(1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)\left[{\left(1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^a-{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^a\right]}$ 。

对上式利用洛必达法则易得上述极限。证毕。

当$a>1,b<1$ 时，函数先增加后趋于平缓，呈J形；当$a<1,b>1$ 时，随着$\sigma$ 的增大，曲线呈浴盆型，p越大，函数越凸。

定理3.3 若$X_1,X_2,\cdots ,X_n$ 是来自KwMO-R的简单随机样本，则第i阶次序统计量为$f_{\left(i\right)}\left(x\right)={\displaystyle \sum _{r,k=0}^{\infty }{\displaystyle \sum _{j=0}^{n-i}{m}_{r,k,j}}}g{\left(x\right)}^{k+r+1}$ ，其中$m_{r,k,j}=\frac{{\left(-1\right)}^{j}n!\left(\begin{array}{c}n-i\\ j\end{array}\right)\left(r+1\right)s_{r+1}{d}_{j+i-1,k}}{\left(i-1\right)!\left(n-i\right)!\left(k+r+1\right)}$ ；记$X_{\left(n\right)}=\max \left\{X_1,X_2,\cdots ,X_n\right\}$ ，则$\underset{n\to \infty }{\lim }P\left(X_{\left(n\right)}\le b_{n}x\right)=\exp \left(-x^2\right)$ $\left(x>0\right)$ ，其中$b_n=F^{-1}\left(1-\frac{1}{n}\right)$ 。

证明：由次序统计量的定义，可得第i阶次序统计量的密度函数为：

$\begin{array}{c}{f}_{\left(i\right)}\left(x;a,b,p,\sigma \right)=\frac{n!}{\left(i-1\right)!\left(n-i\right)!}{\left[F\left(x\right)\right]}^{i-1}\cdot {\left[1-F\left(x\right)\right]}^{n-i}\cdot f\left(x;a,b,p,\sigma \right)\\ =\frac{n!}{\left(i-1\right)!\left(n-i\right)!}{\displaystyle \sum _{j=0}^{n-i}{\left(-1\right)}^j}\left(\begin{array}{c}n-i\\ j\end{array}\right)f\left(x;a,b,p,\sigma \right){\left[F\left(x;a,b,p,\sigma \right)\right]}^{i+j-1}\end{array}$

根据${\left({\displaystyle \sum _{i=0}^{\infty }{a}_i}{u}^i\right)}^n={\displaystyle \sum _{i=0}^{\infty }{c}_{n,i}{u}^i}$ ，$\left(i=1,2,\cdots \right)$ ，这里$c_{n,i}={\left(i{a}_0\right)}^{-1}{\displaystyle \sum _{m=1}^i\left[m\left(n+1\right)-i\right]a_m{c}_{n,i-m}}$ ，$c_{n,0}=a_0^n$ (见参考文献[8])，把(3.1)和(3.2)代入上式，得到：

$\begin{array}{c}\left(x\right)=\frac{n!}{\left(i-1\right)!\left(n-i\right)!}{\displaystyle \sum _{j=0}^{n-i}{\left(-1\right)}^j}\left(\begin{array}{c}n-i\\ j\end{array}\right)\left[{\displaystyle \sum _{r=0}^{\infty }{s}_{r+1}}\left(r+1\right)G{\left(x\right)}^{r}g\left(x\right)\right]{\left[{\displaystyle \sum _{k=0}^{\infty }{s}_k}G{\left(x\right)}^k\right]}^{j+i-1}\\ =\frac{n!}{\left(i-1\right)!\left(n-i\right)!}{\displaystyle \sum _{j=0}^{n-i}{\left(-1\right)}^j}\left(\begin{array}{c}n-i\\ j\end{array}\right)\left[{\displaystyle \sum _{r=0}^{\infty }{s}_{r+1}}\left(r+1\right)G{\left(x\right)}^{r}g\left(x\right)\right]\left[{\displaystyle \sum _{k=0}^{\infty }{d}_{j+i-1,k}}G{\left(x\right)}^k\right]\\ ={\displaystyle \sum _{r,k=0}^{\infty }{\displaystyle \sum _{j=0}^{n-i}\frac{{\left(-1\right)}^{j}n!\left(\begin{array}{c}n-i\\ j\end{array}\right)\left(r+1\right)s_{r+1}{d}_{j+i-1,k}}{\left(i-1\right)!\left(n-i\right)!\left(k+r+1\right)}g{\left(x\right)}^{k+r+1}}}\\ ={\displaystyle \sum _{r,k=0}^{\infty }{\displaystyle \sum _{j=0}^{n-i}{m}_{r,k,j}}}g{\left(x\right)}^{k+r+1}\end{array}$

其中$m_{r,k,j}=\frac{{\left(-1\right)}^{j}n!\left(\begin{array}{c}n-i\\ j\end{array}\right)\left(r+1\right)s_{r+1}{d}_{j+i-1,k}}{\left(i-1\right)!\left(n-i\right)!\left(k+r+1\right)}$ ，$g{\left(x\right)}^{k+r+1}$ 是参数为$k+r+1$ 广义指数分布的密度函数，即$g{\left(x\right)}^{k+r+1}=\frac{\left(k+r+1\right)x}{{\sigma }^2}{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{k+r}{I}_{\left(0,\infty \right)}\left(x\right)$ 。(见参考文献[9])

由(2.3)，(2.4)式有$\underset{x\to 0}{\lim }\left(1-F\left(x\right)\right)={\left[a\left(1-p\right)\exp \left(-\frac{x^2}{2{\sigma }^2}\right)\right]}^b$

$\underset{x\to 0}{\lim }f\left(x\right)=b{\left[a\left(1-p\right)\right]}^b\frac{x}{{\sigma }^2}{\left[\exp \left(-\frac{x^2}{2{\sigma }^2}\right)\right]}^b$

根据洛必达法则有：

$\underset{t\to 0}{\lim }\frac{F\left(tx\right)}{F\left(t\right)}=\underset{t\to 0}{\lim }\frac{xf\left(tx\right)}{f\left(t\right)}=\frac{xb{\left[a\left(1-p\right)\right]}^b\frac{tx}{{\sigma }^2}{\left({\text{e}}^{-\frac{t^2{x}^2}{2{\sigma }^2}}\right)}^b}{b{\left[a\left(1-p\right)\right]}^b\frac{t}{{\sigma }^2}{\text{e}}^{-\frac{t^2}{2{\sigma }^2}}{\left({\text{e}}^{-\frac{t^2}{2{\sigma }^2}}\right)}^b}=x^2$

根据文献[10]定理2.1.2及定理2.4.3可知，当$b_n=F^{-1}\left(1-\frac{1}{n}\right)$ 时，有$\underset{n\to \infty }{\lim }P\left(X_{\left(n\right)}\le b_{n}x\right)=\exp \left(-x^2\right)$ 。

定理3.4 若随机变量$X~\text{KwMO-R}\left(x;a,b,p,\sigma \right)$ 分布，则其Renyi熵为$h_{\gamma }=\frac{1}{1-\gamma }\log \left\{{\left({\sigma }^2\right)}^{-\gamma +1}{\displaystyle \sum _{k=0}^{\infty }{\left[\left(k+1\right)s_{k+1}\right]}^{\gamma }B\left(\gamma ,k\gamma +1\right)}\right\}$ ，其中$\gamma >0$ 且$\gamma \ne 1$ 。

证明：根据Renyi熵$h_{\gamma }=\frac{1}{1-\gamma }\log \left\{{\displaystyle {\int }_{ℝ}{f}^{\gamma }\left(x\right)\text{d}x}\right\}$ 和(3.2)式可得：

${\displaystyle {\int }_0^{+\infty }{f}^{\gamma }\left(x\right)\text{d}x}={\displaystyle {\int }_0^{\infty }{\left[\frac{1}{{\sigma }^2}{\displaystyle \sum _{k=0}^{\infty }\left(k+1\right)s_{k+1}x{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^k\right]}^{\gamma }\text{d}x}$

若$\gamma \ge 1$ 有

$\begin{array}{c}{\displaystyle {\int }_0^{+\infty }{f}^{\gamma }\left(x\right)\text{d}x}={\left({\sigma }^2\right)}^{-\gamma }{\displaystyle \sum _{k=0}^{\infty }{\left(k+1\right)}^{\gamma }{s}_{k+1}^{\gamma }}{\displaystyle {\int }_0^{\infty }x{\text{e}}^{-\frac{\gamma x^2}{2{\sigma }^2}}{\left(1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right)}^{k\gamma }\text{d}x}\\ ={\left({\sigma }^2\right)}^{-\gamma +1}{\displaystyle \sum _{k=0}^{\infty }{\left(k+1\right)}^{\gamma }{s}_{k+1}^{\gamma }}{\displaystyle {\int }_0^1{y}^{\gamma -1}{\left(1-y\right)}^{k\gamma }\text{d}y}\\ ={\left({\sigma }^2\right)}^{-\gamma +1}{\displaystyle \sum _{k=0}^{\infty }{\left[\left(k+1\right)s_{k+1}\right]}^{\gamma }B\left(\gamma ,k\gamma +1\right)}\end{array}$

其中$s_k$ 见(3.1)

故$h_{\gamma }=\frac{1}{1-\gamma }\log \left\{{\left({\sigma }^2\right)}^{-\gamma +1}{\displaystyle \sum _{k=0}^{\infty }{\left[\left(k+1\right)s_{k+1}\right]}^{\gamma }B\left(\gamma ,k\gamma +1\right)}\right\}$ 。

4. 参数估计

设$X_1,X_2,\cdots ,X_n$ 是来自$\text{KwMO-R}\left(x;\theta \right)$ $\theta ={\left(a,b,p,\sigma \right)}^{\prime }$ 的简单随机样本，记$x_{obs}=\left(x_1,\cdots ,x_n\right)$ 为样本观测值，其对数似然函数$l\left(\theta ;x_{obs}\right)$ 为：

$\begin{array}{c}l\left(\theta ;x_{obs}\right)=n\ln \left[ab\left(1-p\right)\right]+{\displaystyle \sum _{i=1}^n\ln g\left(x_i\right)}+\left(a-1\right){\displaystyle \sum _{i=1}^n\ln G\left(x_i\right)}\\ -\left(a+1\right){\displaystyle \sum _{i=1}^n\ln \left[1-p\overline{G}\left(x_i\right)\right]}+\left(b-1\right){\displaystyle \sum _{i=1}^n\ln \left\{1-{\left[\frac{G\left(x_i\right)}{1-p\overline{G}\left(x_i\right)}\right]}^a\right\}}\end{array}$

似然方程为：

$\frac{\partial l}{\partial a}=\frac{n}{a}+{\displaystyle \sum _{i=1}^n\ln }\left[\frac{1-{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}\right]+\left(1-b\right){\displaystyle \sum _{i=1}^n\frac{{\left[\frac{1-e^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{e}^{-\frac{x_i^2}{2{\sigma }^2}}}\right]}^a\ln \left[\frac{1-{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}\right]}{1-{\left[\frac{1-{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}\right]}^a}}=0$ (4.1)

$\frac{\partial l}{\partial b}=\frac{n}{b}+{\displaystyle \sum _{i=1}^n\ln \left\{1-{\left[\frac{1-{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}\right]}^a\right\}}=0$ (4.2)

$\frac{\partial l}{\partial p}=\frac{-n}{1-p}+\left(a+1\right){\displaystyle \sum _{i=1}^n\frac{1-{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}{1-p{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}}=0$ (4.3)

$\begin{array}{l}\frac{\partial l}{\partial \sigma }={\displaystyle \sum _{i=1}^n\frac{g^{\left(\sigma \right)}\left(x_i\right)}{\frac{x_i}{{\sigma }^2}{\text{e}}^{-\frac{x_i^2}{2{\sigma }^2}}}}+\left(a-1\right){\displaystyle \sum _{i=1}^n\frac{G^{\left(\sigma \right)}\left(x_i\right)}{1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}}+p\left(a-1\right){\displaystyle \sum _{i=1}^n\frac{G^{\left(\sigma \right)}\left(x_i\right)}{1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}}}\\ +a\left(1-p\right){\displaystyle \sum _{i=1}^n\frac{G^{\left(\sigma \right)}\left(x_i\right){\left[1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right]}^{a-1}}{\left[1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right]\left\{{\left[1-p{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right]}^a-{\left[1-{\text{e}}^{-\frac{x^2}{2{\sigma }^2}}\right]}^a\right\}}}=0\end{array}$ (4.4)

解方程组(4.1)~(4.4)，其解$\left(\widehat{a},\widehat{b},\widehat{p},\widehat{\sigma }\right)$ 就为参数$\left(a,b,p,\sigma \right)$ 的极大似然估计(MLE)。

下面讨论估计的相合性。

$\text{KwMO-R}\left(x;\theta \right)$ 分布的参数空间记为$\Theta =\left\{\theta =\left(a,b,p,\sigma \right):p\in \left(0,1\right),a\in \left(0,+\infty \right),b\in \left(0,+\infty \right),\sigma \in \left(0,+\infty \right)\right\}$

对任意${\theta }_0\in \Theta$ 及$\epsilon >0$ ，记：

$h\left(x,\theta \right)=\left(\frac{\partial \ln f\left(x\right)}{\partial a},\frac{\partial \ln f\left(x\right)}{\partial b},\frac{\partial \ln f\left(x\right)}{\partial p},\frac{\partial \ln f\left(x\right)}{\partial \sigma }\right)$

$\left|h\left(x,\theta \right)\right|=\left|\frac{\partial \ln f\left(x\right)}{\partial a}\right|+\left|\frac{\partial \ln f\left(x\right)}{\partial b}\right|+\left|\frac{\partial \ln f\left(x\right)}{\partial p}\right|+\left|\frac{\partial \ln f\left(x\right)}{\partial \sigma }\right|$

$H\left(x,{\theta }_0,\epsilon \right)=\sup \left\{\left|h\left(x,\theta \right)\right|:\theta \in B\left({\theta }_0,\epsilon \right)\right\}$ ，这里$B\left({\theta }_0,\epsilon \right)=\left\{\theta :\left|\theta -{\theta }_0\right|\le \epsilon \right\}$

根据文献[11]命题2.4.34，只需要证明对充分小的$\epsilon >0$ ，有$E_{{\theta }_0}\left[H\left(x,{\theta }_0,\epsilon \right)\right]<\infty$ ，MLE就存在。事实上，在$\Theta$ 上，$h\left(x,\theta \right)$ 为$\theta$ 的连续函数，故存在${\theta }_1=\left(a_1,b_1,p_1,{\sigma }_1\right)\in B\left({\theta }_0,\epsilon \right)$ ，使得$\left|h\left(x,{\theta }_1\right)\right|=H\left(x,{\theta }_0,\epsilon \right)$

$E_{{\theta }_0}\left[H\left(X,{\theta }_0,\epsilon \right)\right]=E_{{\theta }_0}\left[{\left|\frac{\partial \ln f\left(X\right)}{\partial a}\right|}_{{\theta }_1}\right]+E_{{\theta }_0}\left[{\left|\frac{\partial \ln f\left(X\right)}{\partial b}\right|}_{{\theta }_1}\right]+E_{{\theta }_0}\left[{\left|\frac{\partial \ln f\left(X\right)}{\partial p}\right|}_{{\theta }_1}\right]+E_{{\theta }_0}\left[{\left|\frac{\partial \ln f\left(X\right)}{\partial \sigma }\right|}_{{\theta }_1}\right]$

$\begin{array}{c}\ln f\left(x\right)=\ln \left[ab\left(1-p\right)\right]+\ln g\left(x\right)+\left(a-1\right)\ln G\left(x\right)-\left(a+1\right)\ln \left(1-p\overline{G}\left(x\right)\right)\\ +\left(b-1\right)\ln \left\{1-{\left(\frac{G\left(x\right)}{1-p\overline{G}\left(x\right)}\right)}^a\right\}\end{array}$

由文献[12]定理2.1.3可知，若$\psi \left(X\right)={\left[\frac{G\left(X\right)}{1-p\overline{G}\left(X\right)}\right]}^a$ ，可以证明$E\left\{{\left[\frac{G\left(X\right)}{1-p\overline{G}\left(X\right)}\right]}^a|X\ge x\right\}=\frac{1}{b+1}\left\{b{\left[\frac{G\left(X\right)}{1-p\overline{G}\left(X\right)}\right]}^a+1\right\}$ ，即$E\left(\psi \left(X\right)\right)<\infty$ 存在。由定理3.1可以证明$E\left(\left|\ln g\left(x\right)\right|\right)<\infty$ ，$E\left(\left|\ln G\left(x\right)\right|\right)<\infty$ 。

故$E_{{\theta }_0}\left[H\left(x,{\theta }_0,\epsilon \right)\right]<\infty$ 成立，似然方程必存在一个强相合解。

5. 数值模拟

由于上述似然方程组不易求出显示解的表达式，考虑用Monte-Carlo模拟算法计算其极大似然估计。

设模型参数的真值分别为$a=4,b=1.2,p=0.2,\sigma =3.5$ 。

1) 确定需要产生的样本容量$n=50,100,500$ ；

2) 产生n个独立随机数$u_1,u_2,\cdots ,u_n~\text{Unif}\left(0,1\right)$ ，计算$x_i=F^{-1}\left(u_i;\theta \right),i=1,2,\cdots ,n$ ，则$x_{obs}=\left(x_1,x_2,\cdots ,x_n\right)$ 为$\text{KwMO-R}\left(x;\theta \right)$ 分布的容量为n的样本。

3) 对给定参数的初始值${\theta }^{\left(0\right)}=\left(a^{\left(0\right)},b^{\left(0\right)},p^{\left(0\right)},{\sigma }^{\left(0\right)}\right)$ ，设定迭代的次数为5000，代入迭代公式(4.1)和(4.4)得参数到$\left(a,b,p,\sigma \right)$ 的均值、偏差、均方误差，见表1。

Table 1. Simulation results of KwMO-R distribution parameters

表1. KwMO-R分布参数的模拟结果

n	Parameter	${\theta }^{\left(0\right)}$	MEAN	Bias	MSN
50	a	4	3.301337	−0.698663	9.4057169
	b	1.2	1.186861	0.186861	2.2464300
	p	0.2	0.3233599	0.1233599	0.1805369
	$\sigma$	3.5	3.406095	−0.093905	5.3471411
100	a	4	3.422929	−0.577071	7.006185
	b	1.2	1.01385	0.01385	1.5863938
	p	0.2	0.3383646	0.1383646	0.1821119
	$\sigma$	3.5	3.524606	0.024606	5.1935504
500	a	4	3.533069	−0.466931	5.811375
	b	1.2	0.9443207	−0.055679	0.4228819
	p	0.2	0.3332891	0.1332891	0.1743206
	$\sigma$	3.5	3.494839	−0.005161	4.667971

从表中可以看出随着样本容量n的增加，虽然p的估计值偏差增加，但是均方误差值是减小的；总体来说随着样本n的增加，$a,b,p,\sigma$ 估计值的偏差和均方误差都是递减渐近趋于0的。

6. 结论

本文针对Rayleigh分布进行扩展，得到一种广义四参数$\text{KwMO-R}\left(x;a,b,p,\sigma \right)$ 分布，该分布的危险函数主要呈J型、浴盆型，因此它可以灵活地拟合一些较复杂的寿命数据；同时对该分布的进行了统计分析，最后得到了未知参数的估计公式，并验证了估计的强相合性；在MC模拟样本下，随着样本容量的增加，未知参数估计值的偏差和均方误差都是趋于0的。

基金项目

2021年湖北省教育厅科学研究计划资助项目(B2021286)；2022年湖北省高等学校优秀中青年科技创新团队计划项目(T2022035)。