基于二阶矩约束的最大熵先验信息解及其在数控系统可靠性评估中的应用

凌 光 戴 怡 李 曦

1.天津工程师范学院,天津,300222 2.华中科技大学,武汉,430074

0 引言

Bayes方法是解决小子样问题的一个重要方法,其核心和关键是如何利用先验信息确定先验分布。根据早期研究,数控系统失效时间服从双参数Weibull分布。目前,双参数先验分布的确定未能得到很好解决：基于1/(θβ)形式的假设的方法[1]虽然简化了双参数相互间的关系,但在一定程度上影响了该方法的精度;独立参数分布的截尾Γ分布和截尾Γ-1分布方法[2]计算复杂、不易操作,将对原分布的统计推断的难度转移到了参数分布计算上。

最大熵方法因充分利用了各种先验信息并尽量避免了引入其他不确定因素而受到重视,许多学者为之作出不懈努力。Mazzuchi等[3]将非参数的 Dirichlet函数与最大熵相结合,提出了MED先验分布;Kazama等[4]将最大熵方法应用到不等式约束规划问题中。Agrawal等[5]提出一种基于矩约束的最大熵双参数联合分布的确定方法,但该研究局限于一阶矩约束。目前,关于一般形式下双参数的最大熵先验信息解的研究还鲜见报道。本文基于双参数二阶矩约束,求取了最大熵先验信息解析解,并将之应用于数控系统可靠性评估中。

1 利用最大熵方法确定联合先验分布

如上所述,数控系统失效时间服从双参数Weibull分布 ,即

f(t;η,m)=m exp(-(t/η)m)/[η(t/η)m-1] t ＞0

它有η、m两个参数,其中,η为特征寿命(分布的0.632分位数);m为形状参数,对密度函数形状有很大影响。

根据二元联合熵的定义和最大熵原理,确定双参数联合先验分布 π(η,m),满足

参数满足的约束条件如下：

式中,ηi、mj分别为两参数η、m的 i阶和 j阶原点矩;k、n分别为两参数原点矩的最高阶数,可以相同,也可以不同。

对联合密度函数π(η,m)的求解是个泛函意义下的条件极值问题。构造如下辅助泛函：

式中,λ0、λ1、…、λk、c1、…、cn为拉格朗日乘子。

本文主要讨论二阶矩约束下双参数联合先验分布 ,对于一般的k 、n,系数λ0、λ1、…、λk、c1、…、cn可通过数值计算[6]求得。故在二阶矩约束下,即当k=n=2时,先验函数形式为

2 双参数二阶矩先验分布的确定和稳健性分析

利用最大熵方法求解双参数联合先验分布过程中,需要参数样本的二阶矩信息。由于试验所得数据为数控系统失效时间,故本文利用自助法,通过构造再生子样来获取两参数η、m的二阶矩估计。详细过程可参见文献[7]。

在应用先验分布前,必须分析其对后验统计推断结果的影响,即分析先验分布的稳健性。一种简单常用的先验分布稳健性判断方法[8-9]就是利用相对似然边际分布函数m*(t|π)。该相对似然函数边际分布值代表以π(η,m)为先验分布时现场样本出现的相对概率。当获得现场样本(t1,t2,…,tN)后,可计算m*(ti|π)(i=1,…,N)。如果所得值都不是非常小(一般不小于10-3),则说明先验分布符合现场数据条件,即可以认为先验分布π(η,m)是稳健的。

对双参数Weibull分布,联合先验分布 π(η,m)的相对似然边际函数为

3 最大熵先验信息解在数控系统可靠性评估中的应用

据Bayes原理,设有数控系统寿命数据T=(t1,t2,…,tN),则其似然函数可表示为

利用已知先验形式,则参数η、m联合后验分布可表示为

借助于数值计算,可以得到后验分布的点估计：

最后根据Weibull分布性质可知,数控系统平均后验无故障工作时间的估计为

设有数控系统失效数据30个,如表1所示,利用这30个历史数据,应用自助法生成200个容量为1000的自助样本组,即求得η和m的参数序列,进一步可以计算两参数的数学期望和方差：

表1 数控系统寿命数据 h

η=3811.6,m=1.4552,σ2η=7920.0,σ2m=0.0011

从而该Weibull分布尺度参数与位置参数的双参数联合先验密度函数为

设有现场数控系统的寿命数据20个,如表2所示。利用该数据对先验分布的稳健性进行验证,计算出其相对似然边际函数值亦如表2所示。从表2中数值可以看出,该先验分布对于所得的现场数据而言,稳健性良好,即利用该分布作为先验分布是可取的。将现场数据形成的似然函数乘以由历史数据求解得到的先验分布,得到由式(8)计算的后验分布,并可求出其点估计：

故而数控系统平均后验无故障工作时间的估计为

表2 现场数控系统寿命数据及其相对似然函数值

4 结语

最大熵方法能够最大限度地利用历史样本信息,适用于产品寿命长、实验数据难以获得的高可靠性产品的评估。本文提出的基于二阶矩约束、通过最大熵原理确定的先验分布,经检验具有良好的稳健性。算例表明,该方法适用于数控系统可靠性评估。

[1] Ferdous J,Uddin M U,Pandey M.Reliability Estimation with Weibull Inter Failure Times[J].Reliability Engineering and System Safety,1995,50：285-296.

[2] 张湘平.小子样统计推断与融合理论在武器系统评估中的应用研究[D].长沙：国防科学技术大学,2003.

[3] Mazzuchi T A,Soofi E S,Soyer R.Computation of Maximum Entropy Dirichlet for Modeling Lifetime Data[J].Computational Statistics&Data Analysis,2000,32：361-378.

[4] Kazama J,Tsujii J.Maximum Entropy Models with Inequality Constraints：a Case Study on Text Categorization[J].Machine Learning,2005,60：159-194.

[5] Agrawal D,Singh J K,Kumar A.Maximum Entropy-based Conditional Probability Distribution Runoff Model[J].Biosystems Engineering,2005,90(1)：103-113.

[6] 夏新涛,陈晓阳,张永振,等.航天轴承摩擦力矩的最大熵概率分布与bootstrap推断[J].宇航学报,2007,28(5)：1395-1400.

[7] 张金槐,刘琦,冯静.Bayes试验分析方法[M].长沙：国防科技大学出版社,2007.

[8] Berger JO.统计决策论及贝叶斯分析[M].贾乃光,译.北京：中国统计出版社,1998.

[9] 余礼军,Ben Heydecker.基于概率加权矩与最大熵原理的交通流统计分布估计方法[J].公路交通科技,2008,25(2)：113-133.