基于贝叶斯方法的逆高斯过程在退化分析中的应用

Application of inverse Gaussian processes based on Bayesian methods in degradation analysis

摘要

退化分析在产品可靠性评估和寿命预测方面起着至关重要的作用，特别是在失效时间难以直接观察的情况下。采用贝叶斯方法对逆高斯（IG）过程模型进行退化分析，基于贝叶斯框架的灵活性，将先验知识和不确定性纳入分析中，从而对产品寿命分布进行更翔实和更可靠的预测。此外，通过仿真模拟对先验分布和样本容量进行了全面的灵敏度分析，并使用OpenBUGS（一个用于贝叶斯分析的开源软件）执行计算。最后，通过一个经典的例子说明贝叶斯方法在逆高斯过程模型的退化分析中的适用性。

Abstract

Degradation analysis is crucial in product reliability assessment and lifetime prediction, especially when the time to failure is difficult to observe directly. In this paper, a Bayesian approach to degradation analysis of inverse Gaussian (IG) process models is used to provide a more informative and reliable prediction of product lifetime distributions based on the flexibility of the Bayesian framework, which incorporates prior knowledge and uncertainty into the analysis. Furthermore, a comprehensive sensitivity analysis of the prior distribution and sample size was performed through simulations, and the calculations were performed using OpenBUGS, an open-source software for Bayesian analysis. Finally, a classical example is used to illustrate the applicability of the Bayesian approach in the degradation analysis of inverse Gaussian process models.

Graphical abstract

关键词

逆高斯过程模型 / 贝叶斯方法 / 随机效应 / OpenBUGS

Key words

inverse Gaussian process model / Bayesian method / random effect / OpenBUGS

退化分析作为一种评估产品寿命和可靠性的先进方法，在传统寿命测试面临长期故障数据难以获得时提供了有效的替代。对于那些设计寿命长、很少发生故障的高可靠性产品，获取足够的直接故障数据以进行统计分析既费时又昂贵，因此，退化分析成了一种可行的替代方案。

考虑到退化建模的研究，一般情况下采用基于随机过程的模型，在这类模型中最常见的包括维纳过程和伽马过程。近年来，逆高斯过程模型提供了一种更加灵活和适应性更强的退化建模方法，能够更好地捕捉和描述退化数据中的复杂模式，极大提升了退化分析的适用性和预测的准确性。Wang等^[1]介绍了一种具有吸引力和灵活性的退化建模模型——逆高斯(IG)退化过程模型。Peng^[2]提出了一个逆正态伽马混合IG过程模型。随后，Peng等^[3]利用贝叶斯方法处理了稀疏退化和演化观测的退化分析。Guo等^[4]考虑了先验与退化增量之间的依赖关系，并提出了一种改进的逆高斯过程，引入退化增量补充估计寿命概率区间。Hao等^[5]研究了在不同随机效应下具有测量误差的逆高斯模型。Peng等^[6]假设观测量经过非线性变换后服从逆高斯分布，总结出一种逆高斯退化模型的改进方法。Sun等^[7]假定时间和测量误差具有相关性，分析了同时带有测量误差及随机影响的逆高斯退化过程模型。Fang等^[8]提出了一种带有关联随机效应的逆高斯退化模型。

传统的频率主义分析方法往往难以有效地整合来自先验知识的信息，而贝叶斯方法可以提供一个更为全面和灵活的分析框架，这种方法的一个关键特点是能够提供模型参数的后验分布。Ye等^[9]引入了一种基于随机漂移逆高斯模型退化分析的贝叶斯方法。Li等^[10]基于构建的贝叶斯框架迭代更新刀具磨损模型参数，提高了剩余使用寿命预测的准确性。Guo等^[11]利用贝叶斯理论对退化过程进行建模和可靠性评估，并研究不同先验分布对评估结果的影响。Zhao等^[12]将专家信息转化为先验分布,基于贝叶斯理论研究参数估计问题。

贝叶斯方法在退化分析中的应用增强了对复杂退化过程的建模能力，提高了寿命预测的准确性和可靠性。许晓东等^[13]在贝叶斯框架下推导剩余寿命预测结果的概率分布函数。王艺玮等^[14]基于贝叶斯理论估计退化模型参数。陈龙等^[15]采用贝叶斯更新逆高斯过程退化模型参数，更准确地预测桥梁剩余使用寿命。杨志远等^[16]分析随机参数对系统可靠性影响，提出基于贝叶斯理论的剩余寿命预测方法。

由此，本文在贝叶斯方法下对建立的4个逆高斯过程性能退化模型在给定的先验分布下进行参数估计及可靠性指标评定，结合实例，根据得到的参数估计结果选择最优模型。

1 逆高斯过程和OpenBUGS

1.1 逆高斯过程

由于在实践中，产品的退化路径在大多数情况下是单调的，因此退化过程可被视为IG过程。IG过程为具有独立、非负增量且遵循IG分布的随机过程。具体来说，具有函数\Lambda \left(t\right)和参数\mu ,\eta的IG过程\left\{Y\left(t\right),t\ge \mathrm{0}\right\}被定义为^[1,9]：

1) Y(\mathrm{0})\equiv \mathrm{0}；

2) 对于任意的t_{\mathrm{2}}>t_{\mathrm{1}}\ge s_{\mathrm{2}}>s_{\mathrm{1}}，有Y\left(t_{\mathrm{2}}\right)-Y\left(t_{\mathrm{1}}\right)和Y\left(s_{\mathrm{2}}\right)-Y\left(s_{\mathrm{1}}\right)是相互独立的退化增量；

3) 对于任意的t>s，退化增量Y\left(t\right)-Y\left(s\right)服从逆高斯分布\mathrm{I}\mathrm{G}\left(\mu \left(\Lambda \left(t\right)-\Lambda \left(s\right)\right),\eta {\left(\Lambda \left(t\right)-\Lambda \left(s\right)\right)}^{\mathrm{2}}\right)。其中：\Lambda \left(t\right)是单调递增函数，\Lambda \left(t\right)=t^q,q>\mathrm{0}。

逆高斯分布为y~\mathrm{I}\mathrm{G}\left(a,b\right)。a,b>\mathrm{0}的概率密度函数(PDF)和累积概率分布函数(CDF)分别为：

f\left(\left.y\right|a,b\right)=\sqrt[]{\frac{b}{\mathrm{2}\mathrm{\pi }{y}^{\mathrm{3}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{b{\left(y-a\right)}^{\mathrm{2}}}{\mathrm{2}{a}^{\mathrm{2}}y}\right],y>\mathrm{0}(1)

\begin{array}{l}F\left(\left.y\right|a,b\right)=\\ \Phi \left[\sqrt[]{\frac{b}{y}}\left(\frac{y}{a}-\mathrm{1}\right)\right]+\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}b}{a}\right)\Phi \left[-\sqrt[]{\frac{b}{y}}\left(\frac{y}{a}+\mathrm{1}\right)\right]\end{array}(2)

其中\Phi \left[\cdot \right]为标准正态分布的CDF。\mathrm{I}\mathrm{G}\left(a,b\right)的均值和方差分别为a和a^{\mathrm{3}}/b。

1.2 贝叶斯估计在OpenBUGS中的实现

OpenBUGS软件是一种使用马尔可夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)方法进行贝叶斯分析的强大工具^[17]。贝叶斯分析可以将先验知识与当前数据结合起来，当无法直接计算模型参数的后验分布时，可以使用MCMC方法来近似这些分布。利用MCMC方法，OpenBUGS软件将未知参数设定为随机变量，随后对概率模型进行求解，得到参数的轨迹图、均值、方差、置信区间以及收敛性诊断结果。这款软件操作灵活简洁，具体的操作步骤如图1所示。

OpenBUGS代表了一个进行贝叶斯分析的综合工具，其提供了一个平台，可以利用贝叶斯方法的优势对各种复杂性进行统计建模。通过MCMC方法的实现，为估计后验分布、适应先验知识的整合和统计推理中的不确定性的全面量化提供了一个稳健的框架。

2 逆高斯过程的贝叶斯分析

2.1 简单的IG过程的贝叶斯分析

2.1.1 模型描述

定义了一个单调递增函数\Lambda \left(t\right)及参数\mu和\eta的简单IG退化过程模型\left\{Y\left(t\right),t\ge \mathrm{0}\right\}。且Y\left(\mathrm{0}\right)\equiv \mathrm{0}，其中函数\Lambda \left(t\right)中的参数记为{\theta }_{\Lambda }。退化过程描述为Y\left(t\right)~\mathrm{I}\mathrm{G}\left(\mu \Lambda \left(t\right),\eta {\Lambda }^{\mathrm{2}}\left(t\right)\right),\mu ,\eta >\mathrm{0}，简单IG过程模型Y\left(t\right)的概率密度函数为：

\begin{array}{l}f\left(\left.y\right|\mu ,\eta ,\Lambda \left(t\right)\right)=\\ \sqrt[]{\frac{\eta {\Lambda }^{\mathrm{2}}\left(t\right)}{\mathrm{2}\mathrm{\pi }{y}^{\mathrm{3}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\eta {\left(y-\mu \Lambda \left(t\right)\right)}^{\mathrm{2}}}{\mathrm{2}{\mu }^{\mathrm{2}}y}\right],y>\mathrm{0}\end{array}(3)

在实际应用中，产品的退化量Y\left(t\right)第一次到达预定阈值D时失效，将首达时间认定为产品的寿命T，则产品寿命定义为：

T=\mathrm{i}\mathrm{n}\mathrm{f}\left\{\left.t\right|Y\left(t\right)\ge D\right\}(4)

产品寿命T_D的累积分布函数为：

\begin{array}{l}{F}_{T_D}\left(t\right)=P\left(Y\left(t\right)>D\right)\\ =\Phi \left[\sqrt[]{\frac{\eta }{D}}\left(\Lambda \left(t\right)-\frac{D}{\mu }\right)\right]-\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}\eta \Lambda \left(t\right)}{\mu }\right)\\ \Phi \left[-\sqrt[]{\frac{\eta }{D}}\left(\Lambda \left(t\right)+\frac{D}{\mu }\right)\right]\end{array}(5)

其中\Phi \left(\cdot \right)为标准正态分布函数。产品寿命T_D的概率密度函数为：

\begin{array}{l}{f}_{{\tau }_D}\left(t\right)=\sqrt[]{\frac{\eta }{D}}\Phi \left[\sqrt[]{\frac{\eta }{D}}\left(\Lambda \left(t\right)-\frac{D}{\mu }\right)\right]{\Lambda }^{\text{'}}\left(t\right)\\ -\frac{\mathrm{2}\eta }{\mu }{\Lambda }^{\text{'}}\left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}\eta \Lambda \left(t\right)}{\mu }\right)\Phi \left[-\sqrt[]{\frac{\eta }{D}}\left(\Lambda \left(t\right)+\frac{D}{\mu }\right)\right]\\ +\sqrt[]{\frac{\eta }{D}}{\Lambda }^{\text{'}}\left(t\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}\eta \Lambda \left(t\right)}{\mu }\right)\Phi \left[-\sqrt[]{\frac{\eta }{D}}\left(\Lambda \left(t\right)+\frac{D}{\mu }\right)\right]\end{array}(6)

产品的平均寿命可以反映产品可靠性，定义为：

E\left(t\right)={\int }_{\mathrm{0}}^{+\mathrm{\infty }}tf\left(t\right)\mathrm{d}t(7)

由产品寿命T_D的累积分布函数可以得到产品的可靠性函数为：

(8)

2.1.2 贝叶斯分析

假设对n个单元进行退化观察的观测值为Y\left(t\right)，设Y\left(t_{ij}\right)表示在t_{ij}时刻对单元i的第j次观测，其中j=\mathrm{1},\mathrm{2},\cdots ,m_i;i=\mathrm{1},\mathrm{2},\cdots ,n。设\mathrm{\Delta }{y}_{ij}=Y\left(t_{ij}\right)-Y\left(t_{i,j-\mathrm{1}}\right)为退化增量，在简单IG过程模型下退化增量\mathrm{\Delta }{y}_{ij}独立且\mathrm{\Delta }{y}_{ij}~\mathrm{I}\mathrm{G}\left(\mu \mathrm{\Delta }{\Lambda }_{ij},\eta \mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}\right)。在不考虑随机效应的情况下，令未知参数\mu ,\eta和{\theta }_{\Lambda }的先验分布为：

\mu ~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right),\eta ~\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right),{\theta }_{\Lambda }~U\left(a_{\Lambda },b_{\Lambda }\right)(9)

其中：\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right)为对数正态分布；\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right)为伽马分布；U\left(a_{\Lambda },b_{\Lambda }\right)为均匀分布。退化观测值在简单IG过程模型下的似然函数为：

L\left(\left.Y\right|\mu ,\eta ,{\theta }_{\Lambda }\right)={\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}f\left(\left.\mathrm{\Delta }{y}_{ij}\right|\mu \mathrm{\Delta }{\Lambda }_{ij},\eta \mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}\right)(10)

根据上面给出的先验分布和似然函数，可以得到参数\mu ,\eta ,{\theta }_{\Lambda }的联合后验分布：

\begin{array}{l}p\left(\left.\mu ,\eta ,{\theta }_{\Lambda }\right|Y\right)\propto \mathrm{\pi }\left(\mu \right)\mathrm{\pi }\left(\eta \right)\mathrm{\pi }\left({\theta }_{\Lambda }\right)L\left(Y\mu ,\eta ,{\theta }_{\Lambda }\right)\\ \propto {\mu }^{-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right\}{\eta }^{\delta -\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\gamma \eta \right\}\\ {\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}\sqrt[]{\eta \mathrm{\Delta }{y}_{ij}^{\mathrm{3}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\eta {\left(\mathrm{\Delta }{y}_{ij}-\mu \mathrm{\Delta }{\Lambda }_{ij}\right)}^{\mathrm{2}}}{\mathrm{2}{\mu }^{\mathrm{2}}\mathrm{\Delta }{y}_{ij}}\right]\end{array}(11)

其中，后验分布以概率核的形式给出，省略了与模型参数无关的比例常数部分。显然，无法得到联合后验分布的具体解析表达式，但可以利用MCMC方法在OpenBUGS软件中生成后验样本，以得出产品的可靠性推断：

\begin{array}{l}R\left(\left.t\right|Y\right)=\\ {\int }_{\mu ,\eta ,{\theta }_{\Lambda }}R\left(\left.t\right|\mu \Lambda \left(t\right),\eta {\Lambda }^{\mathrm{2}}\left(t\right)\right)p\left(\left.\mu ,\eta ,{\theta }_{\Lambda }\right|Y\right)\mathrm{d}\mu \mathrm{d}\eta \mathrm{d}{\theta }_{\Lambda }\end{array}(12)

2.2 随机漂移IG过程模型的贝叶斯分析

2.2.1 模型描述

为了解释产品群体内的异质性，在简单IG过程模型中引入随机效应以生成随机效应模型。Ye等^[9]通过让简单IG模型中的\mu遵循截尾正态分布，引入了随机漂移(random drift，RD)IG过程模型。Peng等^[18]通过修改简单IG过程的参数结构，将RD模型改进为Y\left(t\right)~\mathrm{I}\mathrm{G}\left(\mu \Lambda \left(t\right),{\mu }^{\mathrm{3}}{\Lambda }^{\mathrm{2}}\left(t\right)/\eta \right)，该模型的退化均值为\mu \Lambda \left(t\right)，方差为\eta \Lambda \left(t\right)。本文延用Peng等^[18]的模型Y\left(t\right)~\mathrm{I}\mathrm{G}\left(\mu \Lambda \left(t\right),{\mu }^{\mathrm{3}}{\Lambda }^{\mathrm{2}}\left(t\right)/\eta \right),\mu ~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right),\eta >\mathrm{0}，其中对数正态分布的概率密度函数为：

g_u\left(\left.\mu \right|\omega ,{\kappa }^{\mathrm{2}}\right)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\kappa }^{\mathrm{2}}\mu }}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right](13)

RD-IG过程模型Y\left(t\right)的概率密度函数为：

\begin{array}{l}{f}_{\mathrm{R}\mathrm{D}}\left(\left.y\right|\omega ,\kappa ,\Lambda \left(t\right),\eta \right)=\\ {\int }_{\mu >\mathrm{0}}^{}\sqrt[]{\frac{{\mu }^{\mathrm{3}}{\Lambda }^{\mathrm{2}}\left(t\right)}{\mathrm{2}\mathrm{\pi }\eta y^{\mathrm{3}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{\mu {\left(y-\mu \Lambda \left(t\right)\right)}^{\mathrm{2}}}{\mathrm{2}\eta y}\right\}\\ \frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\kappa }^{\mathrm{2}}}\mu }\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{(\mathrm{l}\mathrm{n}\mu -\omega )}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right]\mathrm{d}\mu ,y>\mathrm{0}\end{array}(14)

具有失效阈值D的可靠性函数为：

\begin{array}{l}{R}_{\mathrm{R}\mathrm{D}}\left(\left.t\right|\omega ,\kappa ,\Lambda \left(t\right),\eta \right)={\int }_{\mu >\mathrm{0}}\left\{\Phi \left[\sqrt[]{\frac{\mu }{\eta D}}\left(D-\mu \Lambda \left(t\right)\right)\right]\right.\\ \left.+\mathrm{e}\mathrm{x}\mathrm{p}\left[\frac{\mathrm{2}{\mu }^{\mathrm{2}}\Lambda \left(t\right)}{\eta }\right]\times \Phi \left[-\sqrt[]{\frac{\mu }{\eta D}}\left(D+\mu \Lambda \left(t\right)\right)\right]\right\}\\ \frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\kappa }^{\mathrm{2}}\mu }}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right]\mathrm{d}\mu \end{array}(15)

2.2.2 贝叶斯分析

给定n个单元在某些离散观测时间点的退化观测值为Y_{\mathrm{R}\mathrm{D}}。RD-IG过程模型退化增量\mathrm{\Delta }{y}_{ij}是独立的，并且\mathrm{\Delta }{y}_{ij}~\mathrm{I}\mathrm{G}\left({\mu }_i\mathrm{\Delta }{\Lambda }_{ij},{\mu }_i^{\mathrm{3}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}/\eta \right),{\mu }_i~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right),其

中\mathrm{\Delta }{\Lambda }_{ij}=\Lambda \left(t_{ij}\right)-\Lambda \left(t_{i,j-\mathrm{1}}\right)。

\mathrm{R}\mathrm{D}模型的先验分布为：

\left\{\begin{array}{c}\omega ~N\left({\mu }_{\omega },{\sigma }_{\omega }^{\mathrm{2}}\right)\\ {\kappa }^{\mathrm{2}}~\mathrm{I}\mathrm{G}\mathrm{a}\left(\alpha ,\beta \right)\end{array},\right.\mathrm{ }\left\{\begin{array}{l}{\theta }_{\Lambda }~U\left(a_{\Lambda },b_{\Lambda }\right)\\ \eta ~\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right)\end{array}\right.(16)

其中：N\left({\mu }_{\omega },{\sigma }_{\omega }^{\mathrm{2}}\right)为正态分布，\mathrm{I}\mathrm{G}\mathrm{a}\left(\alpha ,\beta \right)为逆伽马分布，\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right)为伽马分布，U\left(a_{\Lambda },b_{\Lambda }\right)为均匀分布。

参数\mu含随机效应，对所有的{\mu }_i,i=\mathrm{1},\mathrm{2},\cdots ,n，令其遵循具有相同超参数\omega和\kappa的先验分布{\mu }_i~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right)。这些超参数的先验分布如式(16)所示。退化数据Y_{\mathrm{R}\mathrm{D}}在RD-IG过程模型下的似然函数为：

\begin{array}{l}{L}_{\mathrm{R}\mathrm{D}}\left(\left.Y_{\mathrm{R}\mathrm{D}},\mu \right|\omega ,{\theta }_{\Lambda },\eta \right)=\\ {\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}f\left(\left.\mathrm{\Delta }{y}_{ij}\right|{\mu }_i\mathrm{\Delta }{\Lambda }_{ij},{\mu }_i^{\mathrm{3}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}/\eta \right)g_{\mu }\left(\left.{\mu }_i\right|\omega ,{\kappa }^{\mathrm{2}}\right)\end{array}(17)

其中：\mu =\left({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}},\cdots ,{\mu }_n\right)包括退化路径的所有随机参数，并且所有的{\mu }_i具有相同的超参数为\omega和\kappa的分布。参数\omega ,\kappa ,{\theta }_{\Lambda }和\eta以及随机参数\mu的联合后验分布为：

\begin{array}{l}p\left(\left.\omega ,{\kappa }^{\mathrm{2}},{\theta }_{\Lambda },\eta ,\mu \right|Y_{\mathrm{R}\mathrm{D}}\right)\\ \propto \mathrm{\pi }\left(\omega \right)\mathrm{\pi }\left({\kappa }^{\mathrm{2}}\right)\mathrm{\pi }\left({\theta }_{\Lambda }\right)\mathrm{\pi }\left(\eta \right)L_{\mathrm{R}\mathrm{D}}\left(\left.Y_{\mathrm{R}\mathrm{D}},\mu \right|\omega ,\kappa ,{\theta }_{\Lambda },\eta \right)\\ \propto \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\omega -{\mu }_{\omega }\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }_{\omega }^{\mathrm{2}}}\right]{\left({\kappa }^{\mathrm{2}}\right)}^{-\alpha -\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\beta }{{\kappa }^{\mathrm{2}}}\right){\eta }^{\delta -\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma \eta \right)\\ {\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}\sqrt[]{\frac{{\mu }_i^{\mathrm{3}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}}{\eta }}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{\mu }_i{\left(\mathrm{\Delta }{y}_{ij}-{\mu }_i\mathrm{\Delta }{\Lambda }_{ij}\right)}^{\mathrm{2}}}{\mathrm{2}\eta \mathrm{\Delta }{y}_{ij}}\right\}\end{array}(18)

基于联合后验分布，可以得到产品的可靠性推断：

\begin{array}{l}{R}_{\mathrm{R}\mathrm{D}}\left(\left.t\right|Y_{\mathrm{R}\mathrm{D}}\right)=\\ {\int }_{\omega ,\kappa ,\eta ,{\theta }_{\Lambda }}{R}_{\mathrm{R}\mathrm{D}}\left(\left.t\right|\omega ,\kappa ,\Lambda \left(t\right),\eta \right)p\left(\left.\omega ,\kappa ,\eta ,{\theta }_{\Lambda }\right|Y\right)\mathrm{d}\omega \mathrm{d}\kappa \mathrm{d}\eta \mathrm{d}{\theta }_{\Lambda }\end{array}(19)

方程(19)的计算是通过在OpenBUGS软件使用MCMC方法实现的。

2.3 随机波动IG过程模型的贝叶斯分析

2.3.1 模型描述

通过让简单IG模型中的\eta遵循伽马分布且其他参数保持不变，引入了随机波动(random volatility，RV) IG过程模型（Ye等^[9]）。RV模型定义为：Y_{\mathrm{R}\mathrm{V}}\left(t\right)~\mathrm{I}\mathrm{G}\left(\mu \Lambda \left(t\right),\eta {\Lambda }^{\mathrm{2}}\left(t\right)\right),\eta ~\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right),\mu >\mathrm{0}。

伽马分布的PDF为：

g_{\eta }\left(\left.\eta \right|\delta ,\gamma \right)=\frac{{\gamma }^{\delta }{\eta }^{\delta -\mathrm{1}}}{\mathrm{\Gamma }\left(\delta \right)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma \eta \right),\delta >\mathrm{0},\eta >\mathrm{0}(20)

其中\mathrm{\Gamma }\left(\delta \right)={\int }_{\mathrm{0}}^{\mathrm{\infty }}{t}^{\delta -\mathrm{1}}{e}^{-t}\mathrm{d}t是伽马函数。RV-IG过程模型Y\left(t\right)的概率密度函数为：

\begin{array}{l}{f}_{\mathrm{R}\mathrm{V}}\left(\left.y\right|\mu ,\Lambda \left(t\right),\delta ,\gamma \right)=\\ {\int }_{\eta >\mathrm{0}}^{}\sqrt[]{\frac{\eta {\Lambda }^{\mathrm{2}}\left(t\right)}{\mathrm{2}\mathrm{\pi }{y}^{\mathrm{3}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{\eta {\left[y-\mu \Lambda \left(t\right)\right]}^{\mathrm{2}}}{\mathrm{2}{\mu }^{\mathrm{2}}y}\right\}\frac{{\gamma }^{\delta }{\eta }^{\delta -\mathrm{1}}}{\mathrm{\Gamma }\left(\delta \right)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma \eta \right)\mathrm{d}\eta \end{array}(21)

具有失效阈值D的可靠性函数为：

\begin{array}{l}{R}_{\mathrm{R}\mathrm{V}}\left(\left.t\right|\mu ,\Lambda (t),\delta ,\gamma \right)={\int }_{\eta >\mathrm{0}}^{}\left\{\Phi \left[\sqrt[]{\frac{\eta }{D}}\left(\frac{D}{\mu }-\Lambda \left(t\right)\right)\right]\right.+\\ \left.\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}\eta \Lambda (t)}{\mu }\right)\times \Phi \left[-\sqrt[]{\frac{\eta }{D}}\left(\frac{D}{\mu }+\Lambda \left(t\right)\right)\right]\right\}\\ \frac{{\gamma }^{\delta }{\eta }^{\delta -\mathrm{1}}}{\mathrm{\Gamma }\left(\delta \right)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma \eta \right)\mathrm{d}\eta \end{array}(22)

2.3.2 贝叶斯分析

给定n个单元在某些离散观测时间点的退化观测值为Y_{\mathrm{R}\mathrm{V}}。RV-IG过程模型退化增量\mathrm{\Delta }{y}_{ij}=Y_{\mathrm{R}\mathrm{V}}\left(t_{ij}\right)-Y_{\mathrm{R}\mathrm{V}}\left(t_{i,j-\mathrm{1}}\right)是独立的并且\mathrm{\Delta }{y}_{ij}~\mathrm{I}\mathrm{G}\left(\mu \mathrm{\Delta }{\Lambda }_{ij},{\eta }_i\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}\right)，{\eta }_i~\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right),\mathrm{\Delta }{\Lambda }_{ij}=\Lambda \left(t_{ij}\right)-\Lambda \left(t_{i,j-\mathrm{1}}\right)。参数\mu ,{\theta }_{\Lambda },\delta和\gamma的先验分布为：

\left\{\begin{array}{c}\delta ~\mathrm{G}\mathrm{a}\left(a_{\delta },b_{\delta }\right)\\ \gamma ~\mathrm{G}\mathrm{a}\left(a_{\gamma },b_{\gamma }\right)\end{array}\right.,\mathrm{ }\left\{\begin{array}{l}\mu ~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right)\\ {\theta }_{\Lambda }~U\left(a_{\Lambda },b_{\Lambda }\right)\end{array}\right.(23)

由退化数据Y_{\mathrm{R}\mathrm{V}}得到似然函数：

\begin{array}{l}{L}_{\mathrm{R}\mathrm{V}}\left(\left.Y_{\mathrm{R}\mathrm{V}},\eta \right|\mu ,{\theta }_{\Lambda },\delta ,\gamma \right)=\\ {\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}{g}_{\eta }\left(\left.{\eta }_i\right|\delta ,\gamma \right)f\left(\left.\mathrm{\Delta }{y}_{ij}\right|\mu \mathrm{\Delta }{\Lambda }_{ij},{\eta }_i\dot{\mathrm{\Delta }}{\eta }_{ij}^{\mathrm{2}}\right)\end{array}(24)

其中\eta =\left({\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\cdots ,{\eta }_n\right)。基于先验分布和似然函数，得到参数\mu ,{\theta }_{\Lambda },\delta ,\gamma及随机参数\eta的联合后验分布：

\begin{array}{l}p\left(\left.\mu ,{\theta }_{\Lambda },\delta ,\gamma ,\eta \right|Y_{\mathrm{R}\mathrm{V}}\right)\\ \propto \mathrm{\pi }\left(\mu \right)\mathrm{\pi }\left({\theta }_{\Lambda }\right)\mathrm{\pi }\left(\delta \right)\mathrm{\pi }\left(\gamma \right)L_{\mathrm{R}\mathrm{V}}\left(\left.Y_{\mathrm{R}\mathrm{V}},\eta \right|\mu ,{\theta }_{\Lambda },\delta ,\gamma \right)\\ \propto \frac{\mathrm{1}}{\mu }\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right]{\delta }^{a_s-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-b_{\delta }\delta \right\}{\gamma }^{a_i-\mathrm{1}}\\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-b_y\gamma \right){\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}\frac{{\gamma }^{\delta }{\eta }_i^{\delta -\mathrm{1}}}{\mathrm{\Gamma }\left(\delta \right)}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma {\eta }_i\right)\\ \sqrt[]{{\eta }_i\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\eta }_i{\left(\mathrm{\Delta }{y}_{ij}-\mu \mathrm{\Delta }{\Lambda }_{ij}\right)}^{\mathrm{2}}}{\mathrm{2}{\mu }^{\mathrm{2}}\mathrm{\Delta }{y}_{ij}}\right]\end{array}(25)

基于联合后验分布，产品总体的可靠度表达式为：

\begin{array}{l}{R}_{\mathrm{R}\mathrm{V}}\left(\left.t\right|Y_{\mathrm{R}\mathrm{V}}\right)=\\ {\int }_{\mu ,{\theta }_{\Lambda },\delta ,\gamma }{R}_{\mathrm{R}\mathrm{V}}\left(\left.t\right|\mu ,\Lambda (t),\delta ,\gamma \right)p\left(\left.\mu ,{\theta }_{\Lambda },\delta ,\gamma \right|Y_{\mathrm{R}\mathrm{V}}\right)\mathrm{d}\mu \mathrm{d}{\theta }_{\Lambda }\mathrm{d}\delta \mathrm{d}\gamma \end{array}(26)

2.4 随机漂移-波动率IG过程模型的贝叶斯分析

2.4.1 模型描述

考虑到退化率较高的单元会受到退化方差影响，Ye等^[9]引入随机漂移-波动率IG过程(random drift-volatility，RDV)模型。该模型为：Y_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(t\right)~\mathrm{I}\mathrm{G}\left(\mu \Lambda \left(t\right),\eta {\mu }^{\mathrm{2}}{\Lambda }^{\mathrm{2}}\left(t\right)\right),\mu ~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right),\eta >\mathrm{0}。

与2.2节类似，RDV模型的PDF为：

\begin{array}{l}{f}_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(\left.y\right|\omega ,\kappa ,\Lambda \left(t\right),\eta \right)={\int }_{\mu >\mathrm{0}}^{}\sqrt[]{\frac{\eta {\mu }^{\mathrm{2}}{\Lambda }^{\mathrm{2}}\left(t\right)}{\mathrm{2}\mathrm{\pi }{y}^{\mathrm{3}}}}\\ \mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{\eta {\left[y-\mu \Lambda \left(t\right)\right]}^{\mathrm{2}}}{\mathrm{2}y}\right\}\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\kappa }^{\mathrm{2}}}\mu }\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right]\mathrm{d}\mu \end{array}(27)

具有失效阈值D的可靠性函数为：

\begin{array}{l}{R}_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(\left.t\right|\omega ,\kappa ,\Lambda \left(t\right),\eta \right)=\\ {\int }_{\mu >\mathrm{0}}^{}\left\{\Phi \left[\sqrt[]{\frac{\eta }{D}}\left(D-\mu \Lambda \left(t\right)\right)\right]+\mathrm{e}\mathrm{x}\mathrm{p}\left\{\mathrm{2}\eta \mu \Lambda \left(t\right)\right\}\right.\times \Phi \left.\left[-\sqrt[]{\frac{\eta }{D}}\left(D+\mu \Lambda \left(t\right)\right)\right]\right\}\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }{\kappa }^{\mathrm{2}}}\mu }\\ \mathrm{e}\mathrm{x}\mathrm{p}\left\{-\frac{{\left(\mathrm{l}\mathrm{n}\mu -\omega \right)}^{\mathrm{2}}}{\mathrm{2}{\kappa }^{\mathrm{2}}}\right\}\mathrm{d}\mu \end{array}(28)

2.4.2 贝叶斯分析

给定n个单元在离散观测时间点的退化观测值为Y_{\mathrm{R}\mathrm{D}\mathrm{V}}。RDV-IG过程模型退化增量\mathrm{\Delta }{y}_{ij}=Y_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(t_{ij}\right)-Y_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(t_{i,j-\mathrm{1}}\right)是独立的，并且\mathrm{\Delta }{y}_{ij}~\mathrm{I}\mathrm{G}\left({\mu }_i\mathrm{\Delta }{\Lambda }_{ij},\eta {\mu }_i^{\mathrm{2}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}\right),

{\mu }_i~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right),\mathrm{\Delta }{\Lambda }_{ij}=\Lambda \left(t_{ij}\right)-\Lambda \left(t_{i,j-\mathrm{1}}\right)。

RDV模型的先验分布为：

\left\{\begin{array}{c}\omega ~N\left({\mu }_{\omega },{\sigma }_{\omega }^{\mathrm{2}}\right)\\ {\kappa }^{\mathrm{2}}~\mathrm{I}\mathrm{G}\mathrm{a}\left(\alpha ,\beta \right)\end{array},\mathrm{ }\right.\left\{\begin{array}{l}{\theta }_{\Lambda }~U\left(a_{\Lambda },b_{\Lambda }\right)\\ \eta ~\mathrm{G}\mathrm{a}\left(\delta ,\gamma \right)\end{array}\right.(29)

参数\mu含随机效应，对所有的{\mu }_i,i=\mathrm{1},\mathrm{2},\cdots ,n，令其遵循具有相同超参数\omega和\kappa的先验分布{\mu }_i~\mathrm{L}\mathrm{N}\left(\omega ,{\kappa }^{\mathrm{2}}\right)。这些超参数的先验分布如式(29)所示。退化数据Y_{\mathrm{R}\mathrm{D}\mathrm{V}}在RDV-IG过程模型下的似然函数为：

\begin{array}{l}{L}_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(\left.Y_{\mathrm{R}\mathrm{D}\mathrm{V}},\mu \right|\omega ,{\theta }_{\Lambda },\eta \right)=\\ {\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}f\left(\left.\mathrm{\Delta }{y}_{ij}\right|{\mu }_i\mathrm{\Delta }{\Lambda }_{ij},\eta {\mu }_i^{\mathrm{2}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}\right)g_{\mu }\left(\left.{\mu }_i\right|\omega ,{\kappa }^{\mathrm{2}}\right)\end{array}(30)

其中\mu =\left({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}},\cdots ,{\mu }_n\right)包括退化路径的所有随机参数，并且所有的{\mu }_i具有相同的超参数为\omega和\kappa的分布。参数\omega ,\kappa ,{\theta }_{\Lambda }和\eta及随机参数\mu的联合后验分布为：

\begin{array}{l}p\left(\left.\omega ,{\kappa }^{\mathrm{2}},{\theta }_{\Lambda },\eta ,\mu \right|Y_{\mathrm{R}\mathrm{D}\mathrm{V}}\right)\\ \propto \mathrm{\pi }\left(\omega \right)\mathrm{\pi }\left({\kappa }^{\mathrm{2}}\right)\mathrm{\pi }\left({\theta }_{\Lambda }\right)\mathrm{\pi }\left(\eta \right)L_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(\left.Y_{\mathrm{R}\mathrm{D}\mathrm{V}},\mu \right|\omega ,\kappa ,{\theta }_{\Lambda },\eta \right)\\ \propto \mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\left(\omega -{\mu }_{\omega }\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }_{\omega }^{\mathrm{2}}}\right]{\left({\kappa }^{\mathrm{2}}\right)}^{-\alpha -\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\beta }{{\kappa }^{\mathrm{2}}}\right){\eta }^{\delta -\mathrm{1}}\\ \mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma m\right){\displaystyle \prod _{i=\mathrm{1}}^n}{\displaystyle \prod _{j=\mathrm{1}}^{m_i}}\sqrt[]{\eta {\mu }_i^{\mathrm{2}}\mathrm{\Delta }{\Lambda }_{ij}^{\mathrm{2}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{\eta {\left(\mathrm{\Delta }{y}_{ij}-{\mu }_i\mathrm{\Delta }{\Lambda }_{ij}\right)}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }{y}_{ij}}\right]\end{array}(31)

基于联合后验分布，得到产品的可靠性推断：

R_{\mathrm{R}\mathrm{D}\mathrm{V}}{\left(\left.t\right|Y\right)}_{\mathrm{R}\mathrm{D}\mathrm{V}}={\int }_{\omega ,\kappa ,\eta ,{\theta }_{\Lambda }}{R}_{\mathrm{R}\mathrm{D}\mathrm{V}}\left(\left.t\right|\omega ,\kappa ,{\theta }_{\Lambda },\eta \right)p\left(\left.\omega ,\kappa ,\eta ,{\theta }_{\Lambda }\right|Y\right)\mathrm{d}\omega \mathrm{d}\kappa \mathrm{d}\eta \mathrm{d}{\theta }_{\Lambda }(32)

3 模型检验

模型比较准则是为了确定众多模型中的最优模型，常用以下三种准则：最小信息量准则(Akaike information criterion，AIC)、贝叶斯信息(BIC)准则和偏差信息准则(DIC)准则。

1) AIC准则

AIC准则是一种衡量统计模型拟合优度的指标，其核心思想是在对模型进行评估时，既要考虑模型的拟合优度，又要考虑模型的复杂度，以防止过拟合。其表达式如下：

\mathrm{A}\mathrm{I}\mathrm{C}=-\mathrm{2}\mathrm{l}\mathrm{n}\left(L\left(\left.y\right|\theta \right)\right)+\mathrm{2}N(33)

其中：\theta是未知参数；L\left(\left.y\right|\theta \right)是似然函数；N是未知参数的个数。

2) BIC准则

BIC准则指同时考虑模型的拟合优度和模型的复杂度，以便找到最佳的模型平衡点。其表达式如下：

\mathrm{B}\mathrm{I}\mathrm{C}=-\mathrm{2}\mathrm{l}\mathrm{n}\left(L\left(\left.y\right|\theta \right)\right)+N\mathrm{l}\mathrm{n}\left(n\right)(34)

其中：\theta是未知参数；L\left(\left.y\right|\theta \right)是似然函数；N是未知参数的个数；n是样本容量。

3) DIC准则

DIC准则是一种在贝叶斯框架下用于模型选择和模型复杂度评估准则，DIC准则特别适用于从后验分布中抽样得到的模型。其表达式如下：

\mathrm{D}\mathrm{I}\mathrm{C}=\overline{D}+N_D(35)

其中\overline{D}和N_D的关系如下：

\overline{D}=E_{\left.\theta \right|y}\left[-\mathrm{2}\mathrm{l}\mathrm{n}\left(L\left(\left.y\right|\theta \right)\right)\right](36)

N_D=E_{\left.\theta \right|y}\left[-\mathrm{l}\mathrm{n}\left(L\left(\left.y\right|\overline{\theta }\right)\right)\right](37)

其中：\overline{\theta }为\theta的后验均值；L\left(\left.y\right|\overline{\theta }\right)为已确定后验参数情况下的似然函数；其他参数含义与AIC准则中的参数含义一样。

4 模拟研究

令\left(Y_i,T_i,m_i\right),i=\mathrm{1},\mathrm{2},\cdots ,n为按照上面描述的IG过程所生成的随机退化路径。Y_i，T_i和m_i分别表示每条退化路径的退化数据、观测时间点和观测次数。令样本容量n分别等于5、10、30、50、80和100生成组逐渐增加的退化路径。对于每个退化路径\left(Y_i,T_i,m_i\right)，观测时间点的个数m_i=\mathrm{16}，观测时间点T_i均为C_t=c\left(\mathrm{0.25},\mathrm{0.5},\mathrm{0.75},\cdots ,\mathrm{4}\right)。通过对特定IG过程的退化增量进行采样生成退化数据Y_i=Y_i\left(t_{ij}\right)-Y_i\left(t_{i,j-\mathrm{1}}\right),t_{i,\mathrm{0}}=\mathrm{0},Y_i\left(t_{i,\mathrm{0}}\right)=\mathrm{0},j=\mathrm{1},\mathrm{2},\cdots ,m_i,i=\mathrm{1},\mathrm{2},\cdots ,n，此外假设\Lambda \left(t\right)=t^q,q>\mathrm{0}^[9]。

根据上面生成的退化数据，选择无信息先验和信息先验作为先验分布估计参数。其中无信息先验为大区间的均匀分布，此均匀分布的区间选择为参数真实值的10倍，以带参数\omega =\mathrm{2},\kappa =\mathrm{10},q=\mathrm{1.2},\eta =\mathrm{20}的RDV模型为例，参数的无信息先验分布为\omega ~U\left(\mathrm{0},\mathrm{20}\right),\kappa ~U\left(\mathrm{0},\mathrm{100}\right),q~U\left(\mathrm{0},\mathrm{12}\right),\eta ~U\left(\mathrm{0},\mathrm{200}\right)，信息先验分布以第2节所示的形式给出。具体来说，这些先验分布通过让这些分布的均值等于参数真值，方差等于参数真值的1/4来获得，以带参数\omega =\mathrm{2},\kappa =\mathrm{10},q=\mathrm{1.2},\eta =\mathrm{20}的RDV模型为例，参数的信息先验分布为\omega ~N\left(\mathrm{2,0.5}\right),\kappa ~\mathrm{I}\mathrm{G}\mathrm{a}\left(\mathrm{10},\mathrm{2.5}\right),q~U\left(\mathrm{0.25},\left.\mathrm{2.15}\right)\right.，\eta ~\mathrm{G}\mathrm{a}(\mathrm{16},\mathrm{0.8})。

通过OpenBUGS软件基于MCMC模拟生成后验样本，不使用模拟生成的前5 000个样本，使用后10 000个样本。基于MCMC模拟生成的后验样本得到参数的估计，以RDV模型为例，分别给出了无信息先验和信息先验下的参数估计，其结果如表1所示。

表1表明，提出的贝叶斯方法可以非常准确地估计RDV模型的参数，且随着样本容量的增加和信息先验的加入，都能提高估计的精度。特别是对于参数\omega，q，\eta，当样本量从5增加到50时，方差减小很多，对于参数\kappa，信息先验的加入比样本容量的增加对模型的影响更显著，且在其他IG过程模型上也可以证明引入信息先验的优势。

5 算例

用GaAs激光退化数据集证明所提出的贝叶斯方法对IG过程模型的适用性，Wang等^[1]以及Ye等^[9]使用基于最大期望(EM)算法和Bootstrap法的最大似然估计(MLE)对该数据集进行了研究，已经证明，维纳过程和伽马过程都不能很好地适应这个数据集。此外，Ye等^[9]提出的随机漂移IG过程模型与其他IG过程模型具有良好的拟合性。在本节中，使用所提出的贝叶斯方法进一步研究该数据集和IG过程模型，利用贝叶斯方法对退化性能模型进行参数估计的具体步骤如图2所示。

GaAs激光退化数据集为\left(Y_i,T_i,m_i\right),i=\mathrm{1},\mathrm{2},\cdots ,n，其中T_i为第i个样品对应的观测时间，样本容量n=\mathrm{15},m_i=\mathrm{16}。对于所有样品，观测时间点T_i相同，运行过程的退化量如图3所示。

使用所提出的贝叶斯方法将第2节中描述的IG过程模型应用于GaAs激光数据集并假设\Lambda \left(t\right)=t^q,q>\mathrm{0}。为了获得各个参数的估计，IG模型采用无信息先验分布，这些无信息先验以均匀分布的形式给出，通过OpenBUGS软件实现MCMC获得参数的估计结果，对于每个模型，不使用前5 000个样本，使用后10 000个样本，估计结果呈现在表2中。

表2的估计结果表明，贝叶斯方法适用于基于IG过程模型的退化数据分析。另外，RV模型中参数δ和γ的估计结果差异显著，与Wang等^[1]、Ye等^[9]的结果相似。参数的先验信息可以减少估计结果的扩散性，例如由第2节中给出的步骤，可以假定参数δ和γ的先验分布分别为\delta ~\mathrm{G}\mathrm{a}\left(\mathrm{4},\mathrm{0.04}\right)，\gamma ~\mathrm{G}\mathrm{a}\left(\mathrm{4},\mathrm{2}\right)，通过计算得到参数\delta \mathrm{和}\gamma的置信区间分别为\left[\mathrm{49.86},\mathrm{172.2}\right]和\left[\mathrm{0.8955},\right.\left.\mathrm{3.009}\right]，估计结果的改善说明先验信息对最终结果有显著影响。

利用OpenBUGS软件进行MCMC-Gibbs抽样，样本从后验分布中抽取，为消除迭代初值对结果的影响而舍弃前5 000个，利用后10 000个样本用于参数估计，进而得到RD模型的未知参数的后验分布信息，各个未知参数的概率密度函数以及历史迭代分别如图4、图5所示。

从图5可以看出，马尔可夫链样本变化趋于稳定，即说明MCMC-Gibbs抽样收敛，则模型基本收敛。为了验证模型是否合理，可以通过模型比较准则中的DIC准则进行判断，OpenBUGS软件中可以直接获得DIC值，四个模型对应的DIC值如表3所示。

据表3可知，四种模型中RD模型所对应的DIC值最小，因此推荐使用RD模型，结果与Ye等^[9]得到的结论一致。由于选择了RD模型，在模型参数联合后验分布样本的基础上通过式(19)获得GaAs激光器在不同阈值(D=6%、7%、8%、9%)下的可靠性推断，如图6所示。

6 结语

本文采用贝叶斯方法对4种逆高斯过程退化模型进行了研究，利用OpenBUGS软件对模型进行参数估计，再通过DIC准则选择最优模型。案例分析表明：RD模型对数据拟合良好，贝叶斯方法在逆高斯过程模型退化分析中具有灵活性，先验信息的引入可以作为增加样本量的备选方案，且OpenBUGS软件能够准确有效地处理复杂数据。

贝叶斯方法在逆高斯过程模型的处理中面临诸多挑战，如先验知识的选择、计算的复杂性和模型的适用性等。这些挑战激励笔者在未来的工作中继续探索优化贝叶斯方法在逆高斯过程模型中的应用。此外，关于如何将贝叶斯方法与其他统计学习工具进行有机整合，以更好地处理复杂的退化数据，也是未来研究的一个重要方向。

参考文献

[1] WANG X, XU D H. An inverse gaussian process model for degradation data[J]. Technometrics, 2010, 52(2): 188-197.

[2] PENG C Y. Inverse Gaussian processes with random effects and explanatory variables for degradation data[J]. Technometrics, 2015, 57(1): 100-111.

[3] PENG W W, LI Y F, YANG Y J, et al. Bayesian degradation analysis with inverse Gaussian process models under time-varying degradation rates[J]. IEEE Transactions on Reliability, 2017, 66(1): 84-96.

[4] GUO J, WANG C, CABRERA J, et al. Improved inverse Gaussian process and bootstrap: degradation and reliability metrics[J]. Reliability Engineering & System Safety, 2018, 178: 269-277.

[5] HAO S H, YANG J, BERENGUER C. Degradation analysis based on an extended inverse Gaussian process model with skew-normal random effects and measurement errors[J]. Reliability Engineering & System Safety, 2019, 189: 261-270.

[6] PENG W W, ZHU S P, SHEN L J. The transformed inverse Gaussian process as an age- and state-dependent degradation model[J]. Applied Mathematical Modelling, 2019, 75: 837-852.

[7] SUN B, LI Y, WANG Z L, et al. An improved inverse Gaussian process with random effects and measurement errors for RUL prediction of hydraulic piston pump[J]. Measurement, 2021, 173: 108604.

[8] FANG G Q, PAN R, WANG Y K. Inverse Gaussian processes with correlated random effects for multivariate degradation modeling[J]. European Journal of Operational Research, 2022, 300(3): 1177-1193.

[9] YE Z S, CHEN N. The inverse Gaussian process as a degradation model[J]. Technometrics, 2014, 56(3): 302-311.

[10] LI Y F, XIANG Y Y, PAN B S, et al. A hybrid remaining useful life prediction method for cutting tool considering the wear state[J]. The International Journal of Advanced Manufacturing Technology, 2022, 121(5): 3583-3596.

[11] GUO J Y, LI Y F, ZHENG B, et al. Bayesian degradation assessment of CNC machine tools considering unit non-homogeneity[J]. Journal of Mechanical Science and Technology, 2018, 32(6): 2479-2485.

[12] ZHAO Q, JIA X, GUO B. Parameter estimation for the two-parameter exponentiated Weibull distribution based on multiply type-Ⅰ censored data[J]. IEEE Access, 2019, 7: 45485-45493.

[13] 许晓东, 唐圣金, 谢建, 等. 随机退化应力作用下设备剩余寿命预测方法[J]. 兵工学报, 2022, 43(3): 712-719.

[14] 王艺玮, 邓蕾, 郑联语, 等. 基于多通道融合及贝叶斯理论的刀具剩余寿命预测方法[J]. 机械工程学报, 2021, 57(13): 214-224.

[15] 陈龙, 黄天立. 基于贝叶斯更新和逆高斯过程的在役钢筋混凝土桥梁构件可靠度动态预测方法[J]. 工程力学, 2020, 37(4): 186-195.

[16] 杨志远, 赵建民, 李俐莹, 等. 二元相关退化系统可靠性分析及剩余寿命预测[J]. 系统工程与电子技术, 2020, 42(11): 2661-2668.

[17] 沈可, 王芬, 张超, 等. 应用OpenBUGS软件实现网状Meta分析[J]. 湖北医药学院学报, 2013, 32(6): 476-479, 484.

[18] PENG W W, LI Y F, YANG Y J, et al. Inverse Gaussian process models for degradation analysis: a Bayesian perspective[J]. Reliability Engineering & System Safety, 2014, 130: 175-189.

基金资助

国家自然科学基金项目(12361058)

内蒙古自治区直属高校基本科研业务费项目(JY20220165)

内蒙古自治区直属高校基本科研业务费项目(JY20220083)

内蒙古工业大学博士基金(BS2020029)