基于混合Gamma分布单个观测的可靠性函数无偏估计

Unbiased Estimation of Reliability Function Based on One Single Observation from a Mixture of Gamma Distributions

摘要

针对小样本甚至单样本情况下的疲劳寿命分布特征研究，给出了混合Gamma分布的概率密度函数可定义的充要条件，基于混合Gamma分布的单个观测值，得到了Gamma分布的可靠性函数的两个无偏估计，并讨论了所得无偏估计的适当性， 最后进行了数值模拟验证。

Abstract

In order to study the characteristics of fatigue life distributions in small sample or one sample，we gave the sufficient and necessary condition for the definition of probability density function of mixed Gamma distribution. Then we obtained two unbiased estimations of the reliability function based on one single observation from a mixture of two Gamma distributions. Next we discussed the appropriateness of the obtained unbiased estimations. Finally, we did some simulations.

Graphical abstract

关键词

混合分布 / 污染数据 / 可靠性函数 / Gamma分布 / 无偏估计

Key words

mixture distribution / contaminated data / reliability function / Gamma distribution / unbiased estimation

0　 引 言

在传统的研究中，人们一般首先假设样本数据来自某个特定的分布总体，然后基于这个总体进行后续的统计推断。但随着科技发展，我们所面临的数据越来越复杂。在很多情形下，样本数据在不同局部可能有不同的特性，单一的参数分布族无法确切地描述观测数据，人们想到用混合分布模型对复杂数据进行统计分析。混合分布模型最早可追溯到Newcomb^[1]和Pearson^[2]的相关工作。近三十几年来，混合分布模型发展迅速，应用广泛，尤其是在生物分析、资源评价、金融股票等领域^[3]。比较典型的例子有^[4]:l) 渔业研究中，鱼的体长可以度量，但其性别和年龄不易确定，需要研究如何由体长数据确定不同年龄组和不同性别鱼群的比例；2) 沉积岩学中，一个沙样中颗粒大小的分布已知，但不知道其中包含有哪些不同种类的矿石；3) 医学诊断中可对一组病人做临床调查，但他们所患疾病属于哪一种类型是未知的；4) 在RFLP-QTL连锁分析中，可以知道每个个体的标记基因型及性状值，这个性状值是由三种类型的QTL基因型混合而成，并且个体的QTL基因型未知，如何由标记基因型的表型去推断QTL与标记之间的连锁以及QTL效应值。此外，混合分布模型在寿命数据分析^[5]和污染数据分析中也扮演着重要角色。早在1952年，Davis^[6]就注意到在寿命实验中，元件寿命分布函数可能为两个分布函数的混合。

考虑独立同分布的随机变量列X_{\mathrm{1}},\dots ,X_n，具有分布函数

F_{\alpha }(x)=(\mathrm{1}-\alpha )F_{\mathrm{1}}(x)+\alpha F_{\mathrm{2}}(x)(1)

其中F_{\mathrm{1}}(x),F_{\mathrm{2}}(x)都是分布函数。通常我们更关心F_{\mathrm{1}}(x)，认为数据样本应该服从F_{\mathrm{1}}(x)，但却受到了来自分布F_{\mathrm{2}}(x)的数据的污染，称F_{\mathrm{1}}(x)为主分布，\alpha为污染系数，它衡量数据受污染的程度。1964年Huber^[7]考虑了F_{\mathrm{1}}(x)为标准正态分布而F_{\mathrm{2}}(x)是一族关于原点对称的分布的情形。一个特殊的情况就是:F_{\mathrm{1}}(x)和F_{\mathrm{2}}(x)是具有相同形式但含有不同参数的分布，即F_{\mathrm{1}}(x)=F(x;{\theta }_{\mathrm{1}}),F_{\mathrm{2}}(x)=F(x;{\theta }_{\mathrm{2}})^[8]。近十年来，删失数据和截尾数据^[9]的发展已经取得了很大的成果，但是关于污染数据^[10]的研究却很少。污染数据的处理越来越受到重视，有广泛应用背景。

污染数据服从混合分布，对于混合分布常用的参数估计方法有矩估计、极大似然估计等方法。但遗憾的是，在小样本场合，这些估计将会有较大的偏差和方差。小样本估计问题在社会经济系统、生物医学系统和航空航海系统等领域中是普遍存在的。在疲劳寿命实验中，随着科学技术的迅速发展，产品的可靠性越来越高。很多时候由于人力、财力等的限制，重复实验的代价昂贵，这时可以获取的实验数据是非常少的。如果依然采用极大似然估计等方法进行统计推断，一般无法达到期望的效果。因此，适用于小样本情况下的疲劳寿命分布特性研究是很有必要的。早在20世纪60年代，国外学者已经开始研究小样本问题，80至90年代国内外普遍开始重视小样本问题的研究。比小样本情况更极端的一种情形为单样本情况，即只根据一个观测样本进行相关的统计推断。单样本情形下的统计推断研究非常少。(1)式中，在F_{\mathrm{1}}(x)和F_{\mathrm{2}}(x)是具有相同形式但含有不同参数的指数分布的特殊情形下，2017年Nie等^[11]考虑了基于混合指数分布的单个观测值对F_{\mathrm{1}}(x)的可靠性函数的无偏估计。他们假设随机变量T服从参数为\lambda （\lambda >\mathrm{0}）的指数分布，记为T\sim \mathcal{E}(\lambda )，则T的可靠性函数为R(t)=\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{t}{\lambda }),t\ge \mathrm{0}。两个指数分布的混合，记为p\mathcal{E}(\lambda )+(\mathrm{1}-p)\mathcal{E}(\theta \lambda ),p\in \mathbb{R},\lambda >\mathrm{0},\theta >\mathrm{0}。值得注意的是，有时p在[0,1]之间时无法满足实际情况的需要，因此这里的混合分布可以理解为广义的混合分布。在p,\theta已知，\lambda未知的情况下，Nie等基于混合指数分布p\mathcal{E}(\lambda )+(\mathrm{1}-p)\mathcal{E}(\theta \lambda )的单个观测值首先给出了R(t)的两个无偏估计，然后讨论了这两个无偏估计的适当性，最后给出了这两个无偏估计的方差的无偏估计。

寿命数据分析方法是进行可靠性工程研究的基础。常见系统中元件可靠性参数以及工程材料的疲劳寿命和强度分布都可以用Gamma分布很好地描述。指数分布、卡方分布、Erlang分布是Gamma分布的特例。自从Gamma分布被提出，各国研究者就对Gamma分布的参数估计问题产生了极大兴趣并为之付出了不懈努力。极大似然估计、矩估计、贝叶斯估计等为最常用的方法。由于Gamma分布在拟合数据方面比较灵活，适应性较强，因此Gamma分布被广泛地应用于众多统计应用中。常用的寿命概率分布描述寿命数据规律的能力都有限，而一个拟合能力更好的寿命概率分布模型可以让我们更精确的进行统计推断。混合Gamma分布的稠密性是指对于任何非负随机变量，我们总可以选择一个适当的混合Gamma分布把它拟合到任意需要的精度。因此，从理论上讲，选择混合Gamma分布拟合寿命数据总是合理的^[12]。

本文假设随机变量T服从参数为(\alpha ,\beta )的Gamma分布，记为T\sim \mathrm{G}\mathrm{a}(\alpha ,\beta )，其概率密度函数为

其中\alpha >\mathrm{0}是形状参数，\beta >\mathrm{0}是尺度参数，\mathrm{\Gamma }(\cdot )为Gamma函数。Gamma分布与很多分布都有关系:当\alpha =\mathrm{1}时，它是指数分布;当\alpha为正整数时，它是Erlang分布;当\beta =\frac{\mathrm{1}}{\mathrm{2}}时，它是卡方分布。T的可靠性函数为R(t)={\int }_t^{\mathrm{\infty }}f(x;\alpha ,\beta )\mathrm{d}x,t\ge \mathrm{0}。考虑两个Gamma分布的混合p\mathrm{G}\mathrm{a}(\alpha ,\beta )+(\mathrm{1}-p)\mathrm{G}\mathrm{a}(\alpha ,{\theta }^{-\mathrm{1}}\beta )，其中p\in \mathbb{R},\theta >\mathrm{0},\alpha >\mathrm{0},\beta >\mathrm{0}。p(\mathrm{1}-p)=\mathrm{0}或\theta =\mathrm{1}时混合Gamma分布会退化成单个Gamma分布，故本文假定p(\mathrm{1}-p)\ne \mathrm{0}且\theta \ne \mathrm{1}。我们认为数据本身服从的分布为\mathrm{G}\mathrm{a}(\alpha ,\beta )，但受到了来自分布\mathrm{G}\mathrm{a}(\alpha ,{\theta }^{-\mathrm{1}}\beta )数据的影响，我们更想得到分布\mathrm{G}\mathrm{a}(\alpha ,\beta )的可靠性函数R(t)。在p,\theta已知，\alpha ,\beta未知的情况下，本文基于混合Gamma分布的单个观测值首先给出R(t)的两个无偏估计，然后讨论了这两个无偏估计的适当性。假设p\in \mathbb{R},\theta >\mathrm{0},\alpha >\mathrm{0},\beta >\mathrm{0}。p\mathrm{G}\mathrm{a}(\alpha ,\beta )+(\mathrm{1}-p)\mathrm{G}\mathrm{a}(\alpha ,{\theta }^{-\mathrm{1}}\beta )的“密度”记为

概率密度函数应该具有非负性。由于此时p的范围并不仅限于[\mathrm{0,1}]之间，所以此时\psi (x;p,\theta ,\alpha ,\beta )可能为负；只有当\psi (x;p,\theta ,\alpha ,\beta )非负时，它才成为概率密度函数。

1　 \psi (x;p,\theta ,\alpha ,\beta )的非负性

下面根据p的取值范围分三种情况讨论\psi (x;p,\theta ,\alpha ,\beta )的非负性。

1)　 p>\mathrm{1}。\psi (x;p,\theta ,\alpha ,\beta )非负，即要求

p\frac{{\beta }^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\beta x}+(\mathrm{1}-p)\frac{(\beta \cdot {\theta }^{-\mathrm{1}}{)}^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\frac{\beta }{\theta }x}\ge \mathrm{0}

上式等价于

({\theta }^{-\mathrm{1}}-\mathrm{1})\beta x\ge \mathrm{l}\mathrm{n}\frac{p-\mathrm{1}}{p{\theta }^{\alpha }}

解上述不等式得：。因此，当p>\mathrm{1}时，\psi (x;p,\theta ,\alpha ,\beta )\ge \mathrm{0}当且仅当。

2)　 。\psi (x;p,\theta ,\alpha ,\beta )非负，即要求

p\frac{{\beta }^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\beta x}+(\mathrm{1}-p)\frac{(\beta \cdot {\theta }^{-\mathrm{1}}{)}^{\alpha }}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\frac{\beta }{\theta }x}\ge \mathrm{0}

上式等价于

({\theta }^{-\mathrm{1}}-\mathrm{1})\beta x\le \mathrm{l}\mathrm{n}\frac{p-\mathrm{1}}{p{\theta }^{\alpha }}

解上述不等式得：。因此，当时，\psi (x;p,\theta ,\alpha ,\beta )\ge \mathrm{0}当且仅当。

3)　 。显然\psi (x;p,\theta ,\alpha ,\beta )\ge \mathrm{0}。

综上所述，\psi (x;p,\theta ,\alpha ,\beta )成为概率密度函数的充要条件为下面三种情况之一成立:

(d1)　 p>\mathrm{1}且;

(d2)　 且;

(d3) 。

2　 可靠性函数的两个无偏估计

假设p\in \mathbb{R},\theta >\mathrm{0}已知，\alpha >\mathrm{0},\beta >\mathrm{0}未知。设T\sim \mathrm{G}\mathrm{a}(\alpha ,\beta )，T的概率密度函数为f(x;\alpha ,\beta )，可靠性函数为R(t)={\int }_t^{\mathrm{\infty }}f(x;\alpha ,\beta )\mathrm{d}x, t\ge \mathrm{0}。易知R(\mathrm{0})=\mathrm{1}。若条件(d1)、(d2)、(d3)之一成立，则\psi (x;p,\theta ,\alpha ,\beta )是混合Gamma分布p\mathrm{G}\mathrm{a}(\alpha ,\beta )+(\mathrm{1}-p)\mathrm{G}\mathrm{a}(\alpha ,{\theta }^{-\mathrm{1}}\beta )的概率密度函数。假设X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本，我们的目的是基于单个观测X找到R(t)(t>\mathrm{0})的无偏估计。

假设g(x;t)是定义在(\mathrm{0},\mathrm{\infty })\times (\mathrm{0},\mathrm{\infty })上的实值函数且关于x左连续，下面的定理给出g(X;t)是R(t)的无偏估计的必要条件。

定理1 假设条件(d1)、(d2)、(d3)之一成立，X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本。若E(g(X;t))=R(t),\forall t>\mathrm{0}，则g(x;t)满足

I(x>t)=pg(x;t)+(\mathrm{1}-p)g(\theta x;t),\hspace{1em}\forall t,x>\mathrm{0}(2)

证 若E(g(X;t))=R(t),\forall t>\mathrm{0}，则

{\int }_{\mathrm{0}}^{\mathrm{\infty }}I(x>t)f(x;\alpha ,\beta )\mathrm{d}x=R(t)=E(g(X;t))=

\begin{array}{l}{\int }_{\mathrm{0}}^{\mathrm{\infty }}g(x;t)(pf(x;\alpha ,\beta )+(\mathrm{1}-p)f(x;\alpha ,{\theta }^{-\mathrm{1}}\beta ))\mathrm{d}x=\\ p{\int }_{\mathrm{0}}^{\mathrm{\infty }}g(x;t)f(x;\alpha ,\beta )\mathrm{d}x+(\mathrm{1}-p){\int }_{\mathrm{0}}^{\mathrm{\infty }}g(x;t)f(x;\alpha ,{\theta }^{-\mathrm{1}}\beta )\mathrm{d}x=\\ p{\int }_{\mathrm{0}}^{\mathrm{\infty }}g(x;t)f(x;\alpha ,\beta )\mathrm{d}x+(\mathrm{1}-p){\int }_{\mathrm{0}}^{\mathrm{\infty }}g(\theta x;t)f(x;\alpha ,\beta )\mathrm{d}x=\\ {\int }_{\mathrm{0}}^{\mathrm{\infty }}(pg(x;t)+(\mathrm{1}-p)g(\theta x;t))f(x;\alpha ,\beta )\mathrm{d}x\end{array}

从而\forall t>\mathrm{0}，有

{\int }_{\mathrm{0}}^{\mathrm{\infty }}(I(x>t)-(pg(x;t)+(\mathrm{1}-p)g(\theta x;t)))\frac{\mathrm{1}}{\mathrm{\Gamma }(\alpha )}{x}^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\beta x}\mathrm{d}x=\mathrm{0},\forall \beta >\mathrm{0}

由拉氏变换的唯一性，\forall t>\mathrm{0}，有

I(x>t)-(pg(x;t)+(\mathrm{1}-p)g(\theta x;t))=\mathrm{0},\mathrm{a}. \mathrm{e}. x>\mathrm{0}

又I(x>t)和g(x;t)是关于x左连续的，所以

I(x>t)-(pg(x;t)+(\mathrm{1}-p)g(\theta x;t))=\mathrm{0}，\forall t,x>\mathrm{0}

证毕。

引理1 假设条件(d1)、(d2)、(d3)之一成立，X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本。若g(x;t)满足(2)式，则有

E(g(\theta X;t))=R({\theta }^{-\mathrm{1}}t),\hspace{1em}\forall t>\mathrm{0}(3)

证 若g(x;t)满足(2)式，则有

\begin{array}{l}E(g(\theta X;t))={\int }_{\mathrm{0}}^{\mathrm{\infty }}g(\theta x;t)(pf(x;\alpha ,\beta )+(\mathrm{1}-p)f(x;\alpha ,{\theta }^{-\mathrm{1}}\beta ))\mathrm{d}x=\\ p{\int }_{\mathrm{0}}^{\mathrm{\infty }}g(\theta x;t)f(x;\alpha ,\beta )\mathrm{d}x+(\mathrm{1}-p){\int }_{\mathrm{0}}^{\mathrm{\infty }}g({\theta }^{\mathrm{2}}x;t)f(x;\alpha ,\beta )\mathrm{d}x=\\ {\int }_{\mathrm{0}}^{\mathrm{\infty }}(pg(\theta x;t)+(\mathrm{1}-p)g({\theta }^{\mathrm{2}}x;t))f(x;\alpha ,\beta )\mathrm{d}x=\\ {\int }_{\mathrm{0}}^{\mathrm{\infty }}I(\theta x>t)f(x;\alpha ,\beta )\mathrm{d}x={\int }_{\mathrm{0}}^{\mathrm{\infty }}I(x>{\theta }^{-\mathrm{1}}t)f(x;\alpha ,\beta )\mathrm{d}x=R({\theta }^{-\mathrm{1}}t)\end{array}

证毕。

形式上定义

g_{\mathrm{1}}(x;t)\triangleq \sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{k}I（x>{\theta }^{-k}t）,x,t>\mathrm{0}(4)

g_{\mathrm{2}}(x;t)\triangleq \sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{k}I（x>{\theta }^{k+\mathrm{1}}t）,x,t>\mathrm{0}(5)

下面的引理给出g_{\mathrm{1}}(x;t)和g_{\mathrm{2}}(x;t)收敛的充要条件。

引理2p>\frac{\mathrm{1}}{\mathrm{2}}或是g_{\mathrm{1}}(x;t)收敛的充要条件;或\theta >\mathrm{1}是g_{\mathrm{2}}(x;t)收敛的充要条件。

证 首先分情况讨论g_{\mathrm{1}}(x;t)的收敛性。

1) p>\frac{\mathrm{1}}{\mathrm{2}}。由于，因此|g_{\mathrm{1}}(x;t)|\le \sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left|\mathrm{1}-\frac{\mathrm{1}}{p}\right|}^k=\frac{\mathrm{1}}{p-|p-\mathrm{1}|}，此时g_{\mathrm{1}}(x;t)收敛。

2) 。若，则g_{\mathrm{1}}(x;t)=\mathrm{0};若，令N=⌈\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{t}{x}⌉-\mathrm{1}，则，故g_{\mathrm{1}}(x;t)=\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^k收敛。

3) p\le \frac{\mathrm{1}}{\mathrm{2}}且\theta >\mathrm{1}。若，则g_{\mathrm{1}}(x;t)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^k发散;若，N=⌊\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{t}{x}⌋，则，故g_{\mathrm{1}}(x;t)=\sum _{k=N+\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^k发散。

综上可知，p>\frac{\mathrm{1}}{\mathrm{2}}或是g_{\mathrm{1}}(x;t)收敛的充要条件。接下来分情况讨论g_{\mathrm{2}}(x;t)的收敛性。

1) 。由于，因此|g_{\mathrm{2}}(x;t)|\le \sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left|\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right|}^k=\frac{\mathrm{1}}{\mathrm{1}-p-|p|}，此时g_{\mathrm{2}}(x;t) 收敛。

2) \theta >\mathrm{1}。若，则g_{\mathrm{2}}(x;t)=\mathrm{0};若，令N=⌈\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{x}{t}⌉-\mathrm{1}，则，故g_{\mathrm{2}}(x;t)=\sum _{k=\mathrm{0}}^{N-\mathrm{1}}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^k收敛。

3) p\ge \frac{\mathrm{1}}{\mathrm{2}}且。若，则g_{\mathrm{2}}(x;t)=\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^k发散;若，令N=⌊\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{x}{t}⌋，则，故g_{\mathrm{2}}(x;t)=\sum _{k=N}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^k发散。

综上可知，或\theta >\mathrm{1}是g_{\mathrm{2}}(x;t)收敛的充要条件。证毕。

g_{\mathrm{1}}(x;t)和g_{\mathrm{2}}(x;t)的收敛性情况如表1。

假设条件(d1)、(d2)、(d3)之一成立。接下来我们在表1所示条件下证明g_{\mathrm{1}}(X;t)和g_{\mathrm{2}}(X;t)是R(t)的两个无偏估计量，其中X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本。

定理2 若(d1)或成立，则有E(g_{\mathrm{1}}(X;t))=R(t),\forall t>\mathrm{0}; 若(d2)或成立，则有E(g_{\mathrm{2}}(X;t))=R(t),\forall t>\mathrm{0}。

证 假设(d1)、(d2)、、之一成立，X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本，则

E(I(X>t))={\int }_t^{\mathrm{\infty }}pf(x;\alpha ,\beta )+(\mathrm{1}-p)f(x;\alpha ,{\theta }^{-\mathrm{1}}\beta )\mathrm{d}x=

p{\int }_t^{\mathrm{\infty }}f(x;\alpha ,\beta )\mathrm{d}x+(\mathrm{1}-p){\int }_{{\theta }^{-\mathrm{1}}t}^{\mathrm{\infty }}f(x;\alpha ,\beta )\mathrm{d}x=pR(t)+(\mathrm{1}-p)R({\theta }^{-\mathrm{1}}t), \forall t>\mathrm{0}(6)

因此，

R(t)=E\left(\frac{\mathrm{1}}{p}I(X>t)\right)-\frac{\mathrm{1}-p}{p}R({\theta }^{-\mathrm{1}}t)=

\begin{array}{l}E\left(\frac{\mathrm{1}}{p}I(X>t)\right)-\frac{\mathrm{1}-p}{p}\left(E\left(\frac{\mathrm{1}}{p}I(X>{\theta }^{-\mathrm{1}}t)\right)-\frac{\mathrm{1}-p}{p}R({\theta }^{-\mathrm{2}}t)\right)=\\ E\left(\frac{\mathrm{1}}{p}I(X>t)+\frac{\mathrm{1}}{p}\frac{p-\mathrm{1}}{p}I(X>{\theta }^{-\mathrm{1}}t)\right)+{\left(\frac{p-\mathrm{1}}{p}\right)}^{\mathrm{2}}R({\theta }^{-\mathrm{2}}t)=\\ E\left(\frac{\mathrm{1}}{p}I(X>t)+\frac{\mathrm{1}}{p}\frac{p-\mathrm{1}}{p}I(X>{\theta }^{-\mathrm{1}}t)\right)+{\left(\frac{p-\mathrm{1}}{p}\right)}^{\mathrm{2}}\left(E\left(\frac{\mathrm{1}}{p}I(X>{\theta }^{-\mathrm{2}}t)\right)-\frac{\mathrm{1}-p}{p}R({\theta }^{-\mathrm{3}}t)\right)=\\ E\left(\frac{\mathrm{1}}{p}I(X>t)+\frac{\mathrm{1}}{p}\frac{p-\mathrm{1}}{p}I(X>{\theta }^{-\mathrm{1}}t)+\frac{\mathrm{1}}{p}{\left(\frac{p-\mathrm{1}}{p}\right)}^{\mathrm{2}}I(X>{\theta }^{-\mathrm{2}}t)\right)+{\left(\frac{p-\mathrm{1}}{p}\right)}^{\mathrm{3}}R({\theta }^{-\mathrm{3}}t)=\\ \cdots =\\ E\left(\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{k}I(X>{\theta }^{-k}t)\right)+{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{N+\mathrm{1}}R({\theta }^{-(N+\mathrm{1})}t),N=\mathrm{0,1},\mathrm{2},\cdots \end{array}

当(d1)或成立时，\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{N+\mathrm{1}}R({\theta }^{-(N+\mathrm{1})}t)=\mathrm{0}，由于

\left|\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{k}I(x>{\theta }^{-k}t)\right|\le \frac{\mathrm{1}}{p-|p-\mathrm{1}|}

根据控制收敛定理有

R(t)=\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}E\left(\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{k}I(X>{\theta }^{-k}t)\right)=E\left(\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{k}I(X>{\theta }^{-k}t)\right)=E(g_{\mathrm{1}}(X;t))

所以当(d1)或成立时，g_{\mathrm{1}}(X;t)是R(t)的一个无偏估计。

在(6)式中用\theta t代替t，得到：E(I(X>\theta t))=pR(\theta t)+(\mathrm{1}-p)R(t),\forall t>\mathrm{0}，因此

\begin{array}{l}R(t)=E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>\theta t)\right)-\frac{p}{\mathrm{1}-p}R(\theta t)=\\ E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>\theta t)\right)-\frac{p}{\mathrm{1}-p}\left(E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>{\theta }^{\mathrm{2}}t)\right)-\frac{p}{\mathrm{1}-p}R({\theta }^{\mathrm{2}}t)\right)=\\ E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>\theta t)-\frac{\mathrm{1}}{\mathrm{1}-p}\frac{p}{\mathrm{1}-p}I(X>{\theta }^{\mathrm{2}}t)\right)+{\left(-\frac{p}{\mathrm{1}-p}\right)}^{\mathrm{2}}R({\theta }^{-\mathrm{2}}t)=\\ E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>\theta t)-\frac{\mathrm{1}}{\mathrm{1}-p}\frac{p}{\mathrm{1}-p}I(X>{\theta }^{\mathrm{2}}t)\right)+{\left(-\frac{p}{\mathrm{1}-p}\right)}^{\mathrm{2}}\left(E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>{\theta }^{\mathrm{3}}t)\right)-\frac{p}{\mathrm{1}-p}R({\theta }^{\mathrm{3}}t)\right)=\\ E\left(\frac{\mathrm{1}}{\mathrm{1}-p}I(X>\theta t)-\frac{\mathrm{1}}{\mathrm{1}-p}\frac{p}{\mathrm{1}-p}I(X>{\theta }^{\mathrm{2}}t)+\frac{\mathrm{1}}{\mathrm{1}-p}{\left(-\frac{p}{\mathrm{1}-p}\right)}^{\mathrm{2}}I(X>{\theta }^{\mathrm{3}}t)\right)+{\left(-\frac{p}{\mathrm{1}-p}\right)}^{\mathrm{3}}R({\theta }^{\mathrm{3}}t)=\\ \cdots =\\ E\left(\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{\mathrm{1}-p}{\left(-\frac{p}{\mathrm{1}-p}\right)}^{k}I(X>{\theta }^{k+\mathrm{1}}t)\right)+{\left(-\frac{p}{\mathrm{1}-p}\right)}^{N+\mathrm{1}}R({\theta }^{N+\mathrm{1}}t)=\\ E\left(\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{k}I(X>{\theta }^{k+\mathrm{1}}t)\right)+{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{N+\mathrm{1}}R({\theta }^{N+\mathrm{1}}t),N=\mathrm{0,1},\mathrm{2},\cdots \end{array}

当(d2)或成立时，\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p})}^{N+\mathrm{1}}R({\theta }^{N+\mathrm{1}}t)=\mathrm{0}，由于

|\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{k}I(x>{\theta }^{k+\mathrm{1}}t)|\le \frac{\mathrm{1}}{\mathrm{1}-p-|p|}

根据控制收敛定理，有

\begin{array}{l}R(t)=\underset{N\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}E(\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{k}I(X>{\theta }^{k+\mathrm{1}}t))=\\ E(\sum _{k=\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^{k}I(X>{\theta }^{k+\mathrm{1}}t))=E(g_{\mathrm{2}}(X;t))\end{array}

所以当(d2)或成立时，g_{\mathrm{2}}(X;t)是R(t)的一个无偏估计。证毕。

g_{\mathrm{1}}(X;t)和g_{\mathrm{2}}(X;t)都是R(t)的基于单个样本X的无偏估计量，与Nie等^[11]得到的两个无偏估计形式上一致。他们工作中的X服从混合指数分布，这里的X服从混合Gamma分布。

3　 无偏估计的适当性

假设条件(d1)、(d2)、(d3)之一成立，X是来自\psi (x;p,\theta ,\alpha ,\beta )的单个样本，在p\ne \frac{\mathrm{1}}{\mathrm{2}}时g_{\mathrm{1}}(X;t)或g_{\mathrm{2}}(X;t)是可靠性函数R(t)的无偏估计。由于R(t)取值范围为[0,1]，因此要保证g_{\mathrm{1}}(X;t),g_{\mathrm{2}}(X;t)的取值范围也为[0,1]，称此时的估计是适当的。接下来讨论在什么样的条件下g_{\mathrm{1}}(X;t),g_{\mathrm{2}}(X;t)才是适当的。

定理3 若(d1)成立，则g_{\mathrm{1}}(X;t)是R(t)的一个适当的无偏估计;若(d2)成立，则g_{\mathrm{2}}(X;t)是R(t)的一个适当的无偏估计。

证 若(d1)成立，根据表1，此时g_{\mathrm{1}}(x;t)收敛。由于，因此g_{\mathrm{1}}(x;t)满足：。从而g_{\mathrm{1}}(X;t)是R(t)适当的无偏估计。

若(d2)成立，根据表1，此时g_{\mathrm{2}}(x;t)收敛。由于，因此g_{\mathrm{2}}(x;t)满足：。从而g_{\mathrm{2}}(X;t)是R(t)适当的无偏估计。证毕。

注1 由表1可知，若(d1)成立，g_{\mathrm{2}}(x;t)不是收敛的。若(d2)成立，g_{\mathrm{1}}(x;t)不是收敛的。

注2 当时，g_{\mathrm{1}}(X;t)是R(t)的无偏估计，但不是适当的;当时g_{\mathrm{2}}(X;t)是R(t)的无偏估计，但不是适当的。下面说明之。

1) 对于g_{\mathrm{1}}(x;t)，根据定理2，当时g_{\mathrm{1}}(X;t)是R(t)的一个无偏估计。根据\theta的范围分两种情况讨论g_{\mathrm{1}}(X;t)的适当性。

a) 。若，令N=⌈\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{t}{x}⌉-\mathrm{1}，则，故

g_{\mathrm{1}}(x;t)=\sum _{k=\mathrm{0}}^N\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^k=\mathrm{1}-{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{N+\mathrm{1}}

由于，因此当N为偶数时g_{\mathrm{1}}(x;t)大于1。所以g_{\mathrm{1}}(X;t)不是适当的。

b) 。若，令N=⌊\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{t}{x}⌋，则，故

g_{\mathrm{1}}(x;t)=\sum _{k=N+\mathrm{1}}^{\mathrm{\infty }}\frac{\mathrm{1}}{p}{\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^k={\left(\mathrm{1}-\frac{\mathrm{1}}{p}\right)}^{N+\mathrm{1}}

由于，因此当N为偶数时g_{\mathrm{1}}(x;t)小于0。所以g_{\mathrm{1}}(X;t)不是适当的。

2) 对于g_{\mathrm{2}}(x;t)，根据定理2，当时g_{\mathrm{2}}(X;t)是R(t)的一个无偏估计。根据\theta的范围分两种情况讨论g_{\mathrm{2}}(X;t)的适当性。

a) 。若，令N=⌊\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{x}{t}⌋，则，故

g_{\mathrm{2}}(x;t)=\sum _{k=N}^{\mathrm{\infty }}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^k={\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^N

由于，因此当N为奇数时g_{\mathrm{2}}(x;t)小于0。所以g_{\mathrm{2}}(X;t)不是适当的。

b) 。若，令N=⌈\mathrm{l}\mathrm{o}{\mathrm{g}}_{\theta }\frac{x}{t}⌉-\mathrm{1}，则，故

g_{\mathrm{2}}(x;t)=\sum _{k=\mathrm{0}}^{N-\mathrm{1}}\frac{\mathrm{1}}{\mathrm{1}-p}{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^k=\mathrm{1}-{\left(\mathrm{1}-\frac{\mathrm{1}}{\mathrm{1}-p}\right)}^N

由于，因此当N为奇数时g_{\mathrm{2}}(x;t)大于1。所以g_{\mathrm{2}}(X;t)不是适当的。

根据上面的讨论，g_{\mathrm{1}}(X;t)和g_{\mathrm{2}}(X;t)不可能同时是R(t)的适当的无偏估计。可以发现，只有当混合权重(p或\mathrm{1}-p)其中一个为负值的时候，g_{\mathrm{1}}(X;t)或g_{\mathrm{2}}(X;t)才可能是R(t)的适当的无偏估计。这个结论和Nie等^[11]指数分布情形下的结论一致。

4　 数值模拟

假设条件(d1)成立，T~\mathrm{G}\mathrm{a}(\alpha ,\beta )，其概率密度函数为f(x;\alpha ,\beta )，其可靠性函数为R(t)={\int }_t^{\mathrm{\infty }}f(x;\alpha ,\beta )\mathrm{d}x。混合Gamma分布记为p\mathrm{G}\mathrm{a}(\alpha ,\beta )+(\mathrm{1}-p)\mathrm{G}\mathrm{a}(\alpha ,{\theta }^{-\mathrm{1}}\beta )，其概率密度函数为\psi (x;p,\theta ,\alpha ,\beta )，其分布函数为F(x)={\int }_{\mathrm{0}}^x\psi (x;p,\theta ,\alpha ,\beta )\mathrm{d}x。用逆采样方法生成混合Gamma分布的随机观测:首先，假设随机变量U服从区间[0,1]上的均匀分布，利用计算机生成U的一个随机观测u，然后求解x使得F(x)=u。取x为混合Gamma分布的一个随机观测。根据x得到R(t)(t>\mathrm{0})的无偏估计g_{\mathrm{1}}(x;t)，并与真实的可靠性函数R(t)进行比较。

取p=\mathrm{1.1},\theta =\mathrm{0.8},\alpha =\mathrm{4},\beta =\mathrm{1}，当用计算机生成单个随机观测x后，根据(4)式求出g_{\mathrm{1}}(x;t)，并与真实的可靠性函数R(t)作比较。当生成不同的单个随机观测时，g_{\mathrm{1}}(x;t)和R(t)的比较如图1和图2。

图1为x=\mathrm{3.2}时的模拟结果，图2为x=\mathrm{5.4}时的模拟结果。可以看出，图1和图2中g_{\mathrm{1}}(x;t)能够较好地估计R(t)。g_{\mathrm{1}}(x;t)的估计效果取决于混合Gamma分布的随机观测x。若在单个观测的情况下利用经验分布函数去估计R(t)，只能得到只有一个阶梯的阶梯函数(取值为0或1)。由(4)式知g_{\mathrm{1}}(x;t)是一个阶梯函数，且有好几个阶梯，因此在单个观测的情况下，g_{\mathrm{1}}(x;t)能够比经验分布函数更好地估计R(t)。所以，本文的工作是有意义的。

接下来验证g_{\mathrm{1}}(x;t)的无偏性。每进行一次试验，得到一个观测，进而可得一个无偏估计g_{\mathrm{1}}(x;t)。 当进行大量试验时，希望由此得到的大量g_{\mathrm{1}}(x;t)的平均值能精确地估计R(t)。取p=\mathrm{1.1},\theta =\mathrm{0.8}，重复进行500次随机试验，一共生成500个观测x_i,i=\mathrm{1,2},\cdots ,\mathrm{500}，其中x_i为第i次试验得到的观测。计算R(t)的500个无偏估计g_{\mathrm{1}}(x_i;t),i=\mathrm{1,2},\cdots ,\mathrm{500}，其平均值记为

g_{\mathrm{1}}(t)=\frac{\mathrm{1}}{\mathrm{500}}\sum _{i=\mathrm{1}}^{\mathrm{500}}{g}_{\mathrm{1}}(x_i;t)

当\alpha ,\beta取不同的值时，g_{\mathrm{1}}(t)和R(t)的比较如图3和图4。

图3为\alpha =\mathrm{4},\beta =\mathrm{1}时的模拟结果，图4为\alpha =\mathrm{3},\beta =\mathrm{2}时的模拟结果。图3与图1和2均是在p=\mathrm{1.1},\theta =\mathrm{0.8},\alpha =\mathrm{4},\beta =\mathrm{1} 时的模拟结果，不同的是图1和图2是根据单个无偏估计的模拟，而图3是根据500个无偏估计的模拟。通过图3与图1和2的比较可以看出:当用到的无偏估计的数量很多时，多个估计的平均值能比较精确的拟合R(t)，即验证了g_{\mathrm{1}}(x;t)的无偏性。图4是\alpha =\mathrm{3},\beta =\mathrm{2}时利用500个无偏估计的模拟结果，进一步验证了g_{\mathrm{1}}(x;t)的无偏性。

5　 结 语

本文基于混合Gamma分布的单个观测值讨论了Gamma分布可靠性函数的无偏估计及其相关性质。 首先，给出了混合Gamma分布的概率密度函数可定义的充要条件。其次，在p\in \mathbb{R},\theta >\mathrm{0}已知，\alpha >\mathrm{0},\beta >\mathrm{0}未知的情况下，基于混合Gamma分布的单个观测值给出了Gamma分布的可靠性函数R(t)的无偏估计g_{\mathrm{1}}(x;t)和g_{\mathrm{2}}(x;t)。接下来，我们讨论了g_{\mathrm{1}}(x;t)和g_{\mathrm{2}}(x;t)的适当性，发现当p>\mathrm{1}且时g_{\mathrm{1}}(x;t)在[0,1]之间，当且时g_{\mathrm{2}}(x;t)在[0,1]之间。最后，我们利用计算机生成模拟数据，验证了g_{\mathrm{1}}(x;t) 和g_{\mathrm{2}}(x;t)的无偏性。

Nie等^[11]在混合指数分布的情况给出了指数分布可靠性函数的无偏估计及其相关性质，本文中的相关结论是基于混合Gamma分布得到的，是其工作的推广。

参考文献

[1] NEWCOMB S. A generalized theory of the combination of observations so as to obtain the best result[J]. American Journal of Mathematics, 1886, 8(4):343-366.

[2] PEARSON K. Contributions to the mathematical theory of evolution [J]. Philosophical Transactions of the Royal Society of London, A, 1894, 185:71-110.

[3] WILLMOT G E, LIN X S. Risk modelling with the mixed Erlang distribution [J]. Applied Stochastic Models in Business & Industry, 2015, 27(1):2-16.

[4] 王建康, 盖钧镒. 混合分布理论及应用[J]. 生物数学学报, 1995, 10(3):87-92.

[5] WANG J K, GAI J Y. Mixed distribution theory and its application [J]. Journal of Biomathematics, 1995, 10(3):87-92 (Ch).

[6] 朱利平, 卢一强, 茆诗松. 混合指数分布的参数估计[J]. 应用概率统计, 2006, 22(2)：137-150.

[7] ZHU L P, LU Y Q, MAO S S. Parameter estimation of mixed exponential distribution[J]. Journal of Applied Probability and Statistics, 2006, 22(2)：137-150 (Ch).

[8] DAVIS D J F. An analysis of some failure data [J]. Journal of the American Statistical Association, 1952, 47(258): 113-150.

[9] HUBER P J. Robust estimation of a location parameter [J]. The Annals of Mathematical Statistics, 1964, 35(1):73-101.

[10] FROM S G. Optimal spacing of quantiles for the estimation of the mixing parameters in a mixture of two exponential distributions[J]. Communications in Statistics, 1989, 18(6):2201-2223.

[11] YAN J. Survival analysis: Techniques for censored and truncated data [J]. Biometrics, 2010, 99(467):900-901.

[12] 郑祖康, 吴雪明, 饶刚. 污染数据的处理[J]. 应用概率统计, 1998, 14(3):307-312.

[13] ZHENG Z K, WU X M, RAO G. Processing of contaminated data [J]. Journal of Applied Probability and Statistics, 1998, 14(3):307-312 (Ch).

[14] NIE K Y, SINHA B K, HEDAYAT A S. Unbiased estimation of reliability function from a mixture of two exponential distributions based on a single observation [J]. Statistics & Probability Letters, 2017, 127:7-13.

[15] 黄卓, 李苏军, 郭波. 基于混合Gamma分布的通用可靠性寿命数据拟合方法[J]. 航空学报, 2008, 29(2)：379-386.

[16] HUANG Z, LI S J, GUO B. Fitting method for general reliability life data based on mixed Gamma distribution [J]. Journal of Aeronautics, 2008, 29(2)：379-386 (Ch).

基金资助

国家自然科学基金(11171263)