3  基于伽马过程的统计建模

3.1 伽马过程

如果随机过程 \(\{Y(t), t \ge 0\}\) 满足以下性质:

  1. \(Y(0) = 0\);

  2. \(\{Y(t), t \ge 0\}\) 具有平稳和独立增量;

  3. 对任意的两个时间\(t>s >0\), 增量 \(Y_{ts} = Y(t) - Y(s)\) 服从伽马分布伽马分布 \(Ga(\alpha(t-s),\beta)\), 其PDF为 \[ f(y \mid \alpha, \beta) = \frac{\beta^{\alpha(t-s)} y^{\alpha(t-s)-1}}{\Gamma(\alpha(t-s))} \exp\left\{-\beta y \right\}, \] 其中\(\Gamma(\cdot)\) 表示伽马函数, \(\alpha\) 和 \(\beta\) 是未知参数, 且取值大于0. 称随机过程\(\{Y(t), t \ge 0\}\)为齐次伽马过程, 记作 \(Y(t) \sim \mathcal{GP}(\alpha t, \beta)\).

当产品性能退化服从伽马过程\(\mathcal{GP}(\alpha t, \beta)\)且失效阈值为 \(\omega\)时, 定义产品寿命为 \(T=\inf\{t \mid Y(t)\ge\omega\}\). 此时, 寿命\(T\)的CDF可表示为 \[\begin{equation} \begin{aligned} F_{T}(t \mid \alpha,\beta)&=P(T< t)=P(Y_(t)>\omega)=\frac{\Psi(\beta\omega,\alpha t)}{\Gamma(\alpha t)}, \end{aligned} \end{equation}\] 其中 \(\Psi(k,\alpha)\) 是上不完全伽马函数, 定义为 \(\Psi(k,\alpha)=\int_{k}^{\infty}x^{\alpha-1}\exp(-x)\text{d}x\). 本章采用贝叶斯框架进行统计推断, 故将参数\(\alpha\)和\(\beta\)视为随机 变量, 并统一使用条件分布的形式. 尽管\(F_{T}(t \mid \alpha,\beta)\) 可通过常见软件(如 R 或 MATLAB) 直接计算, 但产品的MTTF与寿命分位数缺乏显式解析解, 导致后验推断效率低下, 尤其在实时在线预测场景中计算负担显著. 为此, 可采用 Park 等 (2005) 推荐的方法, 使用两参数 Birnbaum-Saunders (BS) 分布 \(BS(\alpha^\ast,\beta^\ast)\) 来近似 \(F_{T}(t \mid \alpha,\beta)\). BS分布的CDF为 \[\Phi\left(\frac{1}{\alpha^\ast}\left[\sqrt{\frac{x}{\beta^\ast}}-\sqrt{\frac{\beta^\ast}{x}}\right]\right),\] 其中 \(\alpha^\ast=\sqrt{\frac{1}{\beta\omega}}\), \(\beta^\ast=\frac{\beta \omega}{\alpha}\). 此外, 产品的MTTF可近似为 \[\begin{equation} \beta^\ast\left(1+\frac{\left(\alpha^\ast\right)^2}{2}\right)=\frac{1+2\beta\omega}{2\alpha}. \end{equation}\]

本章围绕伽马过程展开建模与统计推断. 第 3.2 节构建伽马过程的贝叶斯分析框架, 并提出多种后验采样算法以适应不同场景下的参数推断需求. 针对产品异质性带来的模型复杂度提升问题, 第 3.3 节提出一种重参数化伽马过程及其变分贝叶斯 (Variational Bayesian, VB) 算法. 该算法通过赋予参数明确物理意义, 显著提升推断效率, 为复杂退化模型的求解提供高效方案.

3.2 伽马过程的贝叶斯分析

在线RUL预测通过退化建模与统计推断实现参数估计. 在伽马过程退化模型领域, Paroissin (2017) 和 Xu 等 (2018) 提出了递推式线性估计方法, 但其在 RUL 预测及置信区间估计中的应用仍存在局限性. 传统离线方法 (如基于贝叶斯框架和极大似然法的策略 (Ling 等, 2019; Wang, 2008; Wang 等, 2021)) 需依赖过去全部的退化数据, 当新增观测时需重新计算, 对存储空间与计算资源需求较高. 为此, 本节介绍了一种低计算量、高效的伽马过程在线 RUL 预测框架. 内容涵盖了以下几个方面:

  1. 贝叶斯推断方法: 推导伽马过程中模型参数的共轭先验, 并以此为基础提出贝叶斯推断方法. 相较于无信息先验与混合先验形式 (Bousquet 等, 2015; Ling 等, 2019), 该方法通过超参数动态更新简化后验计算, 并针对异质性退化场景提供超参数设定准则.

  2. 算法开发: 开发三种后验分布随机样本生成算法, 克服共轭先验结构复杂性与后验采样效率问题. 通过数值模拟实验验证各算法的计算性能与精度优势.

  3. 在线RUL预测算法: 基于共轭先验的解析递推特性, 构建在线 RUL 预测算法. 其流程 ( 详见第 3.2.5 节) 仅需两步:

  • 后验递推更新: 利用当前观测值与历史推断结果更新后验分布;
  • 后验样本生成: 调用预设算法输出参数后验样本. 相比传统方法, 该机制避免全量数据存储与重复计算, 显著提升在线推断效率

本节结构安排如下: 第 3.2.1 节介绍伽马过程的共轭先验及其性质; 第 3.2.2 节介绍三种基于共轭先验的后验样本生成算法; 第 3.2.3 节通过模拟评估算法的估计精度和计算效率; 第 3.2.4 节讨论共轭先验在具有异质效应伽马过程中的扩展; 第 3.2.5 节提出一种在线 RUL 预测算法; 第 3.2.6 节展示新方法在实际案例中的应用.

3.2.1 共轭先验

假设产品性能的退化路径服从伽马过程 \(\mathcal{GP}(\alpha t,\beta)\). 从总体中随机选择 \(n\) 个产品进行测试, 测量时间 为 \(T_1<T_2<\cdots<T_m\). 记第 \(i\) 个产品在时间点 \(T_j\) 的退化值为 \(Y_i(T_j)\), \(i=1,\dots,n,\) \(j=1,\dots,m\). 设 \(Y_{ij}=Y_i(T_j)-Y_i(T_{j-1})\) 且 \(t_j=T_j-T_{j-1}\), 其中 \(Y_i(T_0)=0\) 且 \(T_0=0\), \(i=1,\dots,n,\) \(j=1,\dots,m\). 根据伽马过程的性质 (iii), 有 \(Y_{ij}\sim \textrm{Ga}(\alpha t_j,\beta)\). 令 \(y_{ij}\) 表示 \(Y_{ij}\) 的观测值, 并将观测数据表示为 \(\bm{y}=\{y_{ij}, i=1,\dots,n, j=1,\dots,m\}\). 此时, 基于数据 \(\bm{y}\), 似然函数为 \[\begin{align}\label{like0} L( \bm{y} \mid \alpha,\beta) &=\prod_{i=1}^{n}\prod_{j=1}^{m}\dfrac{\beta^{\alpha t_j}}{\Gamma(\alpha t_j)} y_{ij}^{\alpha t_j-1}\exp\{-\beta y_{ij}\}\nonumber\\ &\propto\dfrac{\beta^{nT_m\alpha}\left[\prod_{i=1}^{n}\prod_{j=1}^{m}y_{ij}^{t_{j}}\right]^\alpha} {\left[\prod_{j=1}^{m}\Gamma(\alpha t_j)\right]^n}\exp\left\{-\beta \sum_{i=1}^{n}\sum_{j=1}^{m}y_{ij}\right\}\nonumber\\ &\propto \dfrac{\beta^{nm\overline{T}_m\alpha}\exp\{-mn\bar{y}_a\beta\}}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{mn}} \left[\prod_{i=1}^{n}\prod_{j=1}^{m}y_{ij}^{\frac{t_{j}}{nm\overline{T}_m}}\right]^{nm\overline{T}_m\alpha}, \end{align}\] 其中\(\overline{T}_m=\frac{T_m}{m}\), \(\bar{y}_a=\frac{1}{mn}\sum_{i=1}^{n}\sum_{j=1}^{m}y_{ij}\).

定理 3.1 基于似然函数 \(\eqref{like0}\), \(\alpha\) 和 \(\beta\) 的共轭先验为 \[\begin{equation}\label{con0} \pi(\alpha,\beta)= C \dfrac{(\beta\omega)^{\delta \overline{T}_m\alpha}}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta}} \exp\{-\delta\lambda\beta\}, \end{equation}\] 其中 \(C\) 是正则化常数, \(\delta\)、\(\omega\) 和 \(\lambda\) 是非负超参数, 分别为先验分布的扩散、形状以及尺度参数.

定理 3.1 的证明见本节附录 3.2.7. 共轭先验 \(\pi(\alpha,\beta)\) 与测量时间相关. 尽管其形式较为复杂, 但可以通过将\(\pi(\alpha, \beta)\) 分解如下形式, 使其更为简洁和易于处理. \[\begin{align}\label{deco0} \pi(\alpha,\beta)&%=\pi(\beta \mid \alpha)\pi(\alpha) \propto \dfrac{(\delta\lambda)^{\delta \overline{T}_m\alpha+1}}{\Gamma\left(1+\delta \overline{T}_m\alpha\right)} \beta^{\delta \overline{T}_m\alpha}\exp\{-\delta\lambda\beta\}\nonumber\\ &\ \ \ \times \dfrac{\left(\dfrac{\omega}{\delta\lambda}\right)^{\delta \overline{T}_m\alpha}\Gamma\left(1+\delta \overline{T}_m\alpha\right)}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta}}. \end{align}\] 可以发现, 在给定 \(\alpha\) 时, \(\beta\)的条件先验 \(\pi(\beta \mid \alpha)\) 服从伽马分布 \(\textrm{Ga}(\delta \overline{T}_m\alpha+1,\delta\lambda)\). 因此, \(\pi(\beta \mid \alpha)\) 的众数、方差和变异系数分别为 \(\overline{T}_m\alpha/\lambda\)、 \((\delta \overline{T}_m\alpha+1)/(\delta\lambda)^2\) 和 \((\delta \overline{T}_m \alpha+1)^{-1/2}\), 其中超参数 \(\lambda\) 是标准的尺度参数, 而超参数 \(\delta\) 则在给定 \(\alpha\) 时决定了方差和变异系数. 当 \(\delta\) 较大时, \(\pi(\beta \mid \alpha)\) 的分布曲将集中于众数附近. 因此, 将 \(\delta\) \(\delta\) 定义为扩散参数. \(\alpha\) 的边际先验形式正比于:
\[\begin{equation*} h(\alpha)= \dfrac{\left(\dfrac{\omega}{\delta\lambda}\right)^{\delta \overline{T}_m\alpha}\Gamma\left(1+\delta \overline{T}_m\alpha\right)}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta}}. \end{equation*}\] 根据 Stirling 公式, 可知: \[ \lim\limits_{u\rightarrow +\infty}\frac{\sqrt{2\pi}u^{u+1/2}\exp\{-u\}}{\Gamma(u+1)}=1. \] 当 \(\alpha\rightarrow+\infty\) 时, 用Stirling 公式对伽马函数做近似处理, 得到以下结果:
\[\begin{align}\label{appro0} h(\alpha)&\equiv\dfrac{\left(\dfrac{\omega}{\delta\lambda}\right)^{\delta \overline{T}_m\alpha} \sqrt{2\pi}(\delta\overline{T}_m\alpha)^{\delta\overline{T}_m\alpha+1/2}\exp\{-\delta\overline{T}_m\alpha\}}{\left[\prod_{j=1}^{m}\left(\sqrt{2\pi}(\alpha t_j)^{\alpha t_j-1/2}\exp\{-\alpha t_j\}\right)^{1/m}\right]^{\delta}}\nonumber\\ &= O\left(\alpha^{\frac{\delta+1}{2}}\exp\left\{-\alpha\delta \overline{T}_m \log\left(\frac{\lambda}{\omega}\frac{\prod_{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\right\}\right), \end{align}\] 其中\(h(\alpha)=O(g(\alpha))\) 表示 \(h(\alpha)\) 和 \(g(\alpha)\) 同阶. 可以证明 \(\log\left(\frac{\prod{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\ge 0\) (详见本节附录 3.2.7). 为了保证 \(\pi(\alpha)\) 是正常的PDF, 一个充分条件是 \(\omega<\lambda\). 同时根据 \(\eqref{appro0}\) 可知 \(\pi(\alpha)\) 的右尾特性类似于伽马分布 \[ \textrm{Ga}\left(\frac{\delta+3}{2},\delta \overline{T}m \left[\log\left(\frac{\lambda}{\omega}\right)+\log\left(\frac{\prod{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\right]\right).\] 在 \(\pi(\alpha)\) 中, \(\omega\) 起到尺度参数的作用, 进一步影响 \(\pi(\beta \mid \alpha)\) 的形状. 因此, \(\omega\) 被称为形状参数. 由于 \(\pi(\beta \mid \alpha)\) 的伽马条件先验和 \(\pi(\alpha)\) 的右尾特性, 这种共轭先验 \(\pi(\alpha,\beta)\) 被称为近似伽马-伽马 (Approximated-gamma-gamma, AGG) 分布, 记为 \(\textrm{AGG}(\delta,\omega,\lambda)\).

图 3.1 展示了当 \(t_j=1\) 且 \(j=1,\dots,m\) 时, 不同 \((\delta,\omega,\lambda)\) 值下 \(\pi(\alpha,\beta)\) 的函数图和等高线图. 以 \(\delta=2\), \(\omega=0.5\), \(\lambda=1.5\) 为基准, 从图中可以看出, 当 \(\delta\) 的值增加到 5 且其他两个超参数保持不变时, 众数的位置几乎相同, 但等高线更集中于众数附近. 增大 \(\lambda\) 的值同样有类似现象, 但众数的位置有所变化. \(\omega\) 则会改变 \(\pi(\alpha,\beta)\) 的形状, 并使其众数的位置发生偏移. 综上所述, 该图展示了 不同参数对 AGG 分布的形状和众数位置有显著影响, 为选择合理的超参数提供了依据.

注1: 当采用等间隔测量(即\(t_j=l\))时, \(\pi(\alpha,\beta)\) 的形式简化为 \[\begin{equation}\label{con1} \pi(\alpha,\beta)=C \dfrac{(\beta\omega)^{\delta l\alpha}}{\left[\Gamma(l\alpha)\right]^{\delta}} \exp{(-\delta\lambda\beta)}. \end{equation}\] 当 \(l=1\) 时, \(\pi(\alpha,\beta)\) 简化为伽马分布 \(Ga(\alpha,\beta)\) 的共轭先验(Damsleth, 1975).

注2: 超参数的取值可以基于先验信息量的强弱来选择. 如图 3.1 所示, 不同的超参数组合会显著影响 \((\alpha,\beta)\)的分布特性: 较大的 \(\delta\)、较小的 \(\omega\) 或较大的 \(\lambda\) 会导致 \((\alpha, \beta)\) 的方差较小, 表示强的先验信息. 在先验知识较少的情况下, 可以选择较小的 \(\delta\)、较大的 \(\omega\) 或较小的 \(\lambda\). 在实际应用中, 可以通过调整 \(\delta\) 控制先验信息的强度. 例如, 根据式 \(\eqref{post0}\), \(\alpha\) 和 \(\beta\) 的后验分布为 \(\textrm{AGG}(\delta_p,\omega_p,\lambda_p)\). 特别地, 对于 \(\omega\) 和 \(\lambda\), 可以使用以下基于观测数据的选择方案:
\[\begin{equation}\label{hyper} \omega=\prod_{i=1}^{n}\prod_{j=1}^{m}y_{ij}^{\frac{t_{j}}{nm\overline{T}_m}}, ~\lambda=\bar{y}_a. \end{equation}\] 这种超参数选择方法属于数据驱动先验的范畴, 其合理性已在贝叶斯统计中得到广泛验证. 例如: 回归系数的先验 (Zellner, 1986), 阈值参数的有效先验 (Hall 等, 2005), 线性退化路径模型的参考先验 (Xu 等, 2012)等. 式 \(\eqref{hyper}\) 所建议的超参数设置具有以下优点: 1. 共轭性保证: 当观测增量中至少有两个值不相等时, 自动满足共轭先验的条件 \(\omega < \lambda\). 2. 众数合理性: \(\omega\) 和 \(\lambda\) 决定了先验分布 \(\pi(\alpha, \beta)\) 的众数位置, 通过数据驱动方式可确保众数值的合理性. 3. 信息量可解释性: 超参数 \(\delta\) 类似于测量次数, 其大小可解释为等效的先验信息量:
\(\delta = 1\): 等效于单个产品的单次测量信息; \(\delta = 0\): 非信息性先验, 表示完全不依赖先验信息; \(\delta = mn\): 与所有观测数据的信息量等效, 表示非常强的先验信息. 基于上述优点, 本章在模拟研究与数据分析中采用自动化策略 \(\eqref{hyper}\) 指定超参数值, 从而简化共轭先验的设定过程.

注3: 在实际工程应用中, 可能已经掌握了一些关于参数 \(\alpha\) 和 \(\beta\) 的先验信息 (如众数、均值、分位数和方差等). 然而, 由于联合先验分布 \(\pi(\alpha,\beta)\) 的形式复杂, 缺乏解析表达式, 直接确定其超参数 \((\delta,\lambda,\omega)\) 可能较为困难. 以下是一种将这些先验信息合理融入到超参数设定中的方法, 分为以下三个步骤:
(1) 设定 \(\delta\) 的值: 如注释 2 所述, 可以通过等效测量信息量, 主观地设定 \(\delta\); (2) 确定 \(\lambda\) 的值: 利用 \(\alpha\) 的信息, 结合 \(\beta\) 的信息以及条件先验分布 \(\pi(\beta \mid \alpha)\) 来确定 \(\lambda\); (3) 确定 \(\omega\) 的值: 根据 \(\alpha\) 的近似分布 \[ Ga\left(\frac{\delta+3}{2},\delta \overline{T}_{m} \left[\log\left(\frac{\lambda}{\omega}\right)+\log\left(\frac{\prod{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\right]\right),\] 结合 \(\alpha\) 的信息设定 \(\omega\). 例如, 对于 \(\alpha\) 和 \(\beta\) 的众数 \(\alpha_M\) 和 \(\beta_M\), 若工程师对其置信度较低, 则可以选择一个相对较小的 \(\delta\) 值. 给定 \(\alpha=\alpha_M\) 时, 可知 \(\pi(\beta \mid \alpha)\) 的众数为 \(\overline{T}_m\alpha_M/\lambda\). 之后可通过求解 \(\beta_M=\overline{T}_m\alpha_M/\lambda\) 得到 \(\lambda=\overline{T}_m\alpha_M/\beta_M\). 由于 \(\alpha\) 的边际先验的近似具有众数为如下形式:
\[\frac{\delta+1}{2\delta \overline{T}_m \left[\log\left(\frac{\lambda}{\omega}\right)+\log\left(\frac{\prod_{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\right]}. \] 给定众数 \(\alpha_M\), 可以得到 \[\omega=\frac{\prod_{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\exp\left\{\log\lambda-\frac{\delta+1}{2\delta \overline{T}_m\alpha_M}\right\}.\]

图 3.1: 共轭先验的PDF特性概览. 左图为 PDF 的完整形态, 右图为其对应的等高线表示.

3.2.2 后验抽样

由于后验分布 \(\pi(\alpha,\beta \mid \bm{y})\) 无法得到解析形式, 直接计算 \(\alpha\) 和 \(\beta\) 的贝叶斯估计变得不可行. 蒙特卡罗方法为处理这种难以解析的后验分布提供了一种可行的替代推断手段. 该方法的核心思想是从参数的联合后验分布中抽取随机样本, 并利用这些样本对参数或其函数进行点估计和区间估计. 本节将介绍三种用于产生AGG分布随机数的算法.

3.2.2.1 Gibbs 采样

Gibbs采样是一种特殊的马尔可夫链蒙特卡罗 (Monte Carlo Markov chain, MCMC) 方法, 其通过从完全条件后验密度 \(\pi(\beta \mid \alpha,\bm{y})\) 和 \(\pi(\alpha \mid \beta,\bm{y})\)迭代抽样来实现. 根据式 \(\eqref{deco0}\), 完全条件后验密度 \(\pi(\beta \mid \alpha,\bm{y})\) 是一个伽马分布 \(Ga(\delta_p\bar{T}_m\alpha,\delta_p\lambda_p)\) ; 而 \(\pi(\alpha \mid \beta,\bm{y})\) 的完全条件后验密度与下式成比例 \[\dfrac{(\beta\omega_p)^{\delta_p \overline{T}_m\alpha}}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta_p}}. \] 可以证明 \(\pi(\alpha \mid \beta,\bm{y})\) 具有对数凹的性质. 因此, 可以利用自适应拒绝采样 (Adaptive synthetic sampling, ARS) 算法生成 \(\pi(\alpha \mid \beta,\bm{y})\) 的随机样本(Gilks 等, 1992). 在获得 \(\alpha\) 和 \(\beta\) 的后验样本后, 可构造任何参数的函数 \(\eta=p(\alpha,\beta)\) (例如: 产品的可靠度、MTTF) 的贝叶斯估计. Gibbs 采样的后验推断过程详见算法 \(\ref{ITR2025-algo1}\).

3.2.2.2 离散网格采样

后验分布 \(\pi(\alpha,\beta \mid \bm{y})\) 可表示为条件分布 \(\pi(\beta \mid \alpha,\bm{y})\) 与边缘分布 \(\pi(\alpha \mid \bm{y})\) 的乘积, 即 \(\pi(\alpha, \beta \mid \bm{y}) = \pi(\beta \mid \alpha, \bm{y}) \pi(\alpha \mid \bm{y})\), 其中 \(\pi(\beta \mid \alpha,\bm{y})\) 为 \(Ga(\delta_p\bar{T}_m\alpha, \delta_p\lambda_p)\), \(\pi(\alpha \mid \bm{y})\) 与以下函数成比例: \[\begin{equation}\label{pmargi} h_p(\alpha)=\dfrac{\left(\dfrac{\omega_p}{\delta_p\lambda_p}\right)^{\delta_p \overline{T}_m\alpha}\Gamma\left(1+\delta_p \overline{T}_m\alpha\right)}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta_p}}. \end{equation}\] 根据上述分解, 当给定从 \(\pi(\alpha \mid \bm{y})\) 生成的随机样本 \(\alpha_0\)时, 可直接从 \(Ga(\delta_p\bar{T}_m\alpha_0,\delta_p\lambda_p)\) 中生成\(\beta\) 的随机样本. 因此, 目前主要挑战在于生成 \(\alpha\) 的随机样本. 由于 \(\pi(\alpha \mid \bm{y})\) 具有复杂性, 本节提出一种简便的离散网格采样法(Discrete grid sampling, DGS): 通过在一系列网格点上构建离散分布来近似 \(\alpha\) 的边缘后验分布, 从而实现 \(\alpha\) 的统计推断. 具体方法如下: 首先, 选取一个区间 \([A_1, A_2]\), 使\(\alpha\)落入该区间的概率接近1, 即, 计算定积分 \(\int_{A_1}^{A_2}\pi(\alpha \mid \bm{y})\text{d}\alpha\)的值是否接近1. 区间 \([A_1, A_2]\)可利用六西格玛法来确定, 具体步骤如下:

  1. 确定 \(\tilde{\alpha}=\mathop{\arg\max}\limits_{\alpha} \log h_p(\alpha)\), 以及 \[I\left(\tilde{\alpha}\right)=-\dfrac{\partial^2\log h_p(\alpha)}{\partial \alpha^2}\biggl|_{\alpha=\tilde{\alpha}}. \]

  2. 根据 (Berger, 2013) 的研究, \(\pi(\alpha \mid \bm{y})\) 可以用正态分布 \(N\left(\tilde{\alpha},\tilde{\sigma}^2\right)\) 近似, 其中 \(\tilde{\sigma}=\sqrt{1/I\left(\tilde{\alpha}\right)}\).

  3. 构建区间 \(A_1=\max{0,\tilde{\alpha}-6\tilde{\sigma}}\) 和 \(A_2=\tilde{\alpha}+6\tilde{\sigma}\). 根据正态分布的性质, 可知 \(\alpha\) 落在区间 \([A_1,A_2]\) 内的概率几乎为 1.

在区间 \([A_1, A_2]\) 内选择 \(M\) 个等间距网格点 \(\left\{A_1=\alpha^{(1)},\alpha^{(2)},\dots,A_2=\alpha^{(M)}\right\}\), 并使用未标准化的后验密度 \(h_p(\alpha)\) 计算每个网格点的概率 \[\begin{equation}\label{disalp} P(\alpha=\alpha^{(s)})=\dfrac{h_p\left(\alpha^{(s)}\right)}{\sum_{i=1}^{M}h_p\left(\alpha^{(i)}\right)}, \quad s=1,\dots,M. \end{equation}\] 通过选择足够大的 \(M\), 可以确保离散分布的近似精度. 通过离散分布采样显著简化了从 \(\pi(\alpha \mid \bm{y})\) 中抽样的难度. 例如, 在 R 语言中, 可以直接使用 函数来模拟随机样本. 上述方法的核心是将复杂的连续分布转换为易操作的离散分布, 显著降低了抽样的计算复杂度. 基于该方法的后验推断过程可见算法 \(\ref{ITR2025-algo2}\).

3.2.2.3 重要性重采样

第三种算法与第二种算法的主要区别在于生成 \(\pi(\alpha \mid \bm{y})\) 后验样本的方法, 该算法采用了重要性重采样 (Sampling Importance Resampling, SIR) 方法. 在 SIR 中, 采样步骤首先从辅助分布 \(g(\alpha)\) 中生成一组随机数, 然后通过重新计算权重, 将其近似为来自 \(\pi(\alpha \mid \bm{y})\) 的样本. 辅助分布 \(g(\alpha)\) 的选择较为灵活, 可以从一组易于抽样的分布中选取. 然而, SIR 方法的效率取决于 \(g(\alpha)\) 和 \(\pi(\alpha \mid \bm{y})\) 的相似程度, 尤其是分布尾部的相似性会影响近似的精度. 从式 \(\eqref{appro0}\) 可知, \(\pi(\alpha \mid \bm{y})\) 的尾部形式 (当 \(\alpha \rightarrow \infty\) 时) 与形状参数 \((\delta_p+3)/2\) 和尺度参数 \[ \nu=\delta_p\overline{T}m \left[\log\left(\frac{\lambda_p}{\omega_p}\right)+\log\left(\frac{\prod{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\right] \] 的伽马分布尾部同阶. 因此, 伽马分布 \(Ga(a, b)\) 是一个合适的辅助分布, 其参数 \(a\) 和 \(b\) 可以通过以下步骤来确定:

  1. 令 \(\tilde{\alpha}=\mathop{\arg\max}\limits_{\alpha} \log h_p(\alpha)\), 并计算 \[ I\left(\tilde{\alpha}\right)=-\dfrac{\partial^2\log h_p(\alpha)}{\partial \alpha^2}\biggl|_{\alpha=\tilde{\alpha}}. \]

  2. 初始化 \(b\) 为 \(b_0=\nu\), 并初始化 \(a\) 为 \(a_0=\tilde{\alpha}b_0\). 确保辅助分布 \(Ga(a_0,b_0)\) 的均值为 \(\tilde{\alpha}\).

  3. 计算精度比 \(R=\frac{b_0^2/a_0}{I\left(\tilde{\alpha}\right)}\), 然后更新参数 \(a=a_0/R\) 和 \(b=b_0/R\). 该步骤不会改变辅助分布的均值, 但可以使 \(Ga(a, b)\) 的方差与 \(\pi(\alpha \mid \bm{y})\) 的渐近方差一致.

通过以上步骤确定了辅助分布 \(Ga(a, b)\), 并从中产生 \(M\) 个随机样本 \(\left\{\alpha^{(1)},\alpha^{(2)},\dots,\alpha^{(M)}\right\}\). 计算权重 \(w_i=h_p\left(\alpha^{(i)}\right)/f_{Ga}\left(\alpha^{(i)} \mid a,b\right)\), \(i=1,\dots,M\), 其中 \(f_{Ga}\left(\alpha^{(i)} \mid a,b\right)\) 表示 \(Ga(a,b)\) 在 \(\alpha^{(i)}\) 处的PDF值. 对权重进行归一化, 得到 \(\tilde{w}_i=w_i/\sum_{j=1}^{M}w_j\). 最后, 从离散分布中生成 \(\alpha\) 的随机样本, 其分布律为 \[\begin{equation}\label{dis2} P\left(\alpha=\alpha^{(i)}\right)=\tilde{w}_i,~i=1,2,\dots,M. \end{equation}\] SIR 算法的后验推断过程可见算法 \(\ref{ITR2025-algo3}\).

3.2.3 模拟实验

在开展模拟研究之前, 首先利用所提方法分析激光退化数据, 其退化路径如图 1.2 所示, 其中 \(n=15\), \(m=16\), \(\omega=10\). 在应用伽马过程拟合数据之前, 首先, 通过Kolmogorov-Smirnov检验来验证数据是否符合伽马过程的 假设(Marsaglia 等, 2003). 由于该数据集中所有 \(t_j\) 的值均为 250 小时, 因此 \(y_{ij} \sim Ga(250\alpha, \beta)\). 使用 R 软件中的 ks.test() 函数进行检验, 得到 p 值为 0.13, 大于显著性水平 0.05. 这表明伽马过程能够用来拟合该数据. 因此, 假设激光器件的退化路径服从伽马过程 \(\mathcal{GP}(\alpha t,\beta)\). 接着, 基于共轭先验 \(\eqref{con0}\) 进行贝叶斯推断, 其中设定 \(\delta=1\) 且 \(\lambda=\bar{y}_a\). 由于数据是等间隔测量的, 因此 \(\omega=\bar{y}_g=\prod_{i=1}^{n}\prod_{1}^{m}y_{ij}^{1/(mn)}\), 即 \(\{y_{ij}, i=1,\dots,15,j=1,\dots,16\}\) 的几何平均值. 在第 3.2.1 节中曾讨论过, \(\delta=1\) 表示先验信息仅相当于一次测量的信息量, 与总计 \(mn\) 次测量数据相比, 先验信息较弱. 因此, \(\alpha\) 和 \(\beta\) 的后验分布为 \(\textrm{AGG}\left(mn+1,\prod_{i=1}^{n}\prod_{1}^{m}y_{ij}^{1/(mn)},\bar{y}_a\right)\). 利用上述三种算法分别得到 \(\alpha\) 和 \(\beta\) 的点估计及其 95% 区间估计, 同时估计 4500 ( 小时) 时设备的可靠度 \(R(4500)\). 估计结果列在表 \(\ref{tbl-tab:dat1}\)中 (其中”GS” 表示基于 Gibbs 采样的算法). 在 Gibbs 采样中, 设定迭代次数 \(K_1\) 为 3000, 其中前 1000 次迭代作为预烧样本舍弃, 并设定稀释(thinning)间隔为2, 最终保留 1000 个有效样本用于后验推断. 在 DGS 中, 离散化区间为 \([0,10]\), 网格点数量设为 10000, 后验推断的样本量同样为 1000. 在 SIR 中, 设定 \(M=10,000\), 并从中重采样 \(K_3=1,000\)个样本. 正如表 \(\ref{tbl-tab:dat1}\) 所示, 三种算法下 \(\alpha\)、\(\beta\) 和 \(R(4500)\) 的贝叶斯点估计及 95%区间估计上均表现出高度一致性.

为了更全面地比较这三种算法, 本节在具有不同信息量的共轭先验下进行了模拟研究. 数据通过伽马过程 \(\mathcal{GP}(\alpha t,\beta)\)生成, 参数设定为 \(\alpha=0.031\) 和 \(\beta=15.35\) (接近表 \(\ref{tbl-tab:dat1}\) 中的估计值). 实验中, 共对 \(n=15\) 个样本进行了测试, 每个样本以 250 小时的间隔进行测量, 总计进行了 \(m=16\) 次测量, 并将失效阈值设定为 10. 采用共轭先验分布 \(\textrm{AGG}(\delta,\bar{y}_g, \bar{y}_a)\), 通过设定\(\delta=0,1,m/4,m/2\) 来评估先验信息量对结果的影响.

随机生成了 \(N=10000\) 个数据集, 并应用所提出的算法进行分析, 得到 \(\alpha\)、\(\beta\)、\(R(4500)\) 和 MTTF 的贝叶斯点估计及其 95% 区间估计. 基于这些估计结果, 计算了相对偏差 (Relative bias, RB) 和 RMSE, 结果分别列于表 \(\ref{tbl-tab:arb}\) 和表 \(\ref{tbl-tab:rmse}\). 在所有情况下, 参数估计均表现出令人满意的效果: \(\alpha\)、\(\beta\) 和 \(R(4500)\) 的贝叶斯估计的 RB 约为 2%, 而 MTTF 的贝叶斯估计的 RB 约为 0.1%. 三种算法在 RB 和 RMSE 上的表现几乎一致, 且先验信息量 \(\delta\) 的变化对估计结果的影响不显著.

在评估参数的区间估计时, 计算了 95% 可信区间的平均长度和CP. 结果分别列在表 \(\ref{tbl-tab:length}\) 和表 \(\ref{tbl-tab:cp}\). 分析发现, 随着先验信息量的增加, 95% 置信区间的长度逐渐减小, 且三种算法在这一趋势上表现相似, 差异不显著. 然而, 在覆盖概率方面, 不同算法表现出一定差异. 对于模型参数 \(\alpha\) 和 \(\beta\), 无论先验信息量 \(\delta\) 的变化如何, 基于DGS和SIR算法的覆盖概率更接近理想的95%水平, 而基于GS算法的覆盖概率则相对较低. 对于 \(R(4500)\) 和 MTTF的估计, 三种算法的覆盖概率均接近名义水平(即, 95%), 表明所提出的后验采样算法在这些指标估计中具有较高的准确性.

在计算效率方面, 三种算法(Gibbs采样、DGS以及SIR)在每个数据集上的平均运行时间分别为0.602 秒、0.00341 秒以及0.00499 秒. 测试环境为运行 Windows 11 操作系统的台式机, 配备 Intel(R) Core(TM) i7-10700 CPU(2.9 GHz) 和 16 GB RAM. 该结果显示, DGS和 SIR 算法的计算效率相近, 并且都显著快于Gibbs算法, 速度提升了超过一百倍. 在实时推断的应用场景中, 计算效率是一个关键因素. 随着新观测数据的不断收集, 需要实时更新后验分布, 同时在保持估计精度的情况下, 尽可能快速完成推断过程. 正如表 \(\ref{tbl-tab:arb}\) 至表 \(\ref{tbl-tab:cp}\) 所示, DGS和 SIR 算法不仅在计算效率上具有明显优势, 而且在估计精度方面表现优异, 充分满足在线推断的需求. 基于上述结果, 在后续章节中, 将主要采用DGS和 SIR 算法用于在线RUL的预测研究.

3.2.4 异质性

产品间的异质性通常源于内在因素和外在因素的差异. 内在因素包括原材料和生产工艺的变化, 而外在因素则涉及操作环境和使用习惯的不同. 这种异质性导致每个产品的性能退化轨迹各不相同. 尽管如此, 由于这些产品源自同一总体, 它们的失效机制通常是一致的. 基于此, 假设第 \(i\) 个产品的性能退化服从伽马过程 \(\mathcal{GP}(\alpha t, \beta_i)\). 这里形状参数\(\alpha\)相同是由于产品失效机制一致, 而不同的尺度参数 \(\beta_i\) 则刻画产品间的异质性. 为简化符号, 假设采用等间隔测量, 即相邻测量时间点的间隔为 \(l\). 试验中测试了 \(n\) 个产品, 在时间点 \(T_m = ml\) 时, 每个产品进行了 \(m\) 次测量. 令 \(Y_{ij}\) 表示第 \(i\) 个产品在时间点 \(T_j = jl\) 的退化值, 其中 \(i=1, \dots, n\), \(j=1, \dots, m\). 定义退化增量为 \(y_{ij} = Y_{ij} - Y_{ij-1}\), 并设 \(Y_{i0} = 0\). 在时间点 \(T_m\)处, 记观测数据为 \(\bm{y_{(m)}} = \{y_{ij}, i=1,\dots,n,~j=1,\dots,m\}\). 由于 \(Y_i(t) \sim \mathcal{GP}(\alpha t, \beta_i)\), 可得 \(y_{ij} \sim Ga(\alpha l, \beta_i)\). 基于观测数据 \(\bm{y_{(m)}}\), 对应的似然函数为 \[\begin{align}\label{like2} L\left( \bm{y_{(m)}} \mid \alpha,\beta_1,\dots,\beta_n\right)&=\prod_{i=1}^{n}\prod_{j=1}^{m}\dfrac{\beta_i^{\alpha l}}{\Gamma(\alpha l)} y_{ij}^{\alpha l-1}\exp\{-\beta_i y_{ij}\}\nonumber\\ &\propto \dfrac{\bar{\beta}_g^{mnl\alpha}}{\left[\Gamma(\alpha l)\right]^{mn}} \bar{y}_{g(m)}^{mnl\alpha}\exp\left\{-\sum_{i=1}^{n}m\bar{y}_{i(m)}\beta\right\}, \end{align}\] 其中 \[\begin{align*} \bar{\beta}_{g}&=\left[\prod_{i=1}^{n}\beta_{i}\right]^{1/n}, \bar{y}_{g(m)}=\left[\prod_{i=1}^{n}\prod_{j=1}^{m}y_{ij}\right]^{\frac{1}{mn}},\\ \bar{y}_{i(m)}&=\frac{1}{m}\sum_{j=1}^{m}y_{ij}, i=1, \dots, n. \end{align*}\]

定理 3.2 基于似然函数, \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的共轭先验为 \[\begin{equation}\label{con2} \begin{aligned} \pi(\alpha,\beta_1,\dots,\beta_n)&= C \dfrac{\left(\bar{\beta}_g\omega\right)^{\delta_1 l\alpha}}{\left[\Gamma(l\alpha)\right]^{\delta_1}} \exp\left\{-\sum_{i=1}^{n}\delta_2\lambda_i\beta_i\right\}, \end{aligned} \end{equation}\] 其中\(C\) 是正则化常数, \(\delta_1\)、\(\delta_2\)、\(\omega\) 和 \(\lambda_i\) 是非负超参数.

定理 3.2 的证明见本节附录 3.2.7. 当 \(\beta_1=\dots=\beta_n=\beta\), \(\delta_1=n\delta_2=\delta\) 且 \(\lambda_1=\dots=\lambda_n=\lambda\) 时, 共轭先验 \(\eqref{con2}\) 简化为 \(\eqref{con1}\). 为了更好地理解共轭先验 \(\eqref{con2}\), 可将 \(\pi(\alpha,\beta_1,\dots,\beta_n)\) 分解为 \[\begin{equation*} \begin{aligned} \pi(\alpha,\beta_1,\dots,\beta_n)&=\prod_{i=1}^{n}\pi(\beta_i \mid \alpha)\pi(\alpha)\\ &\propto\prod_{i=1}^{n}\dfrac{(\delta_2\lambda_i)^{1+\delta_1l\alpha/n}\beta_i^{\delta_1 l\alpha/n}}{\Gamma(1+\delta_1 l\alpha/n)}\exp\{-\delta_2\lambda_i\beta_i\}\\ &\ \ \ \ \times \dfrac{\left[\Gamma(1+\delta_1 l\alpha/n)\right]^n}{\left[\Gamma(l\alpha)\right]^{\delta_1}}\exp\left\{-\alpha\delta_1 l\left[\log\left(\frac{\delta_2}{\omega}\right)+\frac{1}{n}\sum_{i=1}^{n}\log\lambda_i\right]\right\}. \end{aligned} \end{equation*}\] 因此, 给定 \(\alpha\) 时, \(\beta_i\) 的条件分布为 \(Ga(1+\delta_1 l\alpha/n,\delta_2\lambda_i)\), 而 \(\alpha\) 的边缘PDF与下式成比例 \[\begin{align}\label{margi2} g(\alpha)&=\dfrac{\left[\Gamma(1+\delta_1 l\alpha/n)\right]^n}{\left[\Gamma(l\alpha)\right]^{\delta_1}}\nonumber\\ &\quad\times \exp\left\{-\alpha\delta_1 l\left[\log\left(\frac{\delta_2}{\omega}\right)+ \frac{1}{n}\sum_{i=1}^{n}\log\lambda_i\right]\right\}. \end{align}\] 利用Stirling公式, 当 \(\alpha\rightarrow\infty\) 时, 有 \[ g(\alpha)\equiv O\left(\alpha^{\frac{\delta_1+n}{2}} \exp\left\{-A\alpha\right\}\right), \] 其中 \[ A=\delta_1 l\left[\log\left(\frac{n\delta_2}{\delta_1}\right)+\frac{1}{n}\sum_{i=1}^{n}\log\left(\frac{\lambda_i}{\omega}\right)\right]. \] 因此, 当 \(\alpha\rightarrow\infty\)时, \(\pi(\alpha)\)趋于0的速度与伽马分布\(\textrm{Ga}\left(\frac{\delta_1+n+2}{2},K\right)\)同阶. 因此, 称\(\pi(\alpha,\beta_1,\dots,\beta_n)\) 为近似伽马-多元伽马(Approximated-gamma-multivariate gamma, AGMG) 分布, 维度为 \(n\), 记作 \(\textrm{AGMG}_n(\bm\gamma,\omega,\bm\xi)\), 其中 \(\bm\gamma=(\delta_1,\delta_2)^{'}\), \(\bm\xi=(\lambda_1,\dots,\lambda_n)^{'}\).

由式 \(\eqref{post2}\) 可知, 参数 \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的后验分布服从 \(\textrm{AGMG}_n\left(\bm\gamma_{p(m)}, \omega_{p(m)}, \bm\xi_{p(m)}\right)\), 其中 \(\bm\gamma_{p(m)}=(\delta_{1p(m)},\delta_{2p(m)})^{'}\) 和 \(\bm\xi_{p(m)}=(\lambda_{1p(m)},\dots,\lambda_{np(m)})^{'}\). 超参数 \(\omega\) 和 \(\bm\xi\) 的具体选择分别为 \(\bar{y}{g(m)}\) 和 \(\bm{\bar{y}}{(m)}=\left(\bar{y}{1(m)},\dots,\bar{y}{n(m)}\right)^{'}\). 由于后验分布的形式为 \(\textrm{AGMG}_n(\bm\gamma_{p(m)},\bar{y}{g(m)},\bm{\bar{y}}{(m)})\). 在该设置下, 超参数 \(\delta_1\) 和 \(\delta_2\) 与测量次数有类似的解释. 由等式 \(\eqref{con1}\)可知, \(\delta_1\) 和 \(\delta_2\) 主要决定 AGMG 分布的峰度, 从而控制先验信息的权重. 生成AGMG 分布的随机样本可以通过适当修改算法 \(\ref{ITR2025-algo2}\) 和 \(\ref{ITR2025-algo3}\) 来实现. 主要区别在于目标分布变更为 \(\textrm{AGMG}_n\left(\bm\gamma_{p(m)},\omega_{p(m)},\bm\xi_{p(m)}\right)\) 中 \(\alpha\) 的后验边缘分布, 而在给定 \(\alpha\)时, \(\beta_i\)则从伽马分布 \(\textrm{Ga}\left(1+\delta_{1p(m)} l\alpha/n,\delta_{2p(m)}\lambda_{ip(m)}\right)\)中产生随机数, 其中 \(i=1,\dots,n\). 这两种算法 (DGS 和 SIR) 的计算时间与维度 \(n\) 成正比. 为了说明这一点, 本小节在相同的参数设置下, 对 \(n\) 从 2 到 50 的 AGMG 分布应用了这两种算法, 并记录了计算时间, 结果如图 3.2 所示. 从图中可以看出, 计算时间随着 \(n\) 的增加而线性增长. 例如, 当 \(n\) 从 2 增加到 50 时, DGS 的计算时间从 0.00328 秒增加到 0.00895 秒, 而 SIR 的计算时间从 0.00448 秒增加到 0.0101 秒. 这表明即使在较大的 \(n\) 值下, 这两种算法仍然保持着较高的计算效率.

图 3.2: 不同维度下对 AGMG 分布进行采样的两种算法的计算时间.

RUL预测是产品故障诊断与健康管理的核心基础. 在实际应用中, 当产品处于正常工作状态时 ( 即退化量尚未达到失效阈值) , 需要基于实时监测数据动态评估其剩余使用寿命. 具体而言, 假设在时间 \(t_m\) 监测到第 \(i\) 个产品的退化轨迹 \(\{Y_{i1},\dots,Y_{im}\}\) 均满足 \(Y_{ij} < \omega\) ( 其中 \(\omega\) 为失效阈值) , 则其在时间 \(t_m\) 的 RUL 定义为 \[Z_{it_m}=\inf\{z:Y_i(z+t_m)>\omega|Y_{ij}<\omega,j=1,\dots,m\},\] 其分布函数为 \[\begin{align} F_{Z_{it_m}}(z \mid \alpha,\beta_i)&=P(Y_i(z+t_m)>\omega)=P(Y_i(z)>\omega-Y_{im})\nonumber\\ &=\frac{\Psi(\beta(\omega-Y_{im}),\alpha z)}{\Gamma(\alpha z)}. \end{align}\] 由伽马过程的齐次性可知, 最后两个等式成立. 由于\(F_{Z_{it_m}}(z \mid \alpha,\beta_i)\)较为复杂, 可用两参数BS分布\(BS(\alpha_{im}^\ast,\beta_{im}^\ast)\)近似(Park 等, 2005), 其中 \[\alpha_{im}^\ast=\sqrt{\frac{1}{\beta_i(\omega-Y_{im})}},\quad \beta_{im}^\ast=\frac{\beta_i(\omega-Y_{im})}{\alpha}. \] 因此, 在时间点\(t_m\), 第\(i\)个产品RUL的均值可近似为 \[\mu_{im}(\alpha,\beta_i)=\beta_{im}^\ast\left(1+\left(\alpha_{im}^\ast\right)^2/2\right)= \dfrac{1+2\beta_i(\omega-Y_{im})}{2\alpha}.\] 此外, \(Z_{i t_m}\) 的分布 \(\rho\) 分位数可近似为 \[\mu^{\rho}_{im}(\alpha,\beta_i)=\frac{\beta_{im}^\ast}{4}\left[u_\rho\alpha_{im}^\ast +\sqrt{\left(u_\rho\alpha_{im}^\ast\right)^{2}+4}\right]^{2},\] 其中 \(u_\rho\) 是标准正态分布的 \(\rho\) 分位数. 在时间\(t_m\), 第 \(i\) 个产品RUL 的贝叶斯点预测为 \[\begin{equation}\label{point} \tilde{\mu}_{im}=\int_0^{\infty}\int_0^{\infty}\mu_{im}(\alpha,\beta_i)\pi(\alpha,\beta_i \mid \bm{y_{(m)}})\text{d}\alpha\text{d}\beta_i. \end{equation}\] 在时间\(t_m\)时, 第 \(i\) 个产品 RUL 的 \(1-\rho\) 可信水平的贝叶斯区间预测为 \[\begin{equation}\label{inter} \left(\tilde{\mu}_{im}^{\rho/2},\tilde{\mu}_{im}^{1-\rho/2}\right), \end{equation}\] 其中 \[ \tilde{\mu}_{im}^{\rho}=\int_0^{\infty}\int_0^{\infty}\mu^{\rho}_{im}(\alpha,\beta_i)\pi(\alpha,\beta_i \mid \bm{y_{(m)}})\text{d}\alpha\text{d}\beta_i. \] 给定后验样本 \({(\alpha^{(k)},\beta_i^{(k)}), k=1,\dots,K}\), 式 \(\eqref{point}\) 和 \(\eqref{inter}\) 可通过蒙特卡罗方法近似: \[\begin{equation}\label{mcest} \tilde{\mu}_{im}\approx\frac{1}{K}\sum_{k=1}^{K}\mu_{im}(\alpha^{(k)},\beta_i^{(k)}),~~ \tilde{\mu}_{im}^{\rho}\approx\frac{1}{K}\sum_{k=1}^{K}\mu^{\rho}_{im}(\alpha^{(k)},\beta_i^{(k)}). \end{equation}\]

注4: DGS 和 SIR 算法主要用于产生模型参数 \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的后验样本. 基于这些后验样本, 可通过式 \(\eqref{mcest}\) 对 \(n\) 个产品进行 RUL 预测. 算法具备高度灵活性, 既适用于单个产品的预测, 也适用于多个产品的协同预测. 当 \(n=1\) 时, 算法仅利用单个产品的信息, 退化为第 3.2.2 节中的内容. 当 \(n\ge 2\) 时, 算法通过共享多个产品的信息来共同估计参数 \(\alpha\), 并利用其他产品的信息提高参数 \(\beta_i\) 的估计精度.

3.2.5 RUL预测

假设在时间点\(t_{m+1}=(m+1)l\)收集到了 \(n\) 个产品的新退化增量 \(\bm{y}_{(m+1)}=(y_{1m+1},\dots,y_{nm+1})^{'}\). 在获取新的观测数据后, \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的后验分布需要更新. 对于具有共轭先验的贝叶斯推断, 可使用递推公式来实现更新. 由式 \(\eqref{post2}\) 可知, 在时间 \(t_{m+1}=(m+1)l\) 时, \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的后验分布为 \(\textrm{AGMG}_n\left(\bm\gamma_{p(m+1)},\omega_{p(m+1)},\bm\xi_{p(m+1)}\right)\), 其中参数 \(\bm\gamma_{p(m+1)}\)、\(\omega_{p(m+1)}\) 和 \(\bm\xi_{p(m+1)}\) 可实现递推更新. 即 \[\begin{equation}\label{post3} \begin{aligned} \bm\gamma_{p(m+1)}&=\bm\gamma_{p(m)}+(n,1)^{'}, \\ \omega_{(m+1)}&=\omega_{(m)}^{\frac{mn+\delta_1}{(m+1)n+\delta_1}}\left[\prod_{i=1}^{m}y_{im+1}\right]^{\frac{1}{(m+1)n+\delta_1}},\\ \bm\lambda_{(m+1)}&=\frac{m+\delta_2}{m+1+\delta_2}\bm\lambda_{(m)}+\frac{1}{m+1+\delta_2}\bm{y}_{(m+1)}. \end{aligned} \end{equation}\] 由式 \(\eqref{post3}\)可知, 后验分布的更新仅依赖\(t_m\)时刻后验分布中的参数值和新增观测数据, 可显著提升计算和数据存储效率. 本章提出的高效 DGS 和 SIR 算法, 可快速生成 AGMG 分布样本, 并用于多个产品的在线 RUL 预测. 具体过程如下:

3.2.6 实例分析

3.2.6.1 激光退化数据

基于第 \(\ref{sec-ITR2025-sim}\) 节的激光退化数据集 ( 见图 1.2) , 我们重点分析第 1、6 和 10 个激光设备. 观测数据显示, 这三个设备在运行4000小时后均因性能退化量超过阈值 \(\omega=10\)而发生失效. 由于监测间隔为 250小时, 其失效时间区间可精确锁定为相邻监测时间点之间. 例如, 第一个设备的退化值在 3750 至 4000 小时之间超过了阈值. 因此, 可采用线性插值法来估计其失效时间. 第一个设备在 3750 小时时记录的退化值为 9.87, 而在 4000 小时时为 10.94. 通过线性插值法, 其失效时间估计为 \[3750+\dfrac{10-9.87}{10.94-9.87}\times (4000-3750)=3785.75~~ \text{小时}.\] 同理, 第 6 和第 10 个设备的失效时间分别估计为 3506.75 小时和 3351.25 小时.

图 3.3: 每个激光设备的 Kolmogorov-Smirnov 检验 \(p\) 值.
图 3.4: 激光退化数据中 \(\alpha\)、\(\beta_1\)、\(\beta_6\) 和 \(\beta_{10}\) 的在线估计值.

通过 Kolmogorov-Smirnov 检验验证各设备退化数据与伽马过程 \(\mathcal{GP}(\alpha t, \beta_i)\) 的拟合优度. 如图 3.3 所示, 所有设备的检验 p 值均大于 0.05 的显著性水平, 统计上无法拒绝原假设, 表明伽马过程模型能够有效描述退化轨迹的随机特性.

基于上述模型验证结果, 应用第 \(\ref{sec-ITR2025-rul}\) 节提出的在线预测算法进行实时 RUL 预测. 具体实施过程如下: 从 \(t=500\) 小时开始, 仅基于前两次监测数据初始化预测模型; 随着监测数据累积, 通过算法 \(\ref{algo4}\) 递归更新参数后验分布; 计算出每个时间点三台设备的 RUL 点预测值及其 95% 预测区间. 图 3.4 展示了随时间变化的 \(\alpha\)、\(\beta_1\)、\(\beta_6\) 和 \(\beta_{10}\) 的估计值. 从图中可发现, \(\beta_1\)、\(\beta_6\) 和 \(\beta_{10}\) 的估计值呈上升趋势, 且它们之间差异显著, 表面设备间具有较强的异质性. 图 3.5 展示了在各个时间点上三台设备RUL的点预测值 及其 95% 预测区间, 并与每个时间点的实际 RUL 进行了对比. 如图所示, 几乎所有实际 RUL 均被 95% 预测区间覆盖, 且点预测值与实际 RUL 非常接近, 显示了预测的高准确性.

图 3.5: 第 1、6 和 10 个设备的 RUL 点预测值及其 95% 预测区间.

3.2.6.2 列车车轮数据

本节将使用所提出的模型和方法对列车车轮数据 (见图 1.5) 进行分析, 并预测每个测量点上车轮的RUL. 从图 1.5 可以看出, 退化路径呈线性并且单调递增. 有三个车轮在 600 千公里前已经失效. 通过线性插值法计算第5、第9和第11个车轮的失效时间, 结果分别为 523.537、558.861 和 421.508 千公里.

图 3.6: 每个列车车轮的 Kolmogorov-Smirnov 检验 \(p\) 值.

基于 Kolmogorov-Smirnov 检验结果 ( 图 3.6 显示所有车轮退化数据的检验 p 值均大于 0.05) , 验证了伽马过程模型对列车车轮退化数据的适用性. 采用考虑异质性效应的伽马过程模型进行建模, 并从第二次监测时点开始在线 RUL 预测. 如图 3.7 所示, 三个车轮的 RUL 点预测值与实际退化轨迹高度吻合, 且所有时间节点的真实 RUL 均被 95% 预测区间完整覆盖, 表明所提算法在不确定性量化与工程实用性方面具有显著优势, 可为轨道交通关键部件的预防性维护决策提供可靠依据.

图 3.7: 第 5、9 和 11 个车轮的 RUL 点预测值及其 95% 预测区间.

3.2.7 附录

定理 3.1 的证明

基于似然函数 \(\eqref{like0}\) 和先验 \(\eqref{con0}\), \(\alpha\) 和 \(\beta\) 的联合后验密度为 \[\begin{align}\label{post0} \pi(\alpha,\beta|\bm{y})&\propto L(\bm{y} \mid \alpha,\beta)\pi(\alpha,\beta)\nonumber\\ &\propto \dfrac{(\beta\omega_p)^{\delta_p\overline{T}_m\alpha}}{\left[\prod_{j=1}^{m}\left(\Gamma(\alpha t_j)\right)^{1/m}\right]^{\delta_p}} \exp\{-\delta_p\lambda_p\beta\}, \end{align}\] 其中 \[\begin{align*} \delta_p&=mn+\delta,\\ \omega_p&=\omega^{\frac{\delta}{mn+\delta}}\left[\prod_{i=1}^{n}\prod_{j=1}^{m}y_{ij}^{\frac{t_{j}}{nm\overline{T}_m}}\right]^{\frac{mn}{mn+\delta}},\\ \lambda_p&=\frac{mn}{mn+\delta}\bar{y}_a+\frac{\delta}{mn+\delta}\lambda. \end{align*}\] 因此, \(\pi(\alpha,\beta)\) 和 \(\pi(\alpha,\beta \mid \bm{y})\) 均属于同一分布族.

不等式证明

由 \(T_m=\sum_{j=1}^{m}t_j\), 得到 \[\begin{equation*} \begin{aligned} \log\left(\frac{\prod_{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)&= \sum_{j=1}^{m}\frac{t_j}{T_m}\log(t_j)-\log\left(\overline{T}_m\right) \\ &=\frac{\sum_{j=1}^{m}t_j\log(t_j)}{\sum_{j=1}^{m}t_j}-\log\left(\frac{1}{m}\sum_{j=1}^{m}t_j\right) \end{aligned} \end{equation*}\] 令 \(q(x)=x\log(x)\), 则 \(q(x)\) 是凸函数. 利用 Jensen 不等式得 \[\begin{equation*} \begin{aligned} \sum_{j=1}^{m}t_j\log(t_j)&\ge m\cdot \frac{1}{m}\sum_{j=1}^{m}t_j\cdot \log\left(\frac{1}{m}\sum_{j=1}^{m}t_j\right)\\ &=\left(\sum_{j=1}^{m}t_j\right)\cdot \log\left(\frac{1}{m}\sum_{j=1}^{m}t_j\right). \end{aligned} \end{equation*}\] 因此, 可以得出: \(\log\left(\frac{\prod_{j=1}^{m}t_j^{t_j/T_m}}{\overline{T}_m}\right)\ge 0\).

三种后验抽样方法

定理 3.2 的证明

基于似然函数 \(\eqref{like2}\) 和先验 \(\eqref{con2}\), \((\alpha,\beta_1,\dots,\beta_n)^{'}\) 的联合后验密度为 \[\begin{align}\label{post2} \pi(\alpha,\beta_1,\dots,\beta_n \mid \bm{y})&\propto L( \bm{y} \mid \alpha,\beta_1,\dots,\beta_n)\pi(\alpha,\beta_1,\dots,\beta_n)\nonumber\\ & \propto \dfrac{\bar{\beta}_{g}^{(mn+\delta_1)l\alpha}\bar{y}_{g(m)}^{mnl\alpha}\omega^{\delta_1 l\alpha}}{\left[\Gamma(l\alpha)\right]^{mn+\delta_1}}\exp\left\{-\sum_{i=1}^{n}\left(m\bar{y}_{i(m)}+\delta_2\lambda_i\right)\beta_i\right\}\nonumber\\ &\propto \dfrac{\left(\bar{\beta}_{g}\omega_{p(m)}\right)^{\delta_{1p(m)}l\alpha}}{\left[\Gamma(l\alpha)\right]^{\delta_{1p(m)}}} \exp\left\{-\sum_{i=1}^{n}\delta_{2p(m)}\lambda_{ip(m)}\beta_i\right\}, \end{align}\] 其中 \[\begin{align*} \delta_{1p(m)}&=mn+\delta_1, \delta_{2p(m)}=m+\delta_2, \omega_{p(m)}=\omega^{\frac{\delta_1}{\delta_{1p(m)}}}\bar{y}_{g(m)}^{\frac{mn}{\delta_{1p(m)}}},\\ \lambda_{ip(m)}&=\frac{m}{\delta_{2p(m)}}\bar{y}_{i(m)}+\frac{\delta_2}{\delta_{2p(m)}}\lambda_i, i=1,\dots,n. \end{align*}\] 因此, \(\pi(\alpha,\beta_1,\dots,\beta_n)\) 和 \(\pi(\alpha,\beta_1,\dots,\beta_n \mid \bm{y})\) 属于同一分布族.

3.3 重参化伽马过程的变分贝叶斯推断

随着传感器技术的快速发展, 退化数据采集能力的提升与模型复杂度的增加对统计推断方法提出了双重挑战:
一方面需处理高维异构数据, 另一方面需满足实时性计算需求. 传统MCMC方法虽为复杂退化模型提供了严格的贝叶斯推断框架, 但其存在三方面本质性局限(Zaidan 等, 2016):

  • 计算效率的维度诅咒: 即使采用高性能计算集群, 在处理中等规模数据集 ( 如\(n>50\)的异质退化轨迹) 时, 马尔可夫链的混合速度随参数维度呈指数级下降. 以双随机效应伽马模型为例, 单次迭代耗时可达\(O(nm^2)\)量级, 导致实际工程场景中难以获得实时推断结果.
  • 收敛诊断的主观依赖性: 预烧样本长度与抽样间隔的设定高度依赖经验性准则(如Gelman-Rubin统计量或自相关分析). 在退化建模中, 由于参数空间存在强相关性(如形状参数\(\alpha\)与速率参数\(\beta\)), 往往需要超过\(10^4\)次迭代才能达到平稳分布, 显著增加计算资源消耗.
  • 在线更新的结构性障碍: MCMC的批处理特性导致其无法利用增量数据实现递归更新. 当新增监测数据时, 需重新初始化并执行全链采样, 这种”推倒重来”的计算模式严重制约其在PHM系统中的应用.

变分贝叶斯 (Variational Bayesian, VB) 方法通过将贝叶斯推断转化为优化问题, 利用易处理的代理分布近似真实的后验分布, 并以Kullback-Leibler (KL) 散度最小化为目标函数 (Blei 等, 2017; Jaakkola 等, 2000). 相较于MCMC, VB方法在保持精度的同时显著提升计算效率(例如在双随机效应维纳模型中效率提升达87% (Zhou 等, 2020)). 尽管VB方法在复杂模型中展现优势, 但在经典随机效应伽马退化模型中仍面临双重挑战:

  • 参数可解释性缺失: 传统模型通过尺度参数\(\beta\)引入异质性效应, 但\(\beta\)的统计意义(影响退化增量均值\(E[Y(t)]=\alpha t/\beta\))与工程认知(如退化速率、环境应力加速效应)缺乏直接对应, 导致先验分布设定与后验解释存在障碍.
  • 数值稳定性缺陷: 由于伽马分布的共轭结构限制, 传统变分后验需通过数值积分处理随机效应, 这容易出现浮点溢出或精度丢失问题.

针对上述问题, 本节提出重参数化伽马过程 (Reparametrized gamma process, RGa) 及相应的VB算法:

  • 提出了一种RGa模型, 它为模型参数赋予了明确的物理意义, 使得随机效应能够直接且清晰地融入模型, 同时影响退化过程的均值和波动性.
  • 结合高斯-赫尔米特 (Gauss-Hermite, GH) 求积法(Liu 等, 1994)和拉普拉斯近似(Geisser 等, 1990), 介绍了一种具有解析形式的VB推断方法. 同时还基于GH求积法计算的观测信息矩阵, 提出了一种构建区间估计的高效策略, 补充了VB在区间估计方面的局限性.

本节结构如下: 第 3.3.1 节将介绍重参数化伽马过程及产品寿命分布的推导. 第 3.3.2 节将提出基于GH求积法和拉普拉斯近似的VB方法, 并详细介绍区间估计方法. 第 3.3.3 节和第 3.3.4 节将分别进行数值实验和案例研究, 以验证所提方法的有效性.

3.3.1 重参化伽马过程

令 \(\{Y(t), t \geq 0\}\) 表示系统性能的退化过程. 称\(\{Y(t), t \geq 0\}\) 服从RGa过程, 若其满足以下性质: (i) \(Y(0)=0\); (ii) 对于任意的 \(0 \leq t_{i}\), \(i=1,\dots,4\), \(Y(t_{1}) - Y(t_{2})\) 与 \(Y(t_{3}) - Y(t_{4})\) 相互独立; (iii) 每个增量 \(\Delta Y(t) = Y(t)-Y(t+\Delta t)\) 服从具有形状参数 \(\alpha \Delta t\) 和尺度参数 \(\eta/\alpha\) 的伽马分布, 记为 \(\textrm{Ga}(\alpha \Delta t, \eta/\alpha)\), 其中 \(\alpha>0, \eta>0\). 根据伽马分布的可加性, 可以推导出 \(Y(t)\sim \textrm{Ga}(\alpha t, \eta/\alpha)\), 其PDF为 \[\begin{equation}\label{gammpdfx} p_{Y(t)}(y \mid \eta,\alpha) = \frac{1}{\Gamma(\alpha t)}\left(\frac{\eta}{\alpha}\right)^{-\alpha t} y^{\alpha t -1} e^{-\frac{\alpha y}{\eta}}. \end{equation}\] 可算出\(Y(t)\)的均值、方差以及变异系数分别为\(\eta t\), \(\eta^2 t/\alpha\), 和\(1/\sqrt{\alpha t}\). 因此, 参数 \(\eta\) 表示平均退化速率 (或斜率), 而参数 \(\alpha\) 则量化总体样本退化过程的波动程度. 经典伽马过程的一个重要特性是其对退化均值和波动性的依赖性(Tsai 等, 2012). 正如 Ye 等 (2015) 所讨论的, 这种依赖关系在建模时有助于提高拟合精度. 从均值和方差可知, 所提出的RGa过程保持了这一特性, 但参数的解释性更加清晰.

在退化试验中, 由于材料内在差异性和操作环境的多样性, 同批次产品的性能退化路径可能存在差异. 为体现这种差异性, RGa 模型引入随机效应, 假设退化速率参数 \(\eta\) 在不同单元间服从对数正态分布, 其PDF为 \[\begin{equation} \begin{aligned} p(\eta \mid \mu,\sigma^2) =& \frac{1}{\eta \sigma \sqrt{2\pi}} e^{-\frac{(\ln (\eta) -\mu)^2}{2\sigma^2}}, \end{aligned} \label{lognorpdfx} \end{equation}\] 其中位置参数为 \(\mu \in (-\infty, \infty)\), 尺度参数为 \(\sigma > 0\). 参数 \(\eta\) 的均值和方差分别为 \(e^{\mu+\sigma^2/2}\) 和 \(\left(e^{\sigma^2}-1\right)e^{2\mu+\sigma^2}\). 当有关于参数\(\eta\)的先验信息时, 超参数\(\mu\) 和 \(\sigma^2\) 可以通过以下方式确定: \(\mu = \mathbb{E}_{\eta}[\ln(\eta)]\), \(\sigma^2 = \mathbb{E}_{\eta}[(\ln (\eta) - \mu)^2]\) (Krishnamoorthy, 2016). 随机效应采用对数正态分布有如下理由:

  • 对数正态分布的定义域为正实数, 适用于单调递增且为正的退化量. 单调递减的退化过程可通过简单的数学变换处理为单调递增的形式.
  • 对数正态分布属于一类偏态分布, 具有较高的灵活性. 在实例分析中, 该分布对两组真实数据中随机效应变量 \(\eta\) 的经验估计具有良好的拟合效果.
  • 在变分推断框架下, 假设 \(\eta\) 服从对数正态分布时, 参数 \((\mu, \sigma^2)\) 的联合变分后验分布很容易构造. 这一特性放宽了VB推断中参数完全独立的假设, 允许考虑 \((\mu, \sigma^2)\) 之间的相关性. 若 \(\eta\) 服从其他正实数域分布, 则参数 \((\mu, \sigma^2)\) 的联合后验通常难以解析, 增加了推断复杂性.
  • 相比经典伽马过程中通过尺度参数引入个体异质性的方式, RGa 模型中的随机效应设置针对退化速率, 更具物理意义.
  • 所提模型的退化均值仅依赖于 \(\eta\), 而不受形状参数 \(\alpha\) 的影响. 在VB推断中, 假设 \(\alpha\) VB与 \((\mu, \sigma^2)\) 独立有一定的合理性, 间接减少了由近似带来的信息损失.
  • 在固定检测频率 \(\Delta t\) 下, 包含对数正态随机效应的 RGa 模型结构与伽马对数正态模型(Turlapaty, 2020)相似, 其在无线通信系统的建模中有广泛应用.

基于带随机效应的RGa模型, 可得到\(Y(t)\) 的边际PDF (见定理 3.3), 详细数学推导请参见本节附录 3.3.5.

定理 3.3 \(Y(t)\) 的边际PDF为 \[\begin{equation}\label{ungamlnorpdfx2} p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta}) = \frac{\alpha^{\alpha t} y_{t}^{\alpha t -1}}{\Gamma(\alpha t) \sqrt{\pi}} \int_{-\infty}^{+\infty} q(z,y_{t} \mid \boldsymbol{\Theta}) \exp(-z^2) \diff z, \end{equation}\] 其中 \(q(z,y_{t} \mid \boldsymbol{\Theta}) = \exp\left \lbrace {-\alpha t(\mu + z\sqrt{2\sigma^2}) -\alpha y_{t} \exp{[-(\mu +z\sqrt{2\sigma^2})]}}\right \rbrace\).

在式 \(\eqref{ungamlnorpdfx2}\) 中, 积分项难以直接解析计算, 因此需要采用数值方法进行近似. 不同与蒙特卡洛积分通过生成大量随机点来估计积分值, GH求积方法提供了一种高效的近似手段, 其通过选择最优的节点和权重, 能够精确计算特定阶数的多项式积分. 在计算效率和精度方面, GH求积法都远高于蒙特卡洛积分(Seo 等, 2017). 利用GH求积法, 式 \(\eqref{ungamlnorpdfx2}\) 可数值近似为 \[\begin{equation}\label{galnorGH} \int_{-\infty}^{+\infty} q(z,y_{t} \mid \boldsymbol{\Theta}) \exp(-z^2) \diff z \approx \sum_{l=1}^{L} w(z_{l}) q(z_{l},y_{t} \mid \boldsymbol{\Theta}), \end{equation}\] 其中\(z_{l}, l=1,\dots, L\) 是 \(L\) 阶Hermite多项式 \(H_{L}(z_{l})\) 的根 (Davis 等, 2007), 其表达式为 \[H_{L}(z_{l}) = L!\sum_{k=1}^{\lfloor\frac{L}{2}\rfloor}\frac{(-1)^{k}(2z_{l})^{L-2k}}{k!(L-2k)!}, \] 而权重函数 \(w(z_{l})\) (称为 GH 因子) 定义为 \(w(z_{l}) = (2^{n-1}n!\sqrt{\pi})(n^2(H_{L-1}(z_{l})))\). 图 3.8 展示了基于 GH近似的PDF \(\eqref{galnorGH}\) 与四种不同情况下的经验PDF. 结果表明, GH 方法近似的PDF 在所有情况下均表现良好.

图 3.8: 不同设置下的 \(y_{t}\) 的PDF.

接下来对产品寿命分布进行分析. 设退化阈值为 \(\omega\), 产品寿命 \(T\) 的CDF表示为 \[\begin{align}\label{cdfgamlonr} F_{T}(t) = 1 - Pr(Y(t) \leq \omega)= \int_{\omega}^{\infty} p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta}) \diff y_{t}. \end{align}\] 利用 \(p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta})\) 的精确近似, 式 \(\eqref{cdfgamlonr}\) 可通过梯形积分法计算. 在每个检测点\(\tilde{t}\)处, 定义 \(G(y \mid \tilde{t},\boldsymbol{\Theta})=\tilde{p}_{Y(\tilde{t})}(y_{\tilde{t}} \mid \boldsymbol{\Theta})\), 其中 \(\tilde{p}_{Y(\tilde{t})}(y_{\tilde{t}} \mid \boldsymbol{\Theta})\) 是 \(p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta})\) 的GH积分近似. 随后, 使用 R 软件 pracma (Borchers, 2019)包中的 trapz() 函数计算积分 \(\int_{\omega}^{\infty} G(y \mid \tilde{t},\boldsymbol{\Theta}) \diff y\), 该方法与梯形规则积分一致. 为验证积分效果, 图 3.9 比较了基于 GH 的积分框架与不同参数 \(\boldsymbol{\Theta}\) 和退化阈值 \(\omega\) 下寿命 \(T\) 的 BS 分布近似. 结果显示, 两种方法几乎无差异.

图 3.9: 不同设置下寿命 \(T\) 的CDF.

3.3.2 变分贝叶斯推断

在退化试验中, 假设产品的性能特征具有递增的退化路径. 设 \(\mathbf{y}_{n}\) 为第 \(n\) 个样本的观测退化数据, 该样本在检测时间点 \(t_{n,m}\) 处的观测值 \(y_{n,m}\), 其中 \(n=1,\dots,N\), \(m=1,\dots,M_{n}\), \(N\) 表示样本量, \(M_{n}\)为第 \(n\) 个样本的测量次数. 为简化表示, 定义 \(\mathbf{y} = \{\mathbf{y}_{n}\}_{n=1}^{N}\), 其中 \(\mathbf{y}_{n}=\{y_{n,m}\}_{m=1}^{M_{n}}\). 假设 \(\mathbf{y}\) 是来自RGa 过程 \(\{Y(t), t \geq 0\}\) 的一组观测值. 定义增量 \(x_{n,m} = \Delta y_{n,m}=y_{n,m}-y_{n,m-1}\), 时间间隔 \(h_{nm} = \Delta t_{n,m}=t_{n,m}-t{n,m-1}\), 其中\(y_{n,0}=t_{n,0}=0\). 记\(\mathbf{x}_{n} = \{x_{n,m}\}_{m=1}^{M_{n}}\), \(\mathbf{x} = \{\mathbf{x}_{n}\}_{n=1}^{N}\). 为符号方便, 假设\(M_{n} = M\), \(h_{n,m}=h_{c}\), 即, 试验采用等间隔测量, 且所有样本的测量次数相同. 根据RGa过程的性质 (iii), 增量\(\mathbf{x}\)相互独立, 且服从伽马分布. 在给定随机效应 \(\boldsymbol{\eta} = \{\eta_{n}\}_{n=1}^{N}\) 时, 增量 \(x_{n,m}\) 的条件PDF为 \[\begin{equation} \begin{aligned} p_{x \mid \eta_{n}}(x_{n,m} \mid \alpha h_{c}, \eta_{n}) = \frac{x_{n,m}^{\alpha h_{c}-1}}{\Gamma (\alpha h_{c})} \left(\frac{\eta_{n}}{\alpha}\right)^{-\alpha h_{c}}%\\\nonumber e^{-\frac{ \alpha x_{n,m}}{\eta_{n}}}. \end{aligned}\label{gammpdf} \end{equation}\] 随机效应 \(\boldsymbol{\eta}\)为独立同分布的对数正态分布, 其PDF为 \[\begin{equation}\label{lognorpdf} p(\boldsymbol{\eta} \mid \mu,\sigma^2)=\prod_{n=1}^{N} p(\eta_{n} \mid \mu,\sigma^2) = \prod_{n=1}^{N}\frac{1}{\eta_{n} \sigma \sqrt{2\pi}} e^{-\frac{(\ln(\eta_{n}) -\mu)^2}{2\sigma^2}}. \end{equation}\] 令 \(\boldsymbol{\Theta}=(\alpha, \mu, \sigma^2)\), 则 \((\boldsymbol{\eta}, \boldsymbol{\Theta})\) 的联合后验分布为 \[\begin{equation} \label{gamjointpdf} p(\boldsymbol{\eta}, \boldsymbol{\Theta} \mid \mathbf{x}) \propto \prod_{n=1}^{N}\left[p(\eta_{n} \mid \mu,\sigma^2)\prod_{m=1}^{M}p_{x \mid \eta}(x_{n,m} \mid \alpha h_{c}, \eta_{n})\right] p(\alpha)p(\mu,\sigma^2), \end{equation}\] 其中 \(p(\alpha)\) 和 \(p(\mu,\sigma^2)\) 分别是参数 \(\alpha\) 和 \((\mu,\sigma^2)\) 的先验分布. 式 \(\eqref{gamjointpdf}\) 中的先验设置隐含了以下假设: (i) 允许参数 \(\mu\) 和 \(\sigma^2\) 存在一定的相关性, 可选用二元分布作为其先验; (ii) 假设 \(\alpha\) 与 \((\mu, \sigma^2)\) 相互独立. 这一假设基于上一节中关于随机效应设置的讨论, 认为 \(\alpha\) 主要描述形状参数的变化, 与刻画单元间异质性的参数 \((\mu, \sigma^2)\) 无关.

3.3.2.1 变分贝叶斯估计

由于式 \(\eqref{gamjointpdf}\) 中参数的后验分布非常复杂, 直接进行统计推断较为困难, VB方法通过引入一组可处理的变分后验密度 (Variational posterior density, VPD) \(q(\boldsymbol{\Theta},\boldsymbol{\eta})\) 来近似真实的后验密度 \(p(\boldsymbol{\eta}, \boldsymbol{\Theta} \mid \mathbf{x})\), 即 \(p(\boldsymbol{\eta}, \boldsymbol{\Theta} \mid \mathbf{x}) \approx q(\boldsymbol{\eta}, \boldsymbol{\Theta})\), 该近似可通过求解以下约束优化问题确定. \[\begin{equation}\label{optim} q^{*}(\boldsymbol{\Theta},\boldsymbol{\eta}) = \mathop{\textrm{arg max}}\limits_{q(\boldsymbol{\Theta},\boldsymbol{\eta}) \in \mathcal{Q}} -F(q(\boldsymbol{\Theta},\boldsymbol{\eta})), s.t. \int q(\boldsymbol{\Theta},\boldsymbol{\eta}) d\boldsymbol{\Theta} d\boldsymbol{\eta} =1, \end{equation}\] 其中 \(F(q(\boldsymbol{\eta}, \boldsymbol{\Theta}))= E_{\boldsymbol{\Theta},\boldsymbol{\eta}}\log\left( \frac{q(\boldsymbol{\eta}, \boldsymbol{\Theta})}{p(\boldsymbol{\eta}, \boldsymbol{\Theta} \mid \mathbf{x})} \right)\), 即, \(q(\boldsymbol{\eta}, \boldsymbol{\Theta})\)和 \(p(\boldsymbol{\eta}, \boldsymbol{\Theta} \mid \mathbf{x})\)的Kullback-Liebler(KL) 散度. \(\mathcal{Q}\) 表示所有潜在VPD的集合. 通常而言, 集合 \(\mathcal{Q}\) 的选择直接决定了 VPD 的复杂性, 从而影响优化问题 \(\eqref{optim}\) 的求解难度, 并进一步影响 VB 推断的计算效率. 因此, 合理选择 \(\mathcal{Q}\) 至关重要, 它应足够简单, 以确保优化问题 \(\eqref{optim}\) 可行, 同时又要具备足够的灵活性, 以准确逼近真实的后验分布. 一种常用的选择是均值场变分族 (Mean-field variational family, MFVF), 该方法假设模型参数之间相互独立, 或者在一定条件下表现出独立性. 在 MFVF下, 可以进一步假设:
\[\begin{equation}\label{gamfvf} q(\boldsymbol{\Theta},\boldsymbol{\eta}) = \prod_{n=1}^{N} q(\eta_{n}) q(\alpha) q(\mu,\sigma^2). \end{equation}\] 其中 \(q(\eta_{n})\)、\(q(\alpha)\) 和 \(q(\mu,\sigma^2)\) 分别是 \(\eta_{n}\)、\(\alpha\) 和 \((\mu,\sigma^2)\) 的VPD. 该分解形式在降低计算复杂度的同时, 仍能有效逼近真实的后验分布. 将 MFVF \(\eqref{gamfvf}\) 代入优化目标 \(\eqref{optim}\) 中, 然后使用拉格朗日乘数法求解式\(\eqref{optim}\), 得到 \[\begin{equation}\label{new_lg} F_{1}(q(\boldsymbol{\Theta},\boldsymbol{\eta}))=F(q(\boldsymbol{\Theta},\boldsymbol{\eta}))+\boldsymbol{\lambda}^{T} \left(\int (q(\boldsymbol{\Theta},\boldsymbol{\eta})-1) d\boldsymbol{\Theta} d\boldsymbol{\eta}\right), \end{equation}\] 其中向量 \(\boldsymbol{\lambda}=(\left\{\lambda_{\eta_{n}}\right\}_{n=1}^{N},\lambda_{\alpha},\lambda_{\mu},\lambda_{\sigma^2})\) 是拉格朗日乘数, 表示 \(q(\boldsymbol{\Theta},\boldsymbol{\eta})\) 对目标函数的边际效应.

VB算法可视为传统EM方法和贝叶斯推断在计算时间和计算效率上的一种折中, 通过使用易处理的代理分布来近似复杂的贝叶斯后验. 它利用KL散度将贝叶斯推断转化为数值优化问题, 从随机模拟转向确定性优化. 如 Blei 等 (2017) 所述, VB的优化分为两个迭代步骤: 针对潜在变量的VB后验更新 (变分贝叶斯期望 (VB-E) 步骤) 和针对目标参数的VB后验更新 (变分贝叶斯最大化 (VB-M) 步骤) , 与EM算法的迭代过程相似.

3.3.2.1.1 VB-E: 最优潜变量

对式 \(\eqref{new_lg}\) 中关于 \(q(\boldsymbol{\eta})\) 求导, 可以得到 \(\boldsymbol{\eta}\) 的VPD的对数形式 \[\begin{equation*} \begin{aligned} \ln(q(\boldsymbol{\eta})) =& \sum_{n=1}^{N}\mathbb{E}_{\mu,\sigma^2}\left[\ln(p_{\eta}(\eta_{n} \mid \mu,\sigma^2))\right] + \sum_{n=1}^{N}\sum_{m=1}^{M}\mathbb{E}_{\alpha}\left[p_{x \mid \eta}(x_{n,m} \mid \alpha h_{c}, \eta_{n})\right] + C, \end{aligned} \end{equation*}\] 其中\(C\)表示常数. 为简便符号, 本节后续内容涉及到常数都用\(C\)来表示, 这不影响 各参数的VPD形式. 由独立性假设, \(\eta_{n}\) 的VPD的对数形式为 \[\begin{equation*} \begin{aligned} \ln(q(\eta_{n})) =& \mathbb{E}_{\mu,\sigma^2}\left[\ln(p_{\eta}(\eta_{n} \mid \mu,\sigma^2))\right] + \sum_{m=1}^{M}\mathbb{E}_{q(\alpha)}\left[p_{x \mid \eta}(x_{n,m} \mid \alpha h_{c}, \eta_{n})\right] + C. \end{aligned} \end{equation*}\] 利用式 \(\eqref{gammpdf}\) 和 \(\eqref{lognorpdf}\), 可以得到 \[\begin{equation} \begin{aligned} \ln(q(\eta_{n})) = &(w_{1}w_{2} - w_{0}Mh_{c}-1) \ln(\eta_{n}) - \frac{y_{n,M}w_{0}}{\eta_{n}} - w_{2}(\ln(\eta_{n}))^2/2 +C, \end{aligned}\label{vpdeta3} \end{equation}\] 其中\(y_{n,M} = \sum_{m=1}^{M} x_{n,m}\), \(w_{0} = \mathbb{E}_{\alpha}\left[\alpha\right]\), \(w_{1}=\mathbb{E}_{\mu,\sigma^2}\left[\mu\right]\), 以及 \(w_{2}=\mathbb{E}_{\mu,\sigma^2}\left[\sigma^{-2}\right]\). 则\(\eta_{n}\) 的VPD可表示为 \(q(\eta_{n}) = g_{0}(\eta_{n})/A_{\eta_{n}}\), 其中 \[ g_{0}(\eta_{n}) = \eta_{n}^{w_{1}w_{2} - w_{0}Mh_{c}-1} e^{- \frac{y_{n,M}w_{0}}{\eta_{n}}- \frac{w_{2}(\ln(\eta_{n}))^2}{2}}, \] \(A_{\eta_{n}} = \int_{0}^{+\infty} g_{0}(\eta_{n}) \diff\eta_{n}\)为正则化常数. 令 \(s = \ln(\eta_{n})\sqrt{w_{2}/2}\). 将 \(s\) 代入 \(g_{0}(\eta_{n})\) 得到 \[\begin{equation}\label{g0etan} g_{0}(\eta_{n}) = q(s,y_{n,M} \mid w_{0},w_{1},w_{2})e^{-s^2}, \end{equation}\] 其中 \[\begin{equation*} \begin{aligned} q(s,y_{n,M}& \mid w_{0},w_{1},w_{2}) =\sqrt{2/w_{2}}e^{s\sqrt{2/w_{2}}(w_{1}w_{2}-w_{0}Mh_{c}) -e^{-y_{n,M}w_{0}s\sqrt{2/w_{2}}}}. \end{aligned} \end{equation*}\] 注意到 \(g_{0}(\eta_{n})\) 与式 \(\eqref{galnorGH}\) 中的被积函数具有相同的结构. 因此, \(A_{\eta_{n}}\) 也可以通过第 3.3.1 节中介绍的GH求积方法来计算. 基于式 \(\eqref{vpdeta3}\) 和 \(\eqref{g0etan}\), 任何连续函数 \(f(\eta_{n})\) 的期望可以表示为 \[\begin{align} \mathbb{E}_{\eta_{n}}(f(\eta_{n})) = & \int_{0}^{+\infty} f(\eta_{n}) q(\eta_{n}) \diff \eta_{n}\nonumber\\ =&\frac{1}{A_{\eta_{n}}}\int_{-\infty}^{+\infty} f(\exp(s\sqrt{2/w_{2}}))e^{-s^2} q(s,y_{n,M} \mid w_{0},w_{1},w_{2})\diff s,\label{meta} \end{align}\] 这一积分也可通过GH求积法进行数值计算.

3.3.2.1.2 VB-M: 最优模型参数

对式 \(\eqref{new_lg}\) 中关于 \(q(\boldsymbol{\Theta})\) 求导, 可以得到 \(\boldsymbol{\Theta}\) 的VPD的对数形式为 \[\begin{equation}\label{vpdotr} \begin{split} \ln(q(\alpha,\mu,\sigma^2)) =& \sum_{n=1}^{N}\left[E_{\eta_{n}}\left[\ln(p(\eta_{n} \mid \mu,\sigma^2))\right]\right] \nonumber\\ &+ \sum_{n=1}^{N}\sum_{m=1}^{M}\left[E_{q(\eta_{n})}\left[\ln(p(x_{n,m} \mid \eta_{n},\alpha))\right]\right] \nonumber\\ &+ \ln(p(\alpha)) + \ln(p(\mu,\sigma^2)) +C. \end{split} \end{equation}\] 可发现 \(\alpha\) 和 \((\mu,\sigma^2)\)是相互独立的. 因此, 给定 \(q(\eta_{n}), n=1,2,\dots, N\)后, 可以依次推导出 \(q(\alpha)\) 和 \(q(\mu,\sigma^2)\) 的VPD表达式.

  1. \(\alpha\) 的VPD: \(\alpha\) 的VPD的对数形式为 \[\begin{equation} \begin{aligned} \ln(q(\alpha)) =& \sum_{n=1}^{N}\sum_{m=1}^{M}\left[E_{q(\eta_{n})}\left[\ln(p(x_{n,m} \mid \eta_{n},\alpha))\right]\right] + \ln(p(\alpha)) + C. \end{aligned}\label{vpdalp0} \end{equation}\] 再将 \(\eqref{gammpdf}\) 和 \(\eqref{lognorpdf}\) 代入 \(\eqref{vpdalp0}\), 得到 \[\begin{equation} \begin{aligned} \ln(q(\alpha)) = & - (h_{c}(Mw_{3}-w_{5}) + w_{4}) \alpha + \alpha h_{c} NM \ln(\alpha)\nonumber\\ &- NM \ln(\Gamma(\alpha h_{c})) + \ln(p(\alpha)) + C, \end{aligned}\label{vpdalp1} \end{equation}\] 其中 \[ w_{3}=\sum_{n=1}^{N} E_{\eta_{n}}\left[\ln(\eta_{n})\right], w_{4}= \sum_{n=1}^{N}y_{nM}E_{\eta_{n}}\left[\eta_{n}^{-1}\right], w_{5}=\sum_{n=1}^{N}\sum_{m=1}^{M}\log(x_{nm}), \] 并且 \(E_{\eta_{n}}\left[\ln(\eta_{n})\right]\) 和 \(E_{\eta_{n}}\left[\eta_{n}^{-1}\right]\) 可以通过式 \(\eqref{meta}\)更新. 对\(\eqref{vpdalp1}\) 两边取指数变换, 可得到 \[\begin{equation}\label{vpdalp2} q(\alpha) \propto p(\alpha) \frac{\alpha^{a\alpha}}{\Gamma(\alpha h_{c})^b}\exp\left(-c\alpha\right), \end{equation}\] 其中 \(a = h_{c}NM\), \(b=NM\), \(c=h_{c}(Mw_{3}-w_{5})+w_{4}\). 为方便数学推导, 使用共轭先验可简化计算(Gelman 等, 2014). 可验证: 当先验\(p(\alpha)\)取以下形式时, 其与 \(q(\alpha)\)属于相同的分布族:
    \[\begin{equation}\label{vpdalppri} p(\alpha) \propto \frac{\alpha^{a_{0}}}{\Gamma(\alpha h_{c})^{b_{0}}}\exp(-c_{0}\alpha), \end{equation}\] 其中 \(a_{0}\)、\(b_{0}\) 和 \(c_{0}\) 是先验 \(p(\alpha)\) 的超参数. 将 \(\eqref{vpdalppri}\) 代入 \(\eqref{vpdalp2}\) 并计算正则化常数, 可得 \[\begin{equation}\label{vpdalppos} q(\alpha) = g_{\alpha}(\alpha)/A_{\alpha} \end{equation}\] 其中 \[ g_{\alpha}(\alpha) = \frac{\alpha^{a_{1}\alpha}}{\Gamma(\alpha h_{c})^{b_{1}}}\exp(-c_{1}\alpha),~~A_{\alpha} = \int_{\Omega_{\alpha}} g_{\alpha}(\alpha) \diff \alpha, \] \(a_{1}=a_{0}+a\), \(b_{1}=b_{0}+b\), \(c_{1}=c_{0}+c\). 然而, 由于 \(g_{\alpha}(\alpha)\) 的积分可能在数值上不稳定, 通过 \(\eqref{vpdalppos}\) 计算 \(\alpha\) 的VPD的均值和方差比较困难. 作为替代方案, 可对 \(\eqref{vpdalppos}\) 进行拉普拉斯近似. 此时, \(\alpha\) 的VPD的均值和方差可近似为 \[\begin{equation*} \begin{aligned} m(\alpha) &\approx \arg \min_{\alpha} \left\{a_{1}\log(\alpha)-b_{1} h_{c} \psi(\alpha h_{c}) + a_{1}-c_{1}\right\},\nonumber\\ v(\alpha) &\approx \frac{1}{bh_{c}^{2}\psi'(m(\alpha) h_{c})- \frac{a_{1}}{m(\alpha)}}. \end{aligned} \end{equation*}\]

  2. 参数 \((\mu,\sigma^2)\)的联合VPD: 基于 \(\eqref{vpdotr}\), \((\mu, \sigma^2)\) 的联合VPD的对数形式为 \[\begin{equation} \begin{aligned}\label{vpdmu0} \ln(q(\mu,\sigma^2)) =& \sum_{n=1}^{N}\left[E_{q(\eta_{n})}\left[\ln(p(\eta_{n}|\mu,\sigma^2))\right]\right] + \ln(p(\mu,\sigma^2)) +C. \end{aligned} \end{equation}\] 将 \(\eqref{gammpdf}\) 和\(\eqref{lognorpdf}\)代入 \(\eqref{vpdmu0}\) 并对结果取指数, 可得 \[\begin{equation} \begin{aligned} q(\mu,\sigma^2) \propto p(\mu,\sigma^2)(\sigma^{2})^{-N/2}%\nonumber\\ e^{- \frac{w_{6}-w_{3}^{2}/N}{2\sigma^2} - \frac{N(\mu-w_{3}/N)^2}{2\sigma^2}}, \end{aligned}\label{vpdmu1} \end{equation}\] 其中 \[w_{6} = \sum_{n=1}^{N} E_{q(\eta_{n})}\left[(\ln(\eta_{n}))^2\right], E_{\eta_{n}}\left[(\ln(\eta_{n}))^2\right] \] 可以通过式 \(\eqref{meta}\) 计算. 若 \(p(\mu,\sigma^2)\) 服从正态逆伽马 (Normal inverse gamma, NIGa) 分布, 式 \(\eqref{vpdmu1}\) 恰好是NIGa密度函数的核 (Bolstad 等, 2016). 通常, 随机变量 \((X,Y)\) 服从 \(\textrm{NIGa}((x,y); \rho, \varsigma, \nu, \xi)\), 其PDF为 : \[\begin{equation}\label{niga} p(x,y) = \frac{\sqrt{\varsigma}}{\sqrt{y}\sqrt{2\pi}}\frac{\xi^{\nu}}{\Gamma(\nu)}y^{-\nu-1}e^{-\frac{2\xi + \varsigma (x-\rho)^2}{2y}}, \end{equation}\] 其中 \(x \in (-\infty, +\infty)\), \(y \in (0, +\infty)\), \(\nu>0\), \(\xi>0\), \(\varsigma>0\), 以及 \(\rho\in(-\infty, +\infty)\) 是参数. 根据NIGa分布的定义 (Denison 等, 2002), \(\eqref{niga}\) 可分解为在给定 \(\sigma^2\) 条件下, \(\mu\) 的均值为 \(\rho\)、方差为 \(\sigma^2/\varsigma\) 的条件正态分布和形状参数为 \(\nu\)、尺度参数为 \(\xi\) 的 \(\sigma^2\) 的逆伽马分布. 因此, 可假设参数 \((\mu,\sigma^2)\) 的联合先验为 \[\begin{equation}\label{vpdmupri} p(\mu,\sigma^2) = \textrm{NIGa}((\mu,\sigma^2); \rho_{0}, \varsigma_{0}, \nu_{0}, \xi_{0}). \end{equation}\] 将式 \(\eqref{vpdmupri}\) 代入 \(\eqref{vpdmu1}\) 后, 可得到 \((\mu,\sigma^2)\) 的联合VPD为 \[\begin{equation}\label{vpdmupos} q(\mu,\sigma^2) = \textrm{NIGa}((\mu,\sigma^2); \rho_{1}, \varsigma_{1}, \nu_{1}, \xi_{1}), \end{equation}\] 其中 \[ \begin{aligned} \rho_{1}&=\frac{1}{\varsigma_{1}}\left[\varsigma_{0}\rho_{0}+w_{3}\right], \varsigma_{1} = \varsigma_{0} + N, \nu_{1} = \nu_{0} + \frac{N}{2},\\ \xi_{1} &= \xi_{0} + \frac{1}{2}\left[\varsigma_{0}\rho_{0}^2+w_{6}-\varsigma_{1}\rho_{1}^2\right]. \end{aligned} \] 根据NIGa分布的性质 (Bolstad 等, 2016), \(\mu\) 和 \(\sigma^2\) 的VPD均值和方差为 \[\begin{equation*} \begin{aligned} m(\mu) =&\rho_{1},~~~~v(\mu) = \frac{\xi_{1}}{\varsigma_{1}(\nu_{1}-1)},\nonumber\\ m(\sigma^2) =& \frac{\xi_{1}}{\nu_{1}-1},~~~~ v(\sigma^2) = \frac{\xi_{1}^2}{(\nu_{1}-1)^2(\nu_{1}-2)}, \end{aligned} \end{equation*}\] 并且 \(\sigma^{2}\) 的一阶逆矩等于 \(\frac{\nu_{1}}{\xi_{1}}\).

通过当前模型的VB-E和VB-M步的更新方程, 可对这两步进行迭代, 直到估计方差不再显著改善, 即 \[\begin{equation}\label{differror} \varepsilon(\boldsymbol{\Theta^{(s)}})= ||v(\boldsymbol{\Theta}^{(s)})-v(\boldsymbol{\Theta}^{(s-1)})||_{1}\leq \varepsilon_{0}, \end{equation}\] 其中 \(||\cdot||_{1}\) 为绝对值范数, \(\varepsilon_{0}\) 是确定所需精度的预设值.

注5: VB算法是否成功与高效, 关键在于均值场变分族的假设和模型参数的共轭特性. 均值场变分族假设使得能够基于KL散度和拉格朗日优化算法求得近似贝叶斯后验分布. 共轭参数特性则提供了解析形式的近似贝叶斯后验. 对于非共轭模型参数, 也可以基于自然共轭先验通过拉普拉斯近似得到相应的变分后验. 这不仅减少了数值优化需求, 提高了计算效率, 而且在VB和贝叶斯后验中保持了结构的相似性, 从而提升了精度. 然而, 当模型参数如 \(\eta\) 服从非共轭偏态分布时, 非共轭参数的增加则会降低模型的计算效率和精度.

3.3.2.2 初始值确定

在VB算法中, 参数初始值的选择会影响算法的收敛时间及其局部最优解. 在本节中, 参数 \(\bm{\Theta}\) 的初始值将通过以下步骤确定:

  1. 对于第\(n\)个样品, 分别计算增量 \(\bm{X}_{n}\) 的均值和方差为 \(\bar{x}_{n}=\frac{1}{M} \sum_{m=1}^{M} x_{n,m}\) 和 \(V_{n} = \frac{1}{M-1} \sum_{m=1}^{M} (x_{n,m}-\bar{x}_{n})^2\).

  2. 由于增量 \(\bm{X}_{n}\) 服从 \(Ga(\alpha h_{c}, \eta_{n}/\alpha)\), 其均值为 \(\eta_{n} h_{c}\), 方差为 \(\eta_{n}^2 h_{c}/\alpha\), 因此可估计 \(\eta_{n}\) 和 \(\alpha\) 为 \(\bar{x}_{n}/h_{c}\) 和 \(\bar{x}_{n}^2/(V_{n} h_{c})\).

  3. 由 \(\{\eta_{n}\}_{n=1}^{N}\), 参数 \(\mu\) 和 \(\sigma^2\) 的初始值可以通过以下公式确定: \(\bar{\eta}_{n}=\frac{1}{N}\sum_{n=1}^{N} \ln(\eta_{n})\) 和 \(V_{\ln(\eta_{n})} = \frac{1}{N-1}\sum_{n=1}^{N}(\ln(\eta_{n}) - \bar{\eta}_{n})^2\).

3.3.2.3 置信区间

通过上述VB推断方法, 可得到每个 \(\theta_{i}\in\boldsymbol{\Theta}\) 的VPD均值 \(m(\theta_{i})\) 和方差 \(v(\theta_{i})\). 然后构造 \(\theta_{i}\) 的\(100(1-\gamma)%\) 近似区间估计为 \[\begin{equation}\label{ci1} m(\theta_{i}) \pm Z_{1-\gamma/2} \sqrt{v(\theta_{i})}, \end{equation}\] 其中 \(Z_{\gamma}\) 为标准正态分布的\(\gamma\)分位数. 然而, 式 \(\eqref{ci1}\) 存在不足, 因为VB方法通常会低估后验分布的散度, 从而导致方差估计值比真实值小 (Blei 等, 2017). 为了解决这个问题, Wang 等 (2005) 对VB推断方法下的区间估计进行改进, 其中VPD的方差可通过原始似然函数下的Fisher信息矩阵计算. 当Fisher信息矩阵计算困难时, 则可进一步用观测Fisher信息矩阵替代 (Meeker 等, 1998). 基于此思想, 当前模型的观测Fisher信息矩阵为 \[\begin{equation}\label{ci2} \hat{I}_{\boldsymbol{\Theta}} = -\frac{\partial^2 \log\left(\ell(\boldsymbol{\Theta})\right)}{\partial \boldsymbol{\Theta}^{T} \partial\boldsymbol{\Theta}}, \end{equation}\] 其中 \[\begin{equation*} \begin{aligned} \ell(\boldsymbol{\Theta}) =& - NM \log(\Gamma(\alpha h_{c})) - NM \alpha h_{c} \log(\alpha) \nonumber\\ & + \sum_{n=1}^{N}\sum_{m=1}^{M} (\alpha h_{c}-1)\log(x_{n,m})-\frac{N}{2}\log(\pi)\nonumber\\ & + \sum_{n=1}^{N}\log\left[\int_{-\infty}^{\infty} q(z,y_{n,M},\boldsymbol{\Theta})\exp(-z^2)\diff z\right], \end{aligned} \end{equation*}\] 其中 \(z=\frac{\log(\eta_{n}-\mu)}{\sqrt{2\sigma^2}}\), \[\begin{align}\label{ci4} q(z,y_{n,M},\boldsymbol{\Theta}) = e^{-\alpha h_{c} M\left(\sqrt{2\sigma^2}z+\mu\right) -\alpha y_{n,M} e^{-\sqrt{2\sigma^2}z-\mu}}, \end{align}\] 详细导数可见本节附录 3.3.5. 本节附录 3.3.5 还给出其他两种基准估计方法, 分别为EM算法和贝叶斯方法. 后续会在数值模拟和案例研究中与所提VB推断方法进行比较.

3.3.3 模拟实验

模拟实验将验证 VB 算法的计算精度, 同时比较VB、EM 和 MCMC 三种算法的计算效率. 三种方法的基本设置如下:

  1. VB 算法: 先验分布 \(p(\alpha)\) 和 \(p(\mu, \sigma^2)\) 的超参数使用弱信息先验设定, 具体为 \(a_0 = b_0 = c_0 = 1\), \(\xi_0 = \varsigma_0 = \rho_0 = 1\), \(\nu_0 = 0\), 以便与 EM 算法进行公平比较.
  2. EM 算法: 采用最近两次迭代估计值的绝对值范数作为收敛停止准则.
  3. MCMC 算法: 使用与 VB 相同的超参数设定, 后验采样中包含 3 条马尔科夫链, 每条链迭代 20,000 次, 前半部分作为燃烧期.
  4. 通用设置: VB 和 EM 算法的目标精度为 \(10^{-6}\). 所有算法均通过 R 软件实现, 运行于 Windows 10 桌面电脑 (3.5 GHz Intel Xeon E3-1241 v3 处理器, 16 GB 内存) .

样本量 \(N\) 选为15、30和60, 对于每个\(N\), 测量次数 \(M=\) 10和20. 对于每种 \((N, M)\) 组合, 从参数 \(\boldsymbol{\Theta}=(\alpha, \mu, \sigma^2)=(10, -1, 0.01)\) 的RGa过程中随机生成2,000个退化样本 \(\boldsymbol{y}\). 对每个样本 \(\boldsymbol{y}\), 通过VB、EM和MCMC算法获得 \(\boldsymbol{\Theta}\) 的估计值. 因此, 每种算法均得到2,000个估计值. 基于这些估计值, 计算平均值和MSE. 结果列在表 \(\ref{tbl-tab:ns1}\) 中. 同时还计算了每种算法的平均CPU运行时间(秒, \(\bar{t}\)), 结果如图 3.10 所示. 此外, 还计算了95%可信区间的覆盖率和区间长度 (Interval length, IL), 见表 \(\ref{tbl-tab:ns2}\). 从这两个表格和三种算法的平均CPU运行时间图中, 可以得出以下结果: (1) 所有情况下, 估计值的平均值接近真实值, 且样本量增加时, 平均值与真实值的差异逐渐缩小. (2) 随着样本量增加, 各模拟设置下的 MSE 总体呈下降趋势. VB 和 EM 的 MSE 小于贝叶斯方法, 且样本量增加时, 三者差距逐渐减小. (3) 所有算法的IL随样本量增加而缩小, 同时 基于VB 与贝叶斯方法的 IL 差异也逐渐减小. (4) 在覆盖率方面, VB 略优于 EM, 表现出更接近 0.95 的名义水平, 而贝叶斯方法在小样本 (如 \(N = 15, M = 10, 20\)) 时表现最佳. (5) 在计算效率上, VB 显著优于 EM 和贝叶斯方法, 部分情况下运行时间仅为贝叶斯方法的百分之一或千分之一. 总体而言, 贝叶斯方法提供最高的估计精度, 但计算效率低; 三种算法的估计精度随样本量增加而提高; EM 提高了效率, 但牺牲了部分精度; VB 在精度与效率之间实现了良好平衡, 是贝叶斯与 EM 方法的折中方案.

图 3.10: 在不同算法和样本量下的平均CPU运行时间.

3.3.4 实例分析

3.3.4.1 激光退化数据

以激光退化数据为例 (见图 1.2 ) 进行分析. 拟合前需验证两个目标: 一是检查退化率的单元间变异性是否符合对数正态分布, 二是研究激光数据中 RGa 过程的形状参数与尺度参数的相关性. 为验证第一个目标, 计算经验估计量 \(\hat{\boldsymbol{\eta}} = \{{\hat{\eta}_n = y_{nM} / t_{nM}}\}_{n=1}^{15}\), 并利用式 \(\eqref{lognorpdfx}\) 估计对数正态分布的参数, 得到 \(\hat{\boldsymbol{\eta}} \sim LN(-0.697, 0.045)\). 图 3.11 (a) 展示了带 95% CI 的经验分布函数 和估计的对数正态分布, 表明退化率 \(\hat{\boldsymbol{\eta}}\) 符合对数正态分布的特征. 此外, 通过Anderson-Darling 检验得到的P值为0.633, 进一步说明使用对数正态分布的合理性. 在研究RGa过程形状参数与尺度参数的相关性时, 定义经验估计量为 \[ \hat{\boldsymbol{\alpha}},\hat{\boldsymbol{\beta}}=\left\{\hat{\alpha}_{n}=\frac{\hat{\eta}_{n}^2}{v(\boldsymbol{x}_{n}/\sqrt{h_{c}})}, \hat{\beta}_{n}=\frac{\hat{\eta}_{n}}{\hat{\alpha}_{n}} \right\}_{n=1}^{15}. \] 图 3.11 (b) 展示了 \(\hat{\boldsymbol{\alpha}}^{-1}\) 和 \(\hat{\boldsymbol{\beta}}\) 的散点分布情况以及 回归曲线, 相关系数 \(\rho=0.958\) 表明两者之间存在显著正相关性.

(a) 估计CDF
(b) 形状参数与尺度参数的关系
图 3.11: 模型验证相关结果.

接下来, 采用 VB 算法对 RGa 退化模型进行拟合, 并将其应用于激光数据集, 同时并与 EM 和贝叶斯方法进行比较, 以衡量计算效率和拟合精度. 表 \(\ref{tbl-tab:addlabel}\) 展示了各算法的参数估计值、95% 区间估计以及运行时间. 结果显示, 不同算法的参数估计值相近, 但 VB 算法在 CPU 运行时间上显著优于 EM 算法, 且远低于贝叶斯算法. 图 3.12 (a) 显示了样本平均值、估计的平均退化路径 (\(\hat{\eta} t\)) 以及各算法的 95% 区间估计, 表明所有算法均能有效地拟合样本平均值.

本文提出的模型还与五个基准模型进行了比较,包括伽马过程、逆高斯过程、维纳过程、具有随机效应的同质伽马过程和 Tsai 等 (2012) 所提的RGa过程. 为简化表述, 将本文提出的随机效应伽马过程模型记为 RGaLn, Tsai 等 (2012) 的模型记为 GaGa. 鉴于 EM 和 VB 算法在估计结果上无显著差异, 模型选择基于 EM 算法的结果. 同时 VB 算法作为一种确定性近似方法, 不适用于计算涉及有效样本数量的 BIC, 因此使用 AIC 准则进行评估. 表 \(\ref{tbl-tab:addlabe2}\) 列出了各模型的 AIC 值, 显示随机效应模型 (GaGa 和 RGaLn) 的 AIC 值低于其他模型, 表明在激光退化数据中考虑随机效应是合理的. 然而, RGaLn 和 GaGa 的 AIC 值差异不显著. 图 3.12 (b) 展示了经验 分布函数, 以及阈值 \(\omega=6\) 下 基于GaGa 和 RGaLn 模型的产品寿命分布. 经验分布函数通过中位数秩估计对伪失效时间进行非参数估计(Meeker 等, 1998). 可以看出, GaGa 和 RGaLn 模型均能很好地拟合经验分布函数.

(a) 平均退化路径估计
(b) 寿命估计CDF
图 3.12: 激光数据下寿命分布相关结果.

3.3.4.2 设备-B数据

这里以 7 个设备-B 的退化数据为例 (见图 1.6 中\(150^{\circ}\)C条件下的7个设备-B样本) 进行分析. 图 3.13 (a) 显示了退化率 \(\boldsymbol{\hat{\eta}}\) 的经验分布函数 及其 95% 置信区间, 以及拟合的对数正态分布. 结果表明, 对数正态分布能够很好地拟合退化率分布. 同时Anderson-Darling 检验的 \(P\) 值为 0.883, 进一步验证用对数正态分布刻画样品间变异性的合理性. 图 3.13 (b) 显示了形状参数与尺度参数倒数的回归线及其相关系数 \(\rho\), 结果表明两者之间存在显著的正相关关系.

随后采用 VB、EM 和贝叶斯方法对设备-B 数据进行拟合, 并使用所提 RGa 退化模型进行分析. 拟合结果汇总于表 \(\ref{tbl-tab:addlabe3}\). 从表中可以看出, VB 方法在计算效率上显著优于其他方法, 同时在精度方面表现良好. 图 3.14 (a) 展示了样本平均值及每种算法的估计平均退化路径 (\(\hat{\eta} t\)) . 可以看出, 各算法的估计平均退化路径均能很好地拟合样本平均值.

(a) 估计CDF
(b) 形状参数与尺度参数的关系
图 3.13: 模型验证相关结果.

接下来, 将设备-B数据分别用五种基准模型拟合, 其对数似然值和 AIC 值列于表 \(\ref{tbl-tab:addlabel4}\). 结果显示, 提出的 RGaLn 模型在所有候选模型中具有最小的AIC值, 可认为是最优模型. 图 3.14 (b) 中展示了经验分布函数以及退化阈值 \(\omega = 0.3\) 下基于 GaGa 和 RGaLn 模型估计的寿命分布. 结果表明, RGaLn 模型的拟合效果显著优于 GaGa 模型.

(a) 平均退化路径估计
(b) 寿命估计CDF
图 3.14: 设备-B数据下寿命分布相关结果.

3.3.5 附录

定理 3.3 证明

根据式\(\eqref{gammpdfx}\)和\(\eqref{lognorpdfx}\), \(Y(t)\) 和 \(\eta\) 的联合PDF可表示为 \[\begin{equation} \begin{aligned} p_{Y(t),\eta}(y_{t},\eta \mid \boldsymbol{\Theta}) = p_{Y(t)}(y_{t} \mid \eta,\alpha)p(\eta \mid \mu,\sigma^2), \end{aligned} \label{gamlnorpdfx} \end{equation}\] 其中\(\boldsymbol{\Theta} = (\alpha, \mu, \sigma^2)\). 接下来, 对式 \(\eqref{gamlnorpdfx}\)关于\(\eta\)求积分, 可得到 \(Y(t)\) 的边际PDF:
\[\begin{equation*} \begin{aligned} p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta}) =& \frac{\alpha^{\alpha t y_{t}^{\alpha t -1}}}{\Gamma(\alpha t) \sqrt{2\pi}\sigma} \int_{0}^{+\infty} \eta^{-\alpha t -1} e^{-\frac{\alpha y_{t}}{\eta}-\frac{(\ln(\eta) -\mu)^2}{2\sigma^2}} \diff\eta. \end{aligned} \end{equation*}\]
由于积分函数的复杂性和积分区域的约束, 这一积分的求解可能会导致数值不稳定的结果. 可通过重参数化或变量变换 \(z = (\log(\eta)-\mu)/\sqrt{2\sigma^2}\), 从而得到 \[\begin{equation}\label{ungamlnorpdfx1} p_{Y(t)}(y_{t} \mid \boldsymbol{\Theta}) = \frac{\alpha^{\alpha t} y_{t}^{\alpha t -1}}{\Gamma(\alpha t) \sqrt{\pi}} \int_{-\infty}^{+\infty} q(z,y_{t} \mid \boldsymbol{\Theta}) e^{-z^2} \diff z, \end{equation}\] 其中 \[ q(z,y_{t} \mid \boldsymbol{\Theta}) = \exp\left[-\alpha t(\mu + z\sqrt{2\sigma^2}) -\alpha y_{t} e^{-(\mu +z\sqrt{2\sigma^2})}\right].\]

观测 Fisher 信息矩阵

基于对数似然函数, 对\(\ell(\boldsymbol{\Theta})\)关于 \(\boldsymbol{\Theta}\) 中每个参数求一阶导数可得 \[\begin{align*} \frac{\partial \ell(\boldsymbol{\Theta})}{\partial \alpha} =& - NM h_{c}(\psi(\alpha h_{c})+\log(\alpha)+1) + \sum_{n=1}^{N}\sum_{m=1}^{M} h_{c}\log(x_{n,m}) + \sum_{n=1}^{N}A_{n}(\alpha),\nonumber\\ \frac{\partial \ell(\boldsymbol{\Theta})}{\partial \theta_{d}}=& \sum_{n=1}^{N}A_{n}(\theta_{d}), \end{align*}\] 其中 \[ A_{n}(s) = \frac{\int_{-\infty}^{\infty} \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial s}\exp(-z^2)\diff z}{\int_{-\infty}^{\infty} q(z,y_{n,M},\boldsymbol{\Theta})\exp(-z^2)\diff z}, \theta_{d}\in\{\mu,\sigma^2\}. \] \(q(z,y_{n,M},\boldsymbol{\Theta})\) 中每个元素的一阶导数为 \[\begin{align*} \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha} &= q(z,y_{n,M},\boldsymbol{\Theta})\cdot \log\left(q(z,y_{n,M},\boldsymbol{\Theta})\right)\cdot \alpha^{-1},\nonumber\\ \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu} &= q(z,y_{n,M},\boldsymbol{\Theta})\cdot\left( -\log\left(q(z,y_{n,M},\boldsymbol{\Theta})\right) -\alpha M h_{c}\left(\sqrt{2\sigma^2}z+\mu+1\right)\right) ,\nonumber\\ \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \sigma^2} &= \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu}\cdot \frac{z}{\sqrt{2\sigma^2}} .\nonumber\\ \end{align*}\] 接下来, 对\(\ell(\boldsymbol{\Theta})\)关于 \(\boldsymbol{\Theta}\)求二阶偏导数可得 \[\begin{align*} \frac{\partial \ell(\boldsymbol{\Theta})}{\partial \alpha^2} & = - NM h_{c}(h_{c}\psi'(\alpha h_{c})+1/\alpha)%\\ +\sum_{n=1}^{N} \left[B_{n}(\alpha,\alpha) + A_{n}(\alpha)^2\right], \nonumber\\ \frac{\partial \ell(\boldsymbol{\Theta})}{\partial \theta_{d}\partial \theta_{s}}&=\sum_{n=1}^{N} \left[B_{n}(\theta_{d},\theta_{s}) - A_{n}(\theta_{d})A_{n}(\theta_{s})\right], \end{align*}\] 其中 \[B_{n}(s_{1},s_{2})=\frac{\int_{-\infty}^{\infty} \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial s_{1} \partial s_{2}}\exp(-z^2)\diff z}{\int_{-\infty}^{\infty} q(z,y_{n,M},\boldsymbol{\Theta})\exp(-z^2)\diff z}, \theta_{d},\theta_{s}\in\boldsymbol{\Theta}, (\theta_{d},\theta_{s})\neq (\alpha, \alpha). \] \(q(z,y_{n,M},\boldsymbol{\Theta})\) 对 \(\boldsymbol{\Theta}\) 中每个元素的二阶导数为 \[\begin{align*} \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha^2} &= \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha}\cdot \log\left(q(z,y_{n,M},\boldsymbol{\Theta})\right)\cdot \alpha^{-1},\nonumber\\ \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha \partial \mu} &= \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu \partial \alpha} = \left(\log\left(q(z,y_{n,M},\boldsymbol{\Theta})\right)+ 1 \right)\times\frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu} \cdot\alpha^{-1},\nonumber \\ \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha \partial \sigma^2} & = \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \sigma^2\partial \alpha } =\frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \alpha \partial \mu}\cdot \frac{z}{\sqrt{2\sigma^2}},\\ \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu^2} &= \frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu}\left(-\log\left(q(z,y_{n,M},\boldsymbol{\Theta})\right) -\alpha M h_{c} \left(\sqrt{2\sigma^2}z+\mu+1\right)-1\right)\\ &-q(z,y_{n,M},\boldsymbol{\Theta})\cdot \alpha M h_{c},\nonumber\\ \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu \partial \sigma^2} & = \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \sigma^2 \partial \mu} %\\ = \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu^2} \cdot \frac{z}{\sqrt{2\sigma^2}},\nonumber\\ \frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial (\sigma^2)^2} &= \left(\frac{\partial^2 q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu \partial \sigma^2} -\frac{\partial q(z,y_{n,M},\boldsymbol{\Theta})}{\partial \mu} \cdot \frac{1}{2\sigma^2}\right)\times \frac{z}{\sqrt{2\sigma^2}}. \end{align*}\] 上述推导过程表明, 对式 \(\eqref{ci2}\) 进行数值计算并不困难. 首先, 式 \(\eqref{ci2}\) 中涉及的所有积分项都具有GH求积的结构. 其次, \(\eqref{ci2}\) 中的二阶偏导数项可以通过 \(\eqref{ci4}\) 转化为简单的算术运算, 这些也可以通过GH求积计算. 计算式 \(\eqref{ci2}\) 仅涉及一个低维的3x3矩阵的求逆运算, 这可以通过多种成熟的集成算法完成, 例如R软件中的solve()函数. 在实际操作中, 由于计算精度溢出可能导致矩阵求逆的误差, 这种情况可以通过将参数间的协方差设定为较小值来处理.

基准估计方法

EM算法通常用于对不可解析的似然函数求MLE, 特别是对于具有缺失数据的模型 (McLachlan 等, 2007). 对于所提伽马退化模型, 设隐变量 \(\boldsymbol{\eta}\) 为缺失数据, 然后可以使用EM算法来获得 \(\boldsymbol{\Theta}\) 的MLE, 其流程如下:

  1. E步: 给定第 \(s\) 次迭代后的估计值 \(\boldsymbol{\Theta}^{(s)}\), 计算后续M步所需的 \(\eta_{n}\) 的矩为 \[\begin{equation} \begin{aligned} E_{1n}^{(s)}&=E_{q_{1}(\eta_{n})}[\ln(\eta_{n})]|_{\boldsymbol{\Theta}=\boldsymbol{\Theta^{(s)}}},\nonumber\\ E_{2n}^{(s)}&=E_{q_{1}(\eta_{n})}[1/(\eta_{n})]|_{\boldsymbol{\Theta}=\boldsymbol{\Theta^{(s)}}},\\ E_{3n}^{(s)}&=E_{q_{1}(\eta_{n})}[(\ln(\eta_{n}))^2]|_{\boldsymbol{\Theta}=\boldsymbol{\Theta^{(s)}}},\nonumber \end{aligned}\label{em1} \end{equation}\] 其中 \[\begin{equation} \begin{aligned} q_{1}(\eta_{n}) \propto& \eta_{n}^{\mu^{(s)}/(\sigma^2)^{(s)}-\alpha^{(s)}Mh_{c}-1} e^{-\frac{\alpha^{(s)}y_{n,M}}{\eta_{n}}-\frac{(\ln(\eta_{n}))^2}{2(\sigma^2)^{(s)}}}. \end{aligned}\label{em2} \end{equation}\] 注意, 式 \(\eqref{em2}\) 也可以转换为GH求积, 并用于在E步中计算 \(E_{dn}^{(s)}, d=1,2,3\).

  2. M步: 使用 \(E_{dn}^{(s)}, d=1,2,3\) 更新参数 \(\boldsymbol{\Theta}\) 为 \[\begin{equation} \begin{aligned} \alpha^{(s+1)} &= \mathop{\textrm{arg min}}\limits_{\alpha \in \Omega_{\alpha}}\left\{\left\vert NM\left[h_{c}(\ln(\alpha)+1)-\psi(\alpha h_{c}) \right.\right.\right.\nonumber\\ &\ \ \ \ +\left.\left.\left.h_{c}\left(w_{5}-M\sum_{n=1}^{N}E_{1n}\right)-w_{6}\sum_{n=1}^{N}E_{2n}\right]\right\vert \right\},\nonumber\\ \mu^{(s+1)} &= \frac{1}{N}\sum_{n=1}^{N} E_{1n}^{(s)},\\ (\sigma^2)^{(s+1)}&=\frac{1}{N}\sum_{n=1}^{N}\left[E_{3n}-(\mu^{(s+1)})^2\right].\nonumber \end{aligned}\label{em3} \end{equation}\] EM算法的实现从参数的初始值 \(\boldsymbol{\Theta^{(0)}}\) 开始, 然后在 \(\eqref{em1}\)和\(\eqref{em3}\)之间迭代, 直到满足给定的收敛条件.

贝叶斯方法通过结合数据与先验信息, 推断目标参数的后验分布. 然而, 在贝叶斯分析中, 直接获得每个边际后验的解析形式通常较为困难. 为此, 可采用MCMC方法, 特别是Gibbs抽样算法, 从参数的满条件后验分布中生成样本, 逐步迭代直到马尔可夫链达到平稳时终止. 对于当前的RGa退化模型, 基于MCMC的框架如下:

  1. 在当前后验样本 \(\boldsymbol{\Theta}^{(s)}\) 下, 从 \(\eta_{n}\) 的满条件后验中生成 \(\eta_{n}^{(s+1)}\), 即 \[\begin{equation}\label{gib1} p_{2}(\eta_{n} \mid \boldsymbol{\Theta^{(s)}}) = \prod_{m=1}^{M} p_{x \mid \eta}(x_{n,m} \mid \alpha^{(s)}). \end{equation}\]

  2. 在当前后验样本 \(\boldsymbol{\eta}^{(s+1)}\) 下, 从 \(\alpha\) 的满条件后验中生成 \(\boldsymbol{\Theta}^{(s+1)}\), 即 \[\begin{equation}\label{gib2} p_{2}(\alpha \mid \boldsymbol{\eta}^{(s+1)}) = \prod_{n=1}^{N}\prod_{m=1}^{M} p_{x \mid \eta}(x_{n,m} \mid \eta_{n}^{(s+1)},\alpha) p(\alpha), \end{equation}\] 以及 从\((\mu, \sigma^2)\) 的满条件后验中采样: \[\begin{equation}\label{gib3} p_{2}(\mu,\sigma^2 \mid \boldsymbol{\eta}^{(s+1)}) = \prod_{n=1}^{N} p_{\eta}(\eta_{n}^{(s+1)} \mid \mu,\sigma^2) p(\mu,\sigma^2), \end{equation}\] 其中 \(p(\alpha)\) 和 \(p(\mu,\sigma^2)\) 是 \(\alpha\) 和 \((\mu,\sigma^2)\)的先验, 与VB算法中的设置相同. 贝叶斯方法的实现从设置初始值 \(\boldsymbol{\Theta}^{(0)}\) 和 \(\boldsymbol{\eta}^{(0)}\) 开始, 并基于式 \(\eqref{gib1}\)-\(\eqref{gib3}\) 迭代生成后验样本, 直到所有链达到稳定且充分混合.