【文献阅读】A random-effects Wiener degradation model based on accelerated failuretime

一碗姜汤

已于 2024-01-09 12:35:03 修改

阅读量1.6k

点赞数 16

分类专栏：统计学习文章标签：概率论算法

于 2024-01-07 16:46:25 首次发布

本文链接：https://blog.csdn.net/weixin_55252589/article/details/135436113

版权

统计学习专栏收录该内容

21 篇文章 1 订阅

订阅专栏

一、引言

维纳过程是近几十年来最流行的退化模型之一[5,8,20,26,30]。用 $X (t)$ 表示元器件质量特性的退化数据，例如电池的容量，以下维纳过程经常被用来建模其随时间的演变：

$X(t)=vt+\sigma B(t)\quad (1)$

在这个基本的维纳过程模型中， $B(t)$ 是标准的布朗运动，漂移速率 $v > 0$ 对应于退化过程的平均速率，扩散系数 $\sigma > 0$ 量化了过程波动的幅度。基本的维纳过程模型假设一个种群中的所有单位具有相同的漂移率v和扩散系数σ。然而，这种假设在许多实际应用中都存在一定的局限性。特别是，来自同一批的单位，由于原材料的可变性和制造过程中的波动，可能在其降解路径中表现出异质性(heterogeneities)[2,17]。当使用基本的维纳过程来分析这种退化数据时，忽略了个体间的异质性，这可能会导致不准确的统计推断。为了捕获单位间的异质性，随机效应维纳过程允许降解过程的参数是单位特定的，并遵循整个种群的一定分布，这已被证明是有用的[3,21,24,28,29]。在这种情况下，大多数研究假设扩散系数σ是不变的，而每个单位的漂移率v在一个种群中服从正态分布[3,21,29]。一些研究还考虑了扩散系数的异质性。wang[24]允许漂移和扩散是随机的，并假设扩散 $\sigma^2$ 遵循逆伽马分布和v遵循以σ为条件的高斯分布. Ye et al. [28]假设σ是v的一个线性函数，并指定漂移率的倒数为一个高斯分布。用高斯分布描绘漂移速率的异质性的方法虽然被广泛采用，但在以下方面存在限制：

首先，高斯分布在整个实轴上有支持，这通常与漂移率是单边的非正即负的性质相冲突。为了解决这种不相容性，许多研究假设高斯分布的负部分是可以忽略不计的[28]。但是，当v的均值不显著大于v的标准差时，就不能忽略v出现在负半轴的概率。一些研究尝试对高斯分布进行截断[22]，但它无疑会使模型变得更加复杂。
其次，高斯分布具有对称的钟形概率密度函数（PDF），这使得它对单元间的异质性建模具有限制性。例如，我们检查了来自Meeker和Escobar [11]的15个激光器件的降解数据，并用基本维纳过程eq(1)分别拟合每个降解路径，得到关于漂移率的直方图，如图1所示，一个偏态分布似乎更适合于激光器件的具有单位异质性的漂移速率建模。

图1. 15个激光器件在通过一个基本的维纳工艺单独拟合每条路径时的漂移率直方图。

因此，由于这些限制，具有高斯漂移的维纳过程模型的性能可能并不令人满意。为了克服这一不足，本文提出了一种新的随机效应维纳过程模型，利用了加速测试模型的思想。在加速寿命/降解试验（ALT/ADT）中，两组或多个机组在恶劣运行压力下进行测试，以加速故障/降解过程。因此，可以观察到足够的可靠性信息[4,15,23]，否则在实际测试时间内可能难以获得。由于涉及到正常应力的可靠性，因此将较高应力水平下的模型推断到使用应力水平至关重要。加速失效时间（AFT）是一个被广泛接受的概念，其中在加速条件下的退化路径/寿命相当于同一单元在按加速度因子[12,13]适当缩放时间轴后，在正常使用条件下的退化路径/寿命。（加速因子作用在时间轴上这一条非常巧妙）因此，在转换到相同的基线应力水平后，不同应力水平下的降解路径/寿命具有可比性。

受此启发，退化率的异质性可以看作是由未观察到的随机效应引起的随机加速现象，它可以通过时间轴的随机缩放来建模。（该文章的思想核心！）在此情况下，可以使用逆高斯（IG）分布来表征随机效应，并建立了一种新的随机效应维纳过程模型。作为一个随机效应的模型，IG假设自然地适应了漂移率应该为正的实际要求。此外，由于IG分布的形状是相当灵活的[19]，因此它对具有IG随机效应的维纳过程模型具有广泛的适用性。由于该模型继承了加速测试的思想，因此将该模型扩展到ADT数据的分析也很方便。本文的其余部分组织如下。第2节给出了详细的模型公式和相关性质。第3节提出了一个EM算法的最大似然（ML）估计所提出的模型。在第4节中，我们将所提出的模型扩展到恒定应力的ADT建模。第5节说明了所提出的模型在真实退化数据中的应用。结论见第6节。

二、具有IG漂移率v的维纳过程模型

2.1.模型公式

在本研究中，我们考虑以下随机效应维纳过程模型：

$X(t)=v\Lambda(t)+\kappa B(v\Lambda(t))$ (2)

其中：

$v > 0$ 为单位的漂移率，
$\kappa > 0$ 为扩散参数，
$B(\cdot)$ 为标准布朗运动
$\Lambda(t)=\Lambda(t;\theta)$ 是一个具有参数θ的确定函数，用于捕获可能的非线性退化模式[18,26]。例如， $\Lambda(t)$ 可以是线性的： $\Lambda(t)=t$ ，也可以是幂律形式的： $\Lambda(t)=t^{\theta}$ 。通常， $\Lambda(t)$ 是根据退化物理或经验观测[29]来指定的。为了简化符号，在不至于出现混淆的情况下，我们将省略θ。根据惯例，假设 $\Lambda(0)=0$ 和 $\Lambda(t)$ 是单调递增的。

该模型是基于加速测试的思想而提出的。在ALT中，AFT模型是一种广泛应用于量化加速度效应[7,13,14]的模型。在其最简单的形式，线性AFT模型假设寿命分布在不同的压力可以指定

$F_{elevated \ stress}(t|r)=F_0(rt)$

其中，r表示联合加速度效应， $F_0(\cdot)$ 表示在基线应力水平下的寿命的累积分布函数（CDF）。这种关系意味着，通过适当地缩放时间轴，一个产品在一个较高的应力水平下的寿命相当于同一单位在一个基线应力下的寿命。因此，在不同的应力水平下的寿命是统一的，并在转换为基线应力时变得具有可比性。受到这个想法的启发，退化路径中的单位特定的随机效应可以被建模为导致时间轴的随机缩放的随机加速度效应。具体来说，假设每个单位的时间尺度受一个随机尺度因子r的影响，其中r遵循一个均值为1的连续分布。将该公式应用于基本维纳过程模型eq(1)，该单元的降解过程可以建模为：

$X(t)=\mu rt+\sigma B(rt)$ (3)

其中，μ为种群中单位的名义漂移率。由于随机尺度因子r是隐藏的，我们利用可观测漂移率v=μr来重新表述(3)中的模型，并得到(2)中的模型（这里有其不严谨之处，但是无所谓，没有脱离中心思想）。我们在(2)中注意到了所提模型的以下两个特征。首先，观察到的降解速率和扩散系数分别为 $v$ 和 $\kappa^2 v$ ，这内在地解释了降解速率和过程噪声的非均匀性。其次，随机尺度因子引入了观察到的漂移率与扩散系数之间的依赖关系，这适应了实践中常见的现象，即降解率较大的单元往往具有较大的变化[28]。为了模拟单位特定的漂移率v，一个数学上易处理的选择是采用IG分布 $IG(\mu,\zeta)$ ：（搞了半天一直以为是逆伽马分布，结果应该是逆高斯分布，无语....）

$f(v)=\sqrt\frac{\zeta}{2\pi v^3}\mbox{exp}(-\frac{\zeta(v-\mu)^2}{2\mu^2v})$ (4)

其中， $E[v] = \mu$ ， $Var[v]= \mu^3 /\zeta$ 。

除了数学上的便利性外，使用IG分布来建模单位到单位的可变性有相当大的好处。首先，IG分布支持 $(0,+\infty)$ ，这与v > 0的假设相兼容。相反，现有模型中的高斯分布假设只是一个假设v<0的概率可以忽略不计的近似。其次，根据分布参数（μ，ζ），IG分布具有相当灵活的形状，可以适应多种情况。注意，如果我们让ζ→+∞，所提出的模型退化到常用的维纳过程模型，没有随机效应（方差为0，即恒为均值）。

逆伽马分布（Inverse Gamma Distribution）和逆高斯分布（Inverse Gaussian Distribution）都是统计学中的两种分布。为了帮助你更好地理解它们，我们可以通过图形来描述它们的形态。

逆伽马分布（Inverse Gamma Distribution）: 逆伽马分布是伽马分布的逆分布。它的概率密度函数（PDF）为：

$f(x|\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}x^{-\alpha-1}e^{-\beta/x}$

其中，x>0，α>0 和 β>0，Γ 是伽马函数。

逆伽马分布是右偏的（右偏分布的尾部延伸得更长，因为它的平均值、中位数和众数不在同一点上），其形态与正偏的伽马分布相似，但是它是递减的，并且有一个上限。当x 趋近于0时，其值增长得非常快。

逆高斯分布（Inverse Gaussian Distribution）: 逆高斯分布是一个连续分布，它的概率密度函数（PDF）为：

$f(x|\mu,\lambda)=\sqrt{\frac\lambda{2\pi x^3}}\exp\left(-\frac{\lambda(x-\mu)^2}{2\mu^2x}\right)$

其中，x>0，μ>0 和 λ>0。

逆高斯分布是非对称的，它有一个尖峰在 x=μ，并且是右偏的。当 x 远离 μ 时，分布的密度会迅速减小。

逆伽马分布（Inverse Gamma Distribution）InvGamma(α,β) 是统计学中的一种连续概率分布，通常用于描述随机变量的精度（reciprocal of variance）。逆伽马分布有两个参数：形状参数（shape parameter）α 和尺度参数（scale parameter）β。

形状参数 α：

形状参数决定了分布的形状。当 α 较大时，逆伽马分布的峰值会更尖锐，尾部更重。
当 α=1 时，逆伽马分布就简化为一个指数分布。
α 的值应该是大于0的。

尺度参数 β：

尺度参数决定了分布的尺度或范围。当 β 较大时，分布会更加分散；当 β 较小时，分布会更加集中。
β 的值也应该是大于0的。

理解这两个参数的作用可以帮助我们更好地理解逆伽马分布的性质和应用：

当我们考虑一个随机变量的精度时，逆伽马分布提供了一个工具，使我们能够估计这个精度的不确定性或变异性。
通过调整 α 和 β，我们可以控制分布的形状和尺度，从而适应不同的数据特性和分析需求。

在实际应用中，逆伽马分布经常用于贝叶斯统计、贝叶斯回归分析、金融建模等领域，特别是当我们需要对正态分布的精度（即方差的倒数）进行建模和推断时。

在漂移速率v的条件下，单个单位X(t)的退化为高斯分布：

$X(t)|v\sim N(v\Lambda(t),v\kappa^2\Lambda(t))$

因此，X(t)的无条件分布可以通过对v积分得到：

$\begin{aligned} f_{X(t)}(x)& =\int_0^{+\infty}f(x,v)dv=\int_0^{+\infty}f_{X(t)|v}(x|v)f_v(v)dv\\&=\int_{0}^{+\infty}\frac{1}{\sqrt{2\pi\kappa^{2}\nu\Lambda(t)}}\exp\biggl(-\frac{(x-\nu\Lambda(t))^{2}}{2\kappa^{2}\nu\Lambda(t)}\biggr)\sqrt{\frac{\zeta}{2\pi\nu^{3}}}\exp\biggl(-\frac{\zeta(\nu-\mu)^{2}}{2\mu^{2}\nu}\biggr)d\nu \\ &=\frac{\sqrt{\zeta}}{\pi\kappa\sqrt{\Lambda(t)}}\exp\biggl(\frac{x}{\kappa^{2}}+\frac{\zeta}{\mu}\biggr)\mathcal{K}_{-1}(\sqrt{a(t)b(t)})\sqrt{\frac{a(t)}{b(t)}}, \end{aligned}$

其中： $a(t)=\frac{\zeta}{\mu^2}+\frac{\Lambda(t)}{\kappa^2},b(t)=\zeta+\frac{x^2}{\Lambda(t)\kappa^2},$
$\mathcal{K}_p(z)=\frac{1}{2}\int_0^\infty y^{p-1}\exp\biggl(-\frac{z}{2}(y+y^{-1})\biggr)dy$ 是第二类[1]的修正贝塞尔函数。
此外，我们还有： $\operatorname{E}[X(t)]=\mu\Lambda(t),\operatorname{Var}[X(t)]=\mu\kappa^2\Lambda(t)+\frac{\mu^3}\zeta\Lambda(t)^2.$

2.2.到故障发生时的时间分布

在实际应用中，基于退化建模对退化产品进行可靠性推断是一项常规任务。通常，退化X(t)将对应于给定的故障阈值Df, 生命周期定义为X(t)到Df的第一次命中时间(FHT)： $T_f=\inf\{t:X(t>D_f)\}$ 。

在v的条件下，已知转换后的FHT, Λ(Tf) 遵循倒伽马IG分布 $IG(D_f/v,D_f^2/(v\kappa^2))$ 。Λ(Tf)的无条件地PDF可以被推导出为：

$\begin{array}{rcl}f_{\Lambda(T_f)}(u)&=\frac{D_f}{\pi u\kappa\mu}\\\\\end{array}\sqrt{\frac{\zeta(\kappa^2\zeta+\mu^2u)}{\left(\kappa^2\xi u+D_f^2\right)}}\exp\biggl\{\frac{\zeta}{\mu}+\frac{D_f}{\kappa^2}\biggr\}\mathcal{K}_1\biggl(\sqrt{\left(\frac{\zeta}{\mu^2}+\frac{u}{\kappa^2}\right)\biggl(\xi+\frac{D_f^2}{\kappa^2u}\biggr)}\biggr).$ (5)

假设Λ(t)是可微的，则Tf在正常时间t下的PDF为：

$f_{T_f}(t)=f_{\Lambda(T_f)}(\Lambda(t))\frac{d\Lambda(t)}{dt}.$ (6)

可靠性指标，例如，产品的预期寿命可以根据 $f_{T_f}(t)$ 得到。例如，根据幂律退化趋势 $\Lambda(t)=t^{\theta}$ ， $T_f$ 的均值和方差可以推导出为：

$\begin{aligned} \mathbb{E}[T_{f}]& =\frac{2}{\pi}\sqrt{\frac{\zeta}{\kappa^{2}}}\left(\frac{D_{f}}{\mu}\right)^{\frac{1}{2}+\frac{1}{\theta}}\exp\biggl(\frac{D_{f}}{\kappa^{2}}+\frac{\zeta}{\mu}\biggr)\mathcal{K}_{\frac{1}{2}-\frac{1}{\theta}}\biggl(\frac{D_{f}}{\kappa^{2}}\biggr)\mathcal{K}_{\frac{1}{2}+\frac{1}{\theta}}\biggl(\frac{\zeta}{\mu}\biggr), \\ \operatorname{Var}[T_{f}]& =\frac{2}{\pi}\sqrt{\frac{\zeta}{\kappa^{2}}}\left(\frac{D_{f}}{\mu}\right)^{\frac{1}{2}+\frac{2}{\theta}}\exp\left(\frac{D_{f}}{\kappa^{2}}+\frac{\zeta}{\mu}\right)\mathcal{K}_{\frac{1}{2}-\frac{2}{\theta}}\left(\frac{D_{f}}{\kappa^{2}}\right)\mathcal{K}_{\frac{1}{2}+\frac{2}{\theta}}\left(\frac{\zeta}{\mu}\right)-\mathbb{E}\left[T_{f}\right]^{2}. \end{aligned}$

特别是，对于线性退化情况θ = 1和 $\Lambda(t)=t$ ，我们有

$\begin{aligned}\mathbb{E}[T_f]~&=~D_f\bigg(\frac{1}{\mu}+\frac{1}{\zeta}\bigg),\\[2ex]\operatorname{Var}[T_f]~&=~(D_f\kappa^2+D_f^2)\bigg(\frac{1}{\mu\zeta}+\frac{2}{\xi^2}\bigg)+D_f\kappa^2\bigg(\frac{1}{\mu}+\frac{1}{\xi}\bigg)^2.\end{aligned}$

（不得不说，统计的公式的推导是真的solid，这么复杂的期望方差计算。。。）

2.3.从数据当中进一步推断

假设一个单位的退化是在 $\{t_1,...,t_m\}$ 与观察 $X=(X_1,...,X_m)^T$ 之下。定义 $X_0=0$ 和 $t_0 = 0$ 。然后，增量 $\Delta X_j=X_j-X_{j-1}$ 是独立增量和且服从以v为条件的高斯分布：

$(\Delta X_j|v)\sim N(v\Delta \Lambda_j,v\kappa^2\Delta \Lambda_j))$

其中 $\Delta\Lambda_j=\Lambda(t_j)-\Lambda(t_{j-1})$ 。v和X的联合分布很容易得到为：

$p(\nu,\boldsymbol{X})=\sqrt{\frac{\zeta}{2\pi\nu^{3}}}\exp\biggl\{-\frac{\zeta(\nu-\mu)^{2}}{2\mu^{2}\nu}\biggr\}\cdotp\prod_{j=1}^{m}\frac{1}{\sqrt{2\pi\nu\kappa^{2}\Delta\Lambda_{j}}}\exp\biggl\{-\frac{(\Delta X_{j}-\nu\Delta\Lambda_{j})^{2}}{2\nu\kappa^{2}\Delta\Lambda_{j}}\biggr\}.$

通过对v进行积分，我们可以得到X的无条件分布为：

$p\left(\boldsymbol{X}\right)=\sqrt{\frac{\zeta}{(2\pi)^{m+1}\kappa^{2m}}}\prod_{j=1}^{m}\frac{1}{\sqrt{\Delta\Lambda_{j}}}\exp\left(\frac{X_{m}}{\kappa^{2}}+\frac{\zeta}{\mu}\right)2\mathcal{K}_{p}(\sqrt{ab})(b/a)^{p/2},$

其中： $a=\frac{\zeta}{\mu^2}+\frac{\Lambda_m}{\kappa^2},b=\zeta+\frac{1}{\kappa^2}\sum_{j=1}^m\frac{\Delta X_j^2}{\Delta\Lambda_j},p=-\frac{m+1}{2}.$

在许多情况下，根据退化的观测结果来估计v是很有意义的。例如，在复杂系统的预测和运行状况管理中，可以使用运行中的测量的退化信号来更新v的条件分布，以监测系统的运行状况状态。基于上述分析，v对退化观测X的条件分布为：

$p(\nu|\boldsymbol{X})=\frac{(a/b)^{p/2}}{2\mathcal{K}_p(\sqrt{ab})}\nu^{p-1}\exp\biggl(-\frac12(a\nu+b\nu^{-1})\biggr).$ (7)

这表明，以观察数据为条件，v遵循广义IG分布 $GIG(a,b,p)$ [6]（啊？这难道就是我自己要做的贝叶斯后验分布？）。假设我们有n个单位的退化观测数据 $X=\{X_1,...,X_n\}$ ，其中 $X_i=\{X_{i,1},...,X_{i,m_i}\}$ 是在时间点 $t_{i,1},...,t_{i,m_i}$ 获得的。然后，通过最大化对数似然函数，可以估计出未知的模型参数 $\Theta=(\mu,\zeta,\kappa)$ 和 Λ(·)中的参数θ：

$\ell(\mathbf{\Theta},\theta|\mathbf{X})=\sum_{i=1}^n\ln p(\mathbf{X}_i|\mathbf{\Theta},\theta).$ (8)

由于模型参数涉及到修正的贝塞尔函数 $\mathcal{K}$ ，因此给对数似然函数的直接最大化带来了困难。另一方面，单位特定的退化率v可以看作是缺失的数据，可以使用EM算法来寻找ML估计值。EM算法是一种功能强大的ML估计方法，在退化问题中得到了广泛的应用。下面，我们将详细介绍所提模型的参数推理的EM算法。

三、ML估计的EM算法

3.1.点估计

我们首先假设Λ(·)中的参数θ被确定，并重点关注Θ的估计。对于所提出的模型，我们考虑每个单元的不可观测漂移参数 $v_i$ 为缺失数据。用 $V=\{v_1,...,v_n\}$ 表示。然后，我们可以得到完整数据{V，X}的对数似然函数：

$\ell_{\mathfrak{c}}(\mathbf{\Theta}|\mathbf{V},\mathbf{X})=\ell_{\nu}+\ell_{X},$ (9)

其中： $\ell_{\nu}=-\frac{n}{2}\ln(2\pi)+\frac{n}{2}\ln\zeta-\frac{3}{2}\sum_{i=1}^{n}\ln v_{i}-\sum_{i=1}^{n}\frac{\zeta(\nu_{i}-\mu)^{2}}{2\mu^{2}\nu_{i}},$

$\ell_X=\sum_{i=1}^n\left\{-\frac{m_i}{2}\ln(2\pi)-\frac{m_i}{2}\ln{\kappa^2}v_i-\frac{1}{2}\sum_{j=1}^{m_i}\ln{\Delta}\Lambda_{i,j}-\frac{1}{2\kappa^2\nu_i}\sum_{j=1}^{m_i}\frac{(\Delta X_{i,j}-\nu_i\Delta\Lambda_{i,j})^2}{\Delta\Lambda_{i,j}}\right\}.$

在全数据对数似然条件下，EM算法可以迭代地实现E步和M步。假设我们在第s步迭代中有模型参数Θ的估计Θ(s)，需要计算以下Q函数：

$Q(\Theta|\Theta^{(s)})=\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[\ell_{\mathbf{c}}(\mathbf{\Theta}|\mathbf{V},\mathbf{X})],$ (10)

其中，期望是基于（V|X，Θ(s)）的条件分布。请注意，不同单位的退化观测值和漂移率在统计上是相互独立的。因此

$p(\mathbf{V}|\mathbf{X})=\prod_{i=1}^np(v_i|X_i)$

其中，p(vi|Xi)遵循一个具有参数（ai, bi, pi）的广义IG分布，如eq(7)所示。根据广义IG分布的性质，可以很容易地得到：

$\begin{aligned}\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[\nu_i]~&=~\frac{\sqrt{b_i}\mathcal{K}_{p_i+1}(\sqrt{a_ib_i})}{\sqrt{a_i}\mathcal{K}_{p_i}(\sqrt{a_ib_i})},\\\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[\nu_i^{-1}]~&=~\frac{\sqrt{a_i}\mathcal{K}_{p_i+1}(\sqrt{a_ib_i})}{\sqrt{b_i}\mathcal{K}_{p_i}(\sqrt{a_ib_i})}-\frac{2p_i}{b_i}.\end{aligned}$

一旦得到Q函数，模型参数的估计值将通过Q函数的最大化（m步）来更新。推导出Q-函数对Θ的一阶偏导数，并将其设为零，得到以下结果：

$\begin{aligned} &\widehat{\mu}^{(s+1)}&& =\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[v_{i}], \\ &\widehat{\xi}^{(s+1)}&& =\frac1{\frac1n\sum_{i=1}^n\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[\nu_i^{-1}]-1/\widehat{\mu}^{(s+1)}}, \\ &\widehat{\omega}^{(s+1)}&& =\frac{1}{\sum_{i=1}^{n}m_{i}}\sum_{i=1}^{n}\left\{\Lambda_{i,m_{i}}\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[v_{i}]-2X_{i,m_{i}}\right. +\Bigg[\left.\sum_{j=1}^{m_i}\frac{\Delta X_{i,j}^2}{\Delta\Lambda_{i,j}}\Bigg]\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta}^{(s)}}[\nu_i^{-1}]\Bigg\},\right. \end{aligned}$

为了符号简洁，其中ω≜ $\kappa^2$ 表示。执行迭代过程直到连续两次迭代的估计之间的差值小于给定的界限。然后，可以得到Θ的ML估计值。

如果函数Λ(·)中的θ应该从退化数据中得到估计出来，则可以使用轮廓（profile）似然法。

Profile likelihood method（轮廓似然方法）是统计学中一种用于估计参数的方法，特别是在复杂的模型或多参数模型中。它的核心思想是通过固定某个或某些参数，然后关注其余参数的似然函数（likelihood function）的轮廓或截面。

为了更好地解释轮廓似然方法，我们首先要理解似然函数。在统计学中，似然函数是参数的函数，表示在给定参数值下观察到数据的“可能性”。似然函数值越大，意味着参数值越能够解释观察到的数据。

当模型有多个参数时，直接优化整个似然函数可能会很复杂。轮廓似然方法的思想是：

选择一个或一些参数，并认为它们是固定的或已知的。
对于这些固定的参数值，计算其他参数的似然函数。
重复上述步骤，选择不同的参数组合。

通过这种方法，我们可以得到每个参数在其他参数固定时的“轮廓”或“截面”似然。这些轮廓似然可以用来估计参数的不确定性、进行假设检验或进行模型选择。

轮廓似然方法的主要优势在于，它可以减少参数空间的维度，使得估计和推断变得更为可行，特别是在高维参数空间中。

总的来说，轮廓似然方法是一种在多参数统计模型中估计参数和进行推断的技术，它通过固定某些参数并关注其他参数的似然函数来简化复杂的问题。

具体来说，对于任何给定的θ，我们可以首先利用EM算法估计参数Θ。这里，估计的 $\widehat{\Theta}$ 是 $\theta$ 的函数： $\widehat{\Theta}=\widehat{\Theta}(\theta)$ ，将它们带入eq(8)，我们可以得到作为 $\ell(\theta)=\ell(\widehat{\Theta}(\theta),\theta)$ 的轮廓对数似然函数，然后求其最大化，可以得到θ的估计。

3.2.区间估计

假设我们用EM算法得到ML估计值 $\{\widehat{\Theta},\widehat{\theta}\}$ 。模型参数的置信区间可以通过正态渐近法得到。具体地说，观察到的信息矩阵可以作为EM算法的副产品得到。随后，可以得到估计处的渐近协方差作为观测信息矩阵的逆，并可以推导出相关的区间估计。根据Louis [10]，在 $\{\widehat{\Theta},\widehat{\theta}\}$ 上观察到的信息矩阵可以表示为：

$I(\widehat{\Theta},\widehat{\theta})=-\begin{pmatrix}\frac{\partial^2}{\partial\mu^2}&\frac{\partial^2}{\partial\mu\partial\zeta}&\frac{\partial^2}{\partial\mu\partial\omega}&\frac{\partial^2}{\partial\mu\partial\theta^\mathrm{T}}\\\frac{\partial^2}{\partial\mu\partial\zeta}&\frac{\partial^2}{\partial\zeta^2}&\frac{\partial^2}{\partial\zeta\partial\omega}&\frac{\partial^2}{\partial\zeta\partial\theta^\mathrm{T}}\\\frac{\partial^2}{\partial\mu\partial\omega}&\frac{\partial^2}{\partial\zeta\partial\omega}&\frac{\partial^2}{\partial\omega^2}&\frac{\partial^2}{\partial\omega\partial\theta^\mathrm{T}}\\\frac{\partial^2}{\partial\mu\partial\theta}&\frac{\partial^2}{\partial\zeta\partial\theta}&\frac{\partial^2}{\partial\omega\partial\theta}&\frac{\partial^2}{\partial\theta\partial\theta^\mathrm{T}}\end{pmatrix}\ell(\boldsymbol{\Theta},\boldsymbol{\theta})|_{\boldsymbol{\Theta}=\boldsymbol{\widehat\Theta},\boldsymbol{\theta}={\hat{\theta}}}$

$\begin{array}{rcl}=&\{-\mathrm{E}_{\mathbf{V}|\mathbf{X},\widehat{\mathbf{\Theta}},\widehat{\mathbf{\theta}}}[\nabla\nabla^\mathrm{T}\ell_\mathrm{c}(\mathbf{\Theta},\theta|\mathbf{V},\mathbf{X})]-\mathrm{E}_{\mathbf{V}|\mathbf{X},\widehat{\mathbf{\Theta}},\widehat{\mathbf{\theta}}}[\nabla\ell_\mathrm{c}(\mathbf{\Theta},\theta|\mathbf{V},\mathbf{X})\nabla^\mathrm{T}\ell_\mathrm{c}(\mathbf{\Theta},\theta|\mathbf{V},\mathbf{X})]\\\\&+\mathrm{E}_{\mathbf{V}|\mathbf{X},\widehat{\mathbf{\Theta}},\widehat{\mathbf{\theta}}}[\nabla\ell_\mathrm{c}(\mathbf{\Theta},\theta|\mathbf{V},\mathbf{X})]\mathrm{E}_{\mathbf{V}|\mathbf{X},\widehat{\mathbf{\Theta}},\widehat{\mathbf{\theta}}}[\nabla^\mathrm{T}\ell_\mathrm{c}(\mathbf{\Theta},\theta|\mathbf{V},\mathbf{X})]\}_{\boldsymbol{\Theta}=\widehat{\mathbf{\Theta}},\boldsymbol{\theta}=\widehat{\mathbf{\theta}}},\end{array}$ (11)

其中∇是向量微分算子。基于eq(11)和eq(9)中的完全对数似然值，经过一些代数后，我们可以得到以下结果：

$\frac{\partial^{2}\ell}{\partial\mu^{2}}=\frac{\zeta}{\mu^{4}}\sum_{i=1}^{n}(-3\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[\nu_{i}]+2\mu)+\frac{\zeta^{2}}{\mu^{6}}\sum_{i=1}^{n}\mathrm{Var}[\nu_{i}|\mathbf{X},\mathbf{\Theta},\theta],$

$\frac{\partial^{2}\ell}{\partial\mu\partial\zeta}=\frac{1}{\mu^{3}}\sum\limits_{i=1}^{n}\left(\mathsf{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[v_{i}]-\mu\right)-\frac{\zeta}{2\mu^{3}}\sum_{i=1}^{n}\left(\frac{1}{\mu^{2}}\mathrm{Var}[v_{i}|\mathbf{X},\mathbf{\Theta},\theta]+\mathrm{Cov}[v_{i},v_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\theta]\right)$

$\frac{\partial^2\ell}{\partial\mu\partial\omega}=\frac{\zeta}{2\mu^3\omega^2}\sum_{i=1}^n\left\{\Lambda_{i,m_i}\mathrm{Var}[\nu_i|\mathbf{X},\Theta,\theta]+\sum\limits_{j=1}^{m_i}\frac{\Delta X_{i,j}^2}{\Delta\Lambda_{i,j}}\mathrm{Cov}[\nu_i,\nu_i^{-1}|\mathbf{X},\Theta,\theta]\right\},$

$\frac{\partial^2\ell}{\partial\zeta^2}=-\frac{n}{2\zeta^2}+\frac{1}{4\mu^4}\sum\limits_{i=1}^{n}(\mathrm{Var}[\nu_i|\mathbf{X},\Theta,\theta]+2\mu^2\mathrm{Cov}[\nu_i,\nu_i^{-1}|\mathbf{X},\Theta,\theta]+\mu^4\mathrm{Var}[\nu_i^{-1}|\mathbf{X},\Theta,\theta]),$

$\begin{aligned} \frac{\partial^{2}\ell}{\partial\zeta\partial\omega}& =-\frac{1}{4\mu^{2}\omega^{2}}\sum_{i=1}^{n}\left(\Lambda_{i,m_{i}}\mathrm{Var}[v_{i}|\mathbf{X},\Theta,\theta]+\mu^{2}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\mathrm{Var}[v_{i}^{-1}|\mathbf{X},\Theta,\theta]\right) \\ &+\Bigg(\mu^{2}\Lambda_{i,m_{i}}+\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\Bigg)\mathrm{Cov}[v_{i},v_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\theta]\Bigg), \end{aligned}$

$\begin{aligned} \frac{\partial^2\ell}{\partial\omega^2}& =\frac{1}{2\omega^{2}}\sum_{i=1}^{n}m_{i}-\frac{1}{\omega^{3}}\sum_{i=1}^{n}\left\{\Lambda_{i,m_{l}}\mathbb{E}_{\mathbf{V}|\mathbf{X},\Theta,\theta}[\nu_{i}]-2X_{i,m_{l}}+\sum_{j=1}^{mi}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\mathbb{E}_{\mathbf{V}|\mathbf{X},\Theta,\theta}[\nu_{i}^{-1}]\right\} \\ &+\frac{1}{4\omega^{4}}\sum_{i=1}^{n}\left\{\Lambda_{i,m_{i}}^{2}\operatorname{Var}[\nu_{i}|\mathbf{X},\mathbf{\Theta},\theta]+\left[\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\right]^{2}\operatorname{Var}[\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\theta]\right. \\ &+2\Lambda_{i,m_{i}}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\mathrm{Cov}[\nu_{i},\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}]\Biggr\}, \end{aligned}$

其中

$\begin{gathered} \mathrm{Var}[\nu_i|\mathbf{X},\mathbf{\Theta},\theta] =\mathbb{E}_{\mathbf{V|X},\mathbf{\Theta},\theta}[\nu_{i}^{2}]-\mathbb{E}_{\mathbf{V|X},\mathbf{\Theta},\theta}[\nu_{i}]^{2}, \\ \mathrm{Var}[\nu_i^{-1}|\mathbf{X},\mathbf{\Theta},\theta] =\mathbb{E}_{\mathbf{V|X},\Theta,\theta}[\nu_{i}^{-2}]-\mathbb{E}_{\mathbf{V|X},\Theta,\theta}[\nu_{i}^{-1}]^{2}, \\ \mathrm{Cov}[\nu_i,\nu_i^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}] =1-\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[\nu_i]\mathbb{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[\nu_i^{-1}]. \end{gathered}$

根据广义IG分布的性质，我们有：

$\begin{aligned}\mathbb{E}_{\mathbf{V}|\mathbf{x},\mathbf{\Theta},\theta}[\nu_i^2]&=&\frac{b_i\mathcal{K}_{p_i+2}(\sqrt{a_ib_i})}{a_i\mathcal{K}_{p_i}(\sqrt{a_ib_i})},\\\mathbb{E}_{\mathbf{V}|\mathbf{x},\mathbf{\Theta},\theta}[\nu_i^{-2}]&=&\frac{a_i\mathcal{K}_{p_i-2}(\sqrt{a_ib_i})}{b_i\mathcal{K}_{p_i}(\sqrt{a_ib_i})}.\end{aligned}$

此外，涉及θ的导数为：

$\frac{\partial^2\ell}{\partial\mu\partial\theta}=-\frac{\zeta}{2\mu^3\omega}\sum\limits_{i=1}^n\left(\frac{\partial\Lambda_{i,mi}}{\partial\theta}\mathrm{Var}[\nu_i|\mathbf{X},\Theta,\theta]-\sum\limits_{j=1}^{mi}\frac{\Delta X_{i,j}^2}{\Delta\Lambda_{i,j}^2}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\mathrm{Cov}[\nu_i,\nu_i^{-1}|\mathbf{X},\Theta,\theta]\right),$

$\begin{aligned} \frac{\partial^{2}\ell}{\partial\zeta\partial\theta}& =\frac{1}{4\mu^{2}\omega}\sum_{i=1}^{n}\left(\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}\mathrm{Var}[\nu_{i}|\mathbf{X},\mathbf{\Theta},\theta]-\mu^{2}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\mathrm{Var}[\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}]\right. \\ &+\Bigg(\mu^{2}\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}-\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\Bigg)\mathrm{Cov}[\nu_{i},\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\theta]\Bigg), \end{aligned}$

$\begin{aligned} \frac{\partial^{2}\ell}{\partial\omega\partial\theta}& =\frac{1}{2\omega^{2}}\sum_{i=1}^{n}\left(\frac{\partial\Lambda_{i,mi}}{\partial\theta}\mathbb{E}_{\mathbf{v}|\mathbf{x},\mathbf{\Theta},\mathbf{\theta}}[\nu_{i}]-\sum_{j=1}^{m_i}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\mathbb{E}_{\mathbf{v}|\mathbf{x},\mathbf{\theta}}[\nu_{i}^{-1}]\right) \\ &-\frac{1}{4\omega^{3}}\sum_{i=1}^{n}\left(\Lambda_{i,m_{i}}\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}\mathrm{Var}[\nu_{i}|\mathbf{X},\mathbf{\Theta},\theta]\right. \\ &-\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\mathrm{Var}[\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}] \\ &+\Bigg(\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}}-\Lambda_{i,m_{i}}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\Bigg)\mathrm{Cov}[v_{i},v_{i}^{-1}|\mathbf{X},\boldsymbol{\Theta},\boldsymbol{\theta}]\Bigg), \end{aligned}$

$\begin{aligned} \frac{\partial^2\ell}{\partial\theta\partial\theta^\mathrm{T}}& ~-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{m_{i}}\left(\frac{1}{\Delta\Lambda_{i,j}}\frac{\partial^{2}\Delta\Lambda_{i,j}}{\partial\theta\partial\theta^{\mathrm{T}}}-\frac{1}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta^{\mathrm{T}}}\right) \\ &+\frac{1}{2\omega}\sum_{i=1}^{n}\operatorname{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[\nu_{i}^{-1}]\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{3}}\left(\Delta\Lambda_{i,j}\frac{\partial^{2}\Delta\Lambda_{i,j}}{\partial\theta\partial\theta^{\mathrm{T}}}-2\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta^{\mathrm{T}}}\right) \\ &-\frac{1}{2\omega}\sum_{i=1}^{n}\mathsf{E}_{\mathbf{V}|\mathbf{X},\mathbf{\Theta},\theta}[\nu_i]\frac{\partial^{2}\Lambda_{i,m_{i}}}{\partial\theta\partial\theta^{\mathrm{T}}}+\frac{1}{4\omega^{2}}\sum_{i=1}^{n}\left\{\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta^{\mathrm{T}}}\mathrm{Var}[\nu_{i}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}]\right\} \\ &-2\frac{\partial\Lambda_{i,m_{i}}}{\partial\theta}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta^{\mathrm{T}}}\mathrm{Cov}[v_{i},v_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}] \\ &+\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta}\sum_{j=1}^{m_{i}}\frac{\Delta X_{i,j}^{2}}{\Delta\Lambda_{i,j}^{2}}\frac{\partial\Delta\Lambda_{i,j}}{\partial\theta^{\mathrm{T}}}\mathrm{Var}[\nu_{i}^{-1}|\mathbf{X},\mathbf{\Theta},\mathbf{\theta}]\biggr\}, \end{aligned}$

其中：

$\begin{gathered} \frac{\partial\Delta\Lambda_{i,j}}{\partial\theta} =\begin{cases}\dfrac{\partial\Lambda(t_{i,1};\theta)}{\partial\theta},&j=1\\\dfrac{\partial\Lambda(t_{i,j};\theta)}{\partial\theta}-\dfrac{\partial\Lambda(t_{i,j-1};\theta)}{\partial\theta},&j>1\end{cases} \\ \frac{\partial^2\Delta\Lambda_{i,j}}{\partial\theta\partial\theta^\mathrm{T}} =\begin{cases}\dfrac{\partial^2\Lambda(t_{i,1};\theta)}{\partial\theta\partial\theta^\mathrm{T}},&j=1\\\dfrac{\partial^2\Lambda(t_{i,j};\theta)}{\partial\theta\partial\theta^\mathrm{T}}-\dfrac{\partial^2\Lambda(t_{i,j-1};\theta)}{\partial\theta\partial\theta^\mathrm{T}},&j>1\end{cases} \\ \frac{\partial\Lambda_{i,m_i}}{\partial\theta} =\frac{\partial\Lambda(t_{i,m_i};\theta)}{\partial\theta},\frac{\partial^2\Lambda_{i,m_i}}{\partial\theta\partial\theta^\mathrm{T}}=\frac{\partial^2\Lambda(t_{i,m_i};\theta)}{\partial\theta\partial\theta^\mathrm{T}}. \end{gathered}$

显然，当给出Λ(·)中的参数θ时，我们就有了:

$I(\widehat{\mathbf{\Theta}})=-\begin{pmatrix}\frac{\partial^2}{\partial\mu^2}&\frac{\partial^2}{\partial\mu\partial\zeta}&\frac{\partial^2}{\partial\mu\partial\omega}\\\frac{\partial^2}{\partial\mu\partial\zeta}&\frac{\partial^2}{\partial\xi^2}&\frac{\partial^2}{\partial\zeta\partial\omega}\\\frac{\partial^2}{\partial\mu\partial\omega}&\frac{\partial^2}{\partial\zeta\partial\omega}&\frac{\partial^2}{\partial\omega^2}\end{pmatrix}\ell(\Theta)|_{\Theta=\widehat\Theta}$

因此，通过在EM算法的末端进行ML估计，可以很容易地得到观测信息矩阵和渐近协方差。基于ML估计量的渐近协方差，利用delta方法可以进一步得到其他量的区间估计，如时间t时的退化数据X(t)和失效时间Df的p分位数。

四.扩展到ADT数据

考虑K个应力水平 $\xi_1,...,\xi_K$ 下的恒定应力ADT。假设 $n_k$ 单位在应力 $\xi_k$ 中分配以进行降解测试。各应力水平的退化可以用所提出的具有IG漂移的维纳过程模型来建模。从AFT的角度来看，我们假设一个单位在第k个应力水平上的退化可以被建模为：

$X_{k,i}(t)=\mu\cdotp r_{k,i}h(\xi_k)\Lambda(t)+\sigma B(r_{k,i}h(\xi_k)\Lambda(t)).$

这里， $h(\xi_k)$ 是一个确定的依赖于加速度变量的加速度因子。在这样的公式中，我们本质上假设，通过适当地调整一个因子 $h(\xi_k)$ ，一个单位在高应力下的虚拟年龄等同于在正常使用应力下的虚拟年龄。根据可观察到的降解率重新制定模型，我们考虑以下ADT模型：

$X_{k,i}(t)=\nu_{k,i}\Lambda(t)+\kappa\mathcal{B}(\nu_{k,i}\Lambda(t)),$

其中， $\nu_{k,i}=\mu r_{k,i}h(\xi_k)\sim IG(\mu h(\xi_k),\zeta h(\xi_k))$ ，表示单位i在第k个应力水平下的退化率。我们假设应力 $\xi_k$ 已经用 $\xi_0 = 0$ 标准化，并且 $h(\xi_k)$ 具有指数形式的 $h(\xi_k)=\mbox{exp}(\beta\xi_k)$ ，其中β是一个需要确定的适当参数[9]。因此，ADT设置中的模型参数为 $\widehat{\Theta}=\{\mu,\zeta,\kappa,\beta\}$ 和Λ(·)中可能的未知参数θ。

假设我们有第k应力水平下的 $n_k$ 个器件的退化观测数据 $X_{k,1},...,X_{k,n_k}$ ，其中 $X_{k,i}=(X_{k,i,1},...,X_{k,i,m_{k,i}})$ 在时间点 $t_{k,i,1},...,t_{k,i,m_{k,i}}$ 。模型参数可以通过与单应力情况相似地来估计。下面，我们将简要介绍所提出的EM算法到ADT数据的扩展。

4.1.用EM算法进行点估计

在单应力水平的情况下，不同应力水平下的单元的未观察到的退化漂移 $\nu_{k,i}$ 被视为缺失的数据。因此，完全数据对数似然值为

$\begin{aligned} \ell_{\mathsf{c}}=& \text{const+}\frac{1}{2}\sum\limits_{k=1}^{K}\left\{n_{k}\ln(\zeta h(\xi_{k}))-\sum\limits_{i=1}^{n_{k}}\frac{\zeta h(\xi_{k})(\nu_{k,i}-\mu h(\xi_{k}))^{2}}{(\mu h(\xi_{k}))^{2}\nu_{k,i}}\right. \\ &-\sum_{i=1}^{n_{k}}m_{k,i}\ln\omega-\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{k,i}}\ln\Delta\Lambda_{k,i,j}-\frac{1}{\omega}\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{k,i}}\frac{(\Delta X_{k,i,j}-\nu_{k,i}\Delta\Lambda_{k,i,j})^{2}}{\nu_{k,i}\Delta\Lambda_{k,i,j}}\biggr\}, \end{aligned}$

$\begin{aligned} &=\text{const+}\frac{1}{2}\sum_{\mathrm{k=1}}^{\mathrm{K}}\left\{n_{k}\ln\zeta+\beta n_{k}\xi_{k}-\sum_{i=1}^{n_{k}}\frac{\zeta(\nu_{k,i}-\mu\exp(\beta\xi_{k}))^{2}}{\mu^{2}\exp(\beta\xi_{k})\nu_{k,i}}\right. \\ &-\sum_{i=1}^{n_{k}}m_{k,i}\ln\omega-\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{k,i}}\ln\Delta\Lambda_{k,i,j}-\frac{1}{\omega}\sum_{i=1}^{n_{k}}\sum_{j=1}^{m_{k,i}}\frac{(\Delta X_{k,i,j}-\nu_{k,i}\Delta\Lambda_{k,i,j})^{2}}{\nu_{k,i}\Delta\Lambda_{k,i,j}}\Biggr\}. \end{aligned}$

其中 $\mbox{const} = -\frac12\sum_{k=1}^K\sum_{i=1}^{n_k}(1+m_{k,i})\mathrm{ln}(2\pi)\mathrm{~and~}\omega\overset{\Delta}{\operatorname*{=}}\kappa^2.$ 。 $\nu_{k,i}$ 的条件分布仅依赖于第k个应力水平上的第i个单位的数据，即

$p\left(\nu_{1,1},...,\nu_{1,n_1},...\nu_{\mathrm{K},1},...,\nu_{\mathrm{K},n_K}|\mathbf{X}_{1,1},...,\mathbf{X}_{1,n_1},...,\mathbf{X}_{\mathrm{K},1},...,\mathbf{X}_{\mathrm{K},n_K}\right)=\prod_{k=1}^{K}\prod_{i=1}^{n_k}p\left(\nu_{k,i}|\mathbf{X}_{k,i}\right).$

因此，可以得到e步中要求的vki和−vki1的条件期望，与第3节相似。Mstep可以通过推导Q函数对未知参数Θ的一阶导数并将其设为0来实现。经过一些代数计算，我们得到了以下方程式

$\begin{gathered} \mu =\frac{1}{\sum_{k=1}^{K}n_{k}}\sum_{k=1}^{K}\exp(-\beta\xi_{k})\sum_{i=1}^{n_{k}}\mathbb{E}_{\nu_{k,i}|\boldsymbol{X}_{k,i},\boldsymbol{\Theta},\boldsymbol{\theta}}[\nu_{k,i}], \text{(12)} \\ \frac1{\zeta} =\frac{1}{\sum_{k=1}^{K}n_{k}}\sum_{k=1}^{K}\exp(\beta\xi_{k})\sum_{i=1}^{n_{k}}\mathbb{E}_{\nu_{k,i}|X_{k,i},\Theta,\theta}[\nu_{k,i}^{-1}]-\frac{1}{\mu}, \text{(13)} \end{gathered}$

$\frac1\zeta=\frac1{\sum_{k=1}^Kn_k\xi_k}\sum_{k=1}^K\sum_{i=1}^{n_k}\left(\xi_k\exp(\beta\xi_k)\mathsf{E}_{\nu_{k,i}|X_{k,i},\Theta,\theta}[\nu_{k,i}^{-1}]-\frac1{\mu^2}\xi_k\exp(-\beta\xi_k)\mathsf{E}_{\nu_{k,i}|X_{k,i},\Theta,\theta}[\nu_{k,i}]\right).$

替换“等式”（12）进入（13）和（14），β和ζ可以通过将所得到的两个方程相结合来求解。然后，可以得到μ的估计。参数ω = κ2可以更新为

$\begin{aligned}\omega^{(s+1)} & =\frac{1}{\sum_{k=1}^{K}\sum_{i=1}^{n_{k}}m_{k,i}}\\ & \sum_{k=1}^{K}\sum_{i=1}^{n_{k}}\left\{\left[\sum_{j=1}^{m_{i}}\frac{\Delta X_{k,i,j}^{2}}{\Delta\Lambda_{k,i,j}}\right]\mathsf{E}_{\nu_{k,i}|\mathbf{X}_{k,i},\mathbf{\Theta}^{(s)}}[\nu_{k,i}^{-1}]-2X_{k,i,m_{i}}+\Lambda_{k,i,m_{l}}\mathsf{E}_{\nu_{k,i}|\mathbf{X}_{k,i},\mathbf{\Theta}^{(s)}}[\nu_{k,i}]\right\}\end{aligned}$

在需要估计Λ（·）中的参数θ的情况下，可以通过轮廓似然法进行估计，如第3节所述。

4.2.区间估计

置信区间可以用正规渐近方法构造单应力情况。然而，在实践中，应力水平的数量通常很小，在这种情况下，正态近似是不准确的。因此，我们更倾向于采用自助法来推导区间估计值[25]。

我们实现了基于 $\widehat \Theta$ 估计的的参数自助法，如算法1所述。通过自助法，可以得到ML估计的B个类似物。因此，置信区间可以很容易地通过使用这些B类似物的百分位数来构造。

五、示例

5.1.激光仪退化数据

为了说明该模型的应用，我们将其应用于分析来自Meeker和Escobar [11]的砷化镓激光降解数据。激光装置的工作电流随时间增加，当工作电流超过规定的阈值时，激光装置将失效。该数据集由15个测试样本的退化数据组成，每个单元每次测量…{250,500,4000}。这15个机组的工作电流的增加情况如图2所示。从图中可以看出，15个单元的降解率有明显的离散性。因此，我们初步用等式中给出的基本维纳过程模型来拟合每个退化路径(1)，其中对vi和σi的估计数可以很容易地得到为

$\widehat v_i=\frac{X_{i,m_i}}{t_{i,m_i}},\widehat\sigma_i^2=\frac{1}{m_i}\sum_{j=1}^{m_i}\frac{(\Delta X_{i,j}-\widehat\nu_i\Delta t_{i,j})^2}{\Delta t_{i,j}}.$

这里采用线性退化趋势 $\Lambda(t)=t$ ，因为工作电流似乎随时间线性增加。如图1所示，15个单元的估计漂移率显示出一些右偏度。这表明，所提出的模型具有潜在的应用前景。通过简单的计算，可以发现估计的vi与σi 2之间的相关比为0.5，说明vi与σi 2之间呈正相关。因此，我们施加了 $\sigma^2=\kappa_{ts}\nu$ 的关系，并通过以下固定效应维纳过程模型拟合数据

$X_i(t)=\nu_it+\kappa_\mathrm{ts}B(\nu_it).$

$v_{i}$ 和 $\kappa_{ts}^2$ 的ML估计可以通过一些代数推导出：

$\begin{aligned}\hat{\kappa}_\mathrm{ts}^2~&=~\frac{1}{\sum_{i=1}^nm_i}\sum_{i=1}^n\left(\sum_{j=1}^{mi}\frac{\Delta X_{i,j}^2}{\Delta t_{i,j}}\hat{\nu}_i^{-1}-2X_{i,m_i}+\hat{\nu}_it_{i,m_i}\right),\\\hat{v}_i~&=~\frac{1}{2t_{i,m_i}}\bigg(\sqrt{m_i^2\hat{\kappa}_\mathrm{ts}^4+4t_{i,m_i}\sum_{j=1}^{mi}\frac{\Delta X_{i,j}^2}{\Delta t_{i,j}}}-m_i\hat{\kappa}_\mathrm{ts}^2\bigg),i=1,...,n.\end{aligned}$
$\kappa_{ts}^2$ 的估计值为 5.281×10^3 。 $\nu_i$ 的经验CDF如图3所示。经过对经验CDF的检查，我们发现IG分布 $IG(\mu_{ts},\zeta_{ts})$ 可以给 $\nu_i$ 提供一个很好的拟合。

图3. vi的具有90%置信带（CB）的经验CDF和估计IGCDF，i=1…15

参数 $\mu_{ts}$ 和 $\zeta_{ts}$ 可估计为

$\begin{aligned}\widehat{\mu}_\mathrm{ts}&=\frac{1}{n}\sum_{i=1}^n\widehat{\nu}_i=2.037\times10^{-3},\\\widehat{\xi}_\mathrm{ts}&=(\sum_{i=1}^n\widehat{\nu}_i^{-1}-\widehat{\mu}_\mathrm{ts}^{-1})^{-1}=3.010\times10^{-3}.\end{aligned}$

基于上述分析，我们使用所提出的随机效应维纳过程模型来拟合激光数据。模型参数的估计值列在表1中，其中的标准差（SD）是通过正态渐近法得到的。请注意，这里的估计值 $(\widehat \mu,\widehat \zeta, \widehat \kappa^2)$ 非常接近 $(\widehat \mu_{ts},\widehat \zeta_{ts}, \widehat \kappa_{ts}^2)$ 。实际上，我们可以将用上述两步方法得到的 $(\widehat \mu_{ts},\widehat \zeta_{ts}, \widehat \kappa_{ts}^2)$ 视为一个近似估计量。估计值 $(\widehat \mu_{ts},\widehat \zeta_{ts}, \widehat \kappa_{ts}^2)$ 也可以作为EM算法的初始猜测。

表1. 在拟合激光数据时，该模型中的参数的ML估计值和标准差

图4(a)绘制了预期降解 $E[X(t)]=\mu t$ 的估计值和90%置信带。

图4. 左：预期退化趋势的估计

假设失效阈值为 $D_f=6$ ，即，如果工作电流增加6%，则认为激光装置失效。然后根据eq(5)和eq(6)得到寿命分布。图4(b)给出了故障时间分布与激光器件的90%点态置信带。

图4. 右：Df=6的估计寿命分布和90%点态置信带

为了进行比较，我们还通过以下模型对数据进行了拟合：

(i) 常用的具有高斯漂移的随机效应维纳过程模型，如[29]；
(ii) 具有偏正态漂移的维纳过程模型[16]；
(iii) 具有正态-伽马漂移波动[24]的维纳过程模型。

估算结果列于表2中。

表2. 不同模型拟合激光数据时的性能比较

从对数似然值和AIC值可以看出，所提出的模型比现有的模型具有更好的拟合效果。我们还尝试了幂律Λ(·)，即 $\Lambda(t)=t^{\theta}$ 来捕获激光退化过程中任何可能的非线性退化趋势，相应的估计值如表1所示。值得注意的是，估计的θ相当接近于1，而SD相对较小。此外，具有幂律Λ（·）的对数似然值为74.10，仅略大于线性情况。根据AIC准则，结果表明，线性模型较好。因此，具有IG漂移的线性维纳模型是一个很好的激光降解数据模型。

5.2.红外管ADT数据

为了说明所提出的模型在ADT数据分析中的应用，我们考虑了[27，示例8.10]中红外发光二极管（IRLED）的恒定应力ADT数据。IRLED是一种广泛应用于通信系统的高可靠的光电器件。IRLED的性能主要表现为光功率的变化比。为了估计在50 mA运行电流下的可靠性，采样40个单元，并分为两组。一组25个单元在170 mA下测试，一组15个单元在320 mA测试。降解数据见[27]中的表8.9和表8.10。40个单元的降解路径如图5所示。

图5. 40个irled在两级工作电流下的降解路径

可以注意到，即使在相同的应力水平下，单元的退化也表现出相当大的异质性。为了检验随机效应，我们首先应用第3节中提出的模型来分别拟合每个级别的退化数据。表3给出了估计结果：其中幂律退化趋势 $\Lambda(t;\theta)=t^{\theta}$ 遵循[27]。

表3. 两种应力水平下退化的个别估计结果

可以看出，在两个不同水平下的 $\widehat {\theta_k}$ 彼此接近，这表明在高应力条件下，降解机制没有变化。此外，在两个水平下的 $\kappa_k^2$ 估计是接近的，证明了 $\kappa$ 和θ在不同应力水平上是恒定的假设。正如预期的那样，估计的 $\widehat \mu_k$ 随着运行电流的增加而增加。另一方面， $\widehat \xi_k$ 也随着应力水平的增加而增加。一个简单的计算表明， $\widehat \mu_2/\widehat \mu_1$ 和 $\widehat \zeta_2/\widehat \zeta_1$ 的大小相同。这证明了我们基于加速失效时间原理的模型假设。因此，我们应用所提出的模型来拟合IRLED数据，并使用第4节中描述的EM算法获得以下ML估计值：

$\widehat{\mu}=5.457\times10^{-4},\widehat{\zeta}=1.492\times10^{-3},\widehat{\beta}=6.592,\widehat{\kappa}^2=1.450,\widehat{\theta}=0.741.$

我们已经根据[9]标准化了应力，其中退化率使用幂律关系[27]与工作电流相连。基于估计的模型，我们得到了不同应力水平下预期退化路径的点估计和90%置信区间，如图6所示。

图6. 不同工作电流下的预期退化趋势E[X(t)]和90%置信区间。

置信区间由B = 5000重采样本的参数自助法构造。从图中，我们看到，估计的路径符合两个应力水平的退化观察。利用估计的模型，可以将退化规律外推到正常使用应力中，并在此基础上推导出IRLED的失效时间分布。

图7给出了IRLED在 $D_f=10$ 和 $D_f = 20$ 两个阈值下的寿命的CDF。

图7. 在不同的阈值下，IRLED的寿命分布（实线）和90%置信带（虚线）

CDF的上置信界表明，IRLED在使用应激水平上有相当长的寿命。实际上，在故障阈值Df = 10下，预期寿命的90%区间估计数的下限为3.195×105小时（36.5年）。这意味着IRLED具有超高的可靠性，其故障概率在实际应用中可以忽略不计。

六、结论

本研究提出了一种基于线性加速失效时间原理的退化数据随机效应维纳过程模型。新的维纳过程模型利用IG分布来描述种群中的非均质降解速率，并使扩散系数依赖于漂移速率。与现有模型相比，所提出的具有IG漂移的维纳过程模型具有以下特点：

IG分布支持（0，+∞），形状灵活，比正态分布更适合模拟非均匀漂移率；
扩散系数与漂移率呈正相关，这解释了退化率较大的单元在退化路径上往往具有较大的变异性。

对于具有IG漂移的维纳过程模型，提出了一种EM算法来获得模型参数的点估计和区间估计。将该模型进一步扩展到恒应力ADT数据中。通过应用于激光退化数据集和IRLED数据集，研究了具有IG 漂移率的维纳过程模型的效率和适用性。

参考文献

[1] Handbook of mathematical functions with formulas, graphs, and mathematical ta

bles, chap. 9.6 Modi fi ed BesselFunctions I and K. In: Abramowitz M, Stegun IA,

editors. New York: Dover; 1972. p. 374 – 7 .

[2] Fan M, Zeng Z, Zio E, Kang R. Modeling dependent competing failure processes with

degradation-shock dependence. Reliab Eng Syst Saf 2017;165:422 – 30 .

[3] Feng Q, Peng H, Coit D. A degradation-based model for joint optimization of burn

in, quality inspection, and maintenance: a light display device application. Int J Adv

ManufTechnol 2010;50(5 – 8):801 – 8 .

[4] Hu C-H, Lee M-Y, Tang J. Optimum step-stress accelerated degradation test for

wiener degradation process under constraints. Eur J Oper Res 2015;241(2):412 – 21 .

[5] Huang Z, Xu Z, Wang W, Sun Y. Remaining useful life prediction for a nonlinear

heterogeneous wiener process model with an adaptive drift. IEEE Trans Reliab

2015;64(2):687 – 700. https://doi.org/10.1109/TR.2015.2403433 .

[6] Jørgensen B. Statistical properties of the generalized inverse gaussian distribution.

9. Springer Science & Business Media; 1982 .

[7] Kong Y, Ye Z-S. A cumulative-exposure-based algorithm for failure data from a load

sharing system. IEEE Trans Reliab 2016;65(2):1001 – 13 .

[8] Liao H, Elsayed EA. Reliability inference for fi eld conditions from accelerated de

gradation testing. Nav Res Logist 2006;53(6):576 – 87 .

[9] Lim H, Yum B-J. Optimal design of accelerated degradation tests based on wiener

process models. J Appl Stat 2011;38(2):309 – 25 .

[10] Louis TA. Finding the observed information matrix when using the em algorithm. J

R Stat Soc Ser B (Methodological) 1982;44(2):226 – 33. http://www.jstor.org/

stable/2345828

[11] Meeker WQ, Escobar LA. Statistical methods for reliability data. New York: John

Wiley & Sons, INC.; 1998 .

[12] Meeker WQ, Escobar LA, Lu CJ. Accelerated degradation tests: modeling and ana

lysis. Technometrics 1998;40(2):89 – 99 .

[13] Nelson W. Accelerated life testing-step-stress models and data analyses. IEEE Trans

Reliab 1980;R-29(2):103–8. https://doi.org/10.1109/TR.1980.5220742 .

[14] Newby M. Accelerated failure time models for reliability data analysis. Reliab Eng

Syst Saf 1988;20(3):187 –97 .

[15] Park JI, Bae SJ. Direct prediction methods on lifetime distribution of organic light

emitting diodes from accelerated degradation tests. IEEE Trans Reliab

2010;59(1):74 – 90. https://doi.org/10.1109/TR.2010.2040761 .

[16] Peng CY, Tseng ST. Statistical lifetime inference with skew-wiener linear de

gradation models. IEEE Trans Reliab 2013;62(2):338 – 50. https://doi.org/10.1109/

TR.2013.2257055 .

[17] Peng W, Li YF, Mi J, Yu L, Huang HZ. Reliability of complex systems under dynamic

conditions: a bayesian multivariate degradation perspective. Reliab Eng Syst Saf

2016;153:75 – 87 .

[18] Peng W, Li YF, Yang YJ, Mi J, Huang HZ. Bayesian degradation analysis with in

verse gaussian process models under time-varying degradation rates. IEEE Trans

Reliab 2017;66(1):84 – 96. https://doi.org/10.1109/TR.2016.2635149 .

[19] Seshadri V. The inverse gaussian distribution: statistical theory and applications.

137. Springer Science & Business Media; 2012 .

[20] Si X-S, Chen M-Y, Wang W, Hu C-H, Zhou D-H. Specifying measurement errors for

required lifetime estimation performance. Eur J Oper Res 2013;231(3):631 – 44 .

[21] Si X-S, Wang W, Hu C-H, Zhou D-H, Pecht MG. Remaining useful life estimation

based on a nonlinear di ff usion degradation process. IEEE Trans Reliab

2012;61(1):50 – 67 .

[22] Tang S, Yu C, Wang X, Guo X, Si X. Remaining useful life prediction of lithium-ion

batteries based on the wiener process with measurement error. Energies

2014;7(2):520 – 47 .

[23] Tseng S-T, Wen Z-C. Step-stress accelerated degradation analysis for highly reliable

products. J Qual Technol 2000;32(3):209 .

[24] Wang X. Wiener processes with random e ff ects for degradation data. J Multivariate

Anal 2010;101(2):340 – 51 .

[25] Wang X, Xu D. An inverse gaussian process model for degradation data.

Technometrics 2010;52(2):188 – 97 .

[26] Whitmore G, Schenkelberg F. Modelling accelerated degradation data using wiener

di ff usion with a time scale transformation. Lifetime Data Anal 1997;3(1):27 – 45 .

[27] Yang G. Life cycle reliability engineering. John Wiley & Sons; 2007 .

[28] Ye Z-S, Chen N, Shen Y. A new class of wiener process models for degradation

analysis. Reliab Eng Syst Saf 2015;139:58 – 67 .

[29] Ye Z-S, Wang Y, Tsui K-L, Pecht M. Degradation data analysis using wiener pro

cesses with measurement errors. IEEE Trans Reliab 2013;62(4):772 – 80 .

[30] Zhai Q, Ye Z-S, Yang J, Zhao Y. Measurement errors in degradation-based burn-in.

Reliab Eng Syst Saf 2016;150:126 – 35. https://doi.org/10.1016/j.ress.2016.01.015 .

一碗姜汤

关注

16
点赞
踩
23

收藏

觉得还不错? 一键收藏
0
评论
【文献阅读】A random-effects Wiener degradation model based on accelerated failuretime

注意，如果我们让ζ→+∞，所提出的模型退化到常用的维纳过程模型，没有随机效应（方差为0，即恒为均值）。受到这个想法的启发，退化路径中的单位特定的随机效应可以被建模为导致时间轴的随机缩放的随机加速度效应。其次，随机尺度因子引入了观察到的漂移率与扩散系数之间的依赖关系，这适应了实践中常见的现象，即降解率较大的单元往往具有较大的变化[28]。例如，在复杂系统的预测和运行状况管理中，可以使用运行中的测量的退化信号来更新v的条件分布，以监测系统的运行状况状态。从图中，我们看到，估计的路径符合两个应力水平的退化观察。
复制链接

扫一扫