可靠性数据分析

寿命数据分析

简介

寿命数据分析（Life Data Analysis，LDA）也称为威布尔分析或可靠性分析，是用于评估产品或系统失效时间数据的统计方法。它帮助工程师和决策者了解产品在失效前可能运行多长时间，从而对可靠性、可维修性和性能作出有依据的预测。通过将寿命数据拟合至威布尔、指数或对数正态等概率分布，LDA 可以得到平均失效时间（MTTF）、失效率和随时间变化的可靠度等重要指标。这些信息可用于改进产品设计、保修规划和维修安排。

LDA 广泛应用于航空航天、制造、电子和能源等重视运行时间与安全性的行业。它通过量化失效行为和识别主要失效模式，支持以数据为依据的决策。LDA 不仅能够处理完整失效数据，也能处理删失（不完整）数据，例如试验结束时仍在运行的样品，或因其他原因退出观察的样品。因此，它为提高系统可靠性和优化寿命周期成本提供了灵活而有力的分析框架。

统计学基础

随机变量

随机变量是代表随机试验结果的数值。它为试验的各种可能结果赋予数值，使我们能够利用数学方法分析和描述不确定性。随机变量是概率论和统计学的基础，用于量化结果并计算其发生概率。

随机变量主要分为离散型和连续型两类。离散随机变量只能取可数的数值，例如试验中出现 0、1、2 或 3 次失效；它通常来自计数或分类过程。连续随机变量则可以在某个范围内取任意数值，例如某个部件的准确失效时间。连续变量通常通过测量得到，并以概率密度函数（PDF）描述；离散变量则使用概率质量函数（PMF）。

寿命数据分析只处理连续随机变量。

概率密度函数（PDF）

假设以下数据是 1,000 个轴承的失效时间（天）。

以 100 天为组距绘制直方图：

假设我们要估计运行至 300 天时总体中已经失效的百分比。300 天之前有三个柱，其高度分别为 6、31 和 74，合计得到 111 次失效。

样本总数为 1,000 个，因此轴承在运行 300 天内失效的概率为 0.111。请注意，以这个直方图只能按照 100 天的组距查询 100、200、300 天等时点的失效概率。

现在在直方图上叠加一条函数曲线：

该函数应具有以下特性：

1. 函数轮廓反映并平滑各柱的数值。
2. 函数曲线下的总面积等于 1。

有了这个函数，我们可以查询任意目标时间的失效概率。例如，将该函数从 0 积分至 150 天，即可求得 150 天内的失效概率。

这个函数称为概率密度函数（Probability Density Function，PDF），记为 $f(t)$。

概率密度函数是描述连续随机变量在某一数值附近出现可能性的数学函数。变量落在指定区间内的概率，由 PDF 曲线在该区间下方的面积确定。整条 PDF 曲线下的总面积恒等于 1，表示变量必然在其定义域内取某个数值。

累积分布函数（CDF）

PDF 可以直观显示随机变量的分布，但要直接读取函数曲线截至某个时点的面积并不方便。累积分布函数（CDF）则允许我们直接读取时间 $t$ 时的失效概率。

大多数经典统计教材以 $F(t)$ 表示 CDF；可靠性教材则常用 $Q(t)$。CDF 也称为不可靠度函数。

如下图所示，可直接从 CDF 读取：$Q(t=300\text{ 天})=0.11$。

若从 PDF 读取该数值，必须求出 $f(t)$ 从 $t=0$ 至 $t=300$ 天的曲线下面积，这并非简单的工作。CDF 是 PDF 的积分：

$Q(t_{0}) = P(t \leq t_{0}) = \int_{- \infty}^{t_{0}}{f(t)dt}$

反之，PDF 是 CDF 的导数：

$f(t) = \frac{d(Q(t))}{dt}$

可靠度与不可靠度

不可靠度图（CDF）允许用户直接读取时间 $t$ 时的失效概率。时间 $t$ 时的可靠度为：

$R(t) = 1 - Q(t)$

利用上述关系可以得到可靠度曲线。

失效率

数学上，时间 $t$ 时的失效率 $\lambda(t)$ 定义为：在总体已经存活至时间 $t$ 的条件下，单位时间内发生失效的概率。

根据定义：

其中，$N(t)$ 是时间 $t$ 时的存活数量，$N_{0}$ 是初始总体数量。

条件可靠度

条件可靠度 $R(t/T)$ 是指：某单元已经累积运行至年龄 $T$ 的条件下，在随后任务时间 $t$ 内继续存活的概率。

要使某单元从 A 点存活至 C 点，它必须先从 A 点存活至 B 点，再从 B 点存活至 C 点。因此：

均值与中位数

平均失效时间（MTTF）$\overline{T}$ 定义为：

$\overline{T} = \int_{0}^{\infty}t \bullet f(t)dt$

中位寿命 $\breve{T}$ 是总体中 50% 已经失效时所对应的时间：

$Q\left( \breve{T} \right) = 0.5$

可靠性数据分类

寿命数据分析使用以下主要数据类型：

1. 完整（失效）数据：每个样品均观察至失效，准确失效时间已知。这是最直接的数据类型。
2. 右删失数据（暂停数据）：样品在观察期内没有失效，只知道它至少存活至某个时间。实际试验和现场数据中经常出现这种情况。
3. 区间删失数据：准确失效时间未知，但已知失效发生在某个时间区间内，例如定期检查时发现样品在两次检查之间失效。
4. 左删失数据：已知失效发生在某个时间之前，但没有观察到准确失效时间。

这些数据类型可以共同用于估计寿命分布和可靠性指标。

分布模型

分布模型用于从统计角度描述部件、系统或产品的失效时间。它帮助工程师了解潜在失效行为，例如失效是随机发生、随年龄增加，还是遵循某种退化模式。将合适的分布拟合至观测失效数据后，可以估计失效率、随时间变化的可靠度和平均失效时间等关键指标，从而支持设计、维修、保修预测和风险管理方面的决策。

指数分布（λ）

指数分布通常用于描述随机且无记忆事件之间的时间，特别适用于失效率随时间保持恒定的情况。它结构简单，但在可靠性和生存分析中十分重要，并常被用作基准模型。

概率密度函数（PDF）

指数分布的 PDF 为：

$f(t) = {\lambda e}^{- \lambda t}$

当 $\lambda$ 减小时，指数分布向右伸展。为保持 PDF 曲线下总面积恒等于 1，曲线高度会相应降低。

累积分布函数（CDF）

当 $\lambda$ 减小时，指数分布向右伸展。

失效率函数

平均失效时间（MTTF）

中位寿命

中位寿命 $\breve{T}$ 是总体中 50% 已经失效时所对应的时间：

威布尔分布（β，η）

威布尔分布常用于描述失效时间数据。它具有很高的灵活性，可通过形状参数 $\beta$ 表示多种不同的失效行为。

概率密度函数（PDF）

威布尔分布的 PDF 为：

$f(t;\beta,\eta) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

• 当 $0 < \beta < 1$：$f(t)$ 在起点为无穷大，随后迅速下降。
• 当 $\beta = 1$：呈指数分布形状。
• 当 $\beta > 1$：$f(t)$ 在起点为 0，随后上升至峰值再下降。

当 $\beta$ 保持不变而 $\eta$ 增大时，威布尔分布向右伸展。曲线整体形状保持不变，但峰值降低，以保证 PDF 曲线下总面积仍为 1。

累积分布函数（CDF）

CDF 也称为不可靠度，通常记为 $Q(t)$ 或 $F(t)$：

$Q(t;\beta,\eta) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

• 当 $0 < \beta < 1$：初期迅速上升，之后逐渐趋缓。
• 当 $\beta = 1$：呈指数分布形状。
• 当 $\beta > 1$：呈 S 形，初期增长缓慢，随后快速上升。

当 $\beta$ 保持不变而 $\eta$ 增大时，CDF 曲线向右伸展，整体形状不变。

失效率函数

失效率函数 $\lambda(t)$ 表示单位时间内发生失效的比例：

$\lambda(t) = \frac{f(t)\text{~}}{R(t)}$

代入 $f(t) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

和 $R(t) = e^{- \left( \frac{t}{\eta} \right)^{\beta}}$：

• 当 $0 < \beta < 1$：初始失效率较高，之后随时间下降，适合描述早期失效。
• 当 $\beta = 1$：失效率恒定。
• 当 $\beta > 1$：失效率随时间上升，适合描述磨损失效。

减小 $\eta$ 会使失效提前，因此相同时间 $t$ 下的失效率更高。

平均失效时间（MTTF）

威布尔分布的均值为：

$\overline{T} = \eta \bullet \Gamma\left( \frac{1}{\beta} + 1 \right)$

其中：

$\Gamma(n) = \int_{0}^{\infty}{e^{- x}x^{n - 1}dx}$

中位寿命

正态分布（μ，σ）

正态分布是一种对称的钟形概率分布，在统计学中应用广泛。由于它允许出现负值，而失效时间不可能为负，因此在寿命数据分析中不如威布尔或对数正态分布常用。但当变异较小且失效集中在某个中心值附近时，正态分布仍可能适用。

概率密度函数（PDF）

正态分布的 PDF 为：

$f(t;\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{t - \mu}{\sigma} \right)^{2}}$

$\mu$ 决定分布的中心（均值）或位置；$\sigma$ 控制数据围绕 $\mu$ 的离散程度。

累积分布函数（CDF）

$Q(t) = \frac{1}{\sigma\sqrt{2\pi}}\int_{- \infty}^{t}e^{- \frac{1}{2}\left( \frac{x - \mu}{\sigma} \right)^{2}}\mathbb{d}x$

正态 CDF 也可以用标准正态 CDF 表示：

$Q(t) = \Phi\left( \frac{t - \mu}{\sigma} \right)$

其中，$\Phi(\bullet)$ 为标准正态 CDF。

失效率函数

正态分布的失效率为：

$\lambda(t) = \frac{f(t)}{R(t)}$

其中：

$f(t) = \frac{1}{\sigma\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{t - \mu}{\sigma} \right)^{2}}$

$R(t) = 1 - \frac{1}{\sigma\sqrt{2\pi}}\int_{- \infty}^{t}e^{- \frac{1}{2}\left( \frac{x - \mu}{\sigma} \right)^{2}}\mathbb{d}x$

平均失效时间（MTTF）与中位寿命

由于正态分布对称，中位数等于均值，也就是分布参数 $\mu$：

$\overline{T} = \breve{T} = \mu$

对数正态分布（μ'，σ'）

对数正态分布用于描述其对数值服从正态分布的失效时间。它特别适用于失效机理由疲劳、腐蚀或化学分解等损伤或退化逐步累积所驱动的产品或系统。对数正态模型为右偏分布：多数单元在某一时间附近失效，但少数单元可以运行得更久。其特征是初始失效率较低，随后随时间上升，适合描述失效时间变异较大的磨损行为。

概率密度函数（PDF）

对数正态分布的 PDF 为：

$f(t) = \frac{1}{{t \bullet \sigma}'\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{ln(t) - \mu'}{\sigma'} \right)^{2}\text{~}}$

其中：

• $\mu'$：$t$ 的自然对数均值（对数位置参数）。
• $\sigma'$：$t$ 的自然对数标准差（对数尺度参数）。

$\sigma'$ 控制分布的离散程度和偏度。较大的 $\sigma'$ 使曲线更平、更宽，表示出现很长寿命的可能性更高；较小的 $\sigma'$（小于 1）使曲线更窄、更高，表示失效更集中在中位寿命附近。

$\mu'$ 控制 PDF 的水平位置（时间尺度），但不改变曲线形状。

累积分布函数（CDF）

$Q(t) = \frac{1}{t\sigma'\sqrt{2\pi}}\int_{0}^{t}{e^{- \frac{1}{2}\left( \frac{ln(s) - \mu'}{\sigma'} \right)^{2}\text{~}}ds}$

也可以表示为：

$Q(t) = \Phi\left( \frac{\ln(t) - \mu'}{\sigma'} \right)$

其中，$\Phi(\bullet)$ 为标准正态 CDF。

较大的 $\sigma'$ 在水平方向拉伸 CDF，表示失效分散在更宽的时间范围；较小的 $\sigma'$ 则压缩 CDF 并使其更陡，表示失效集中在较窄的时间区间。

$\mu'$ 控制 CDF 的水平位置（时间尺度），但不改变其形状。

失效率函数

对数正态分布的失效率为：

$\lambda(t) = \frac{f(t)}{R(t)}$

其中：

$f(t) = \frac{1}{{t \bullet \sigma}'\sqrt{2\pi}}e^{- \frac{1}{2}\left( \frac{ln(t) - \mu'}{\sigma'} \right)^{2}\text{~}}$

$R(t) = 1 - \frac{1}{t\sigma^{'\sqrt{2\pi}}}\int_{0}^{t}{e^{- \frac{1}{2}\left( \frac{\ln(s) - \mu'}{\sigma'} \right)^{2}\text{~}}ds}$

平均失效时间（MTTF）

该积分在统计学中的已知结果为：

$\overline{T} = e^{\mu + \frac{1}{2}{\sigma'}^{2}}$

中位寿命

$Q\left( \breve{T} \right) = \Phi\left( \frac{\ln\left( \breve{T} \right) - \mu'}{\sigma'} \right) = 0.5$

由于 $\Phi(0)=0.5$：

参数估计

参数估计的目的是确定最能描述一组观测数据的分布参数值。以下介绍两种常用方法：概率图法和极大似然估计（MLE）。它们常用于利用样本数据拟合威布尔、对数正态、正态和指数等参数分布。

概率图法

概率图法是一种以图形方式估计理论概率分布参数（例如均值、标准差和形状参数）的方法，使所选分布尽可能拟合给定数据。

以下使用数据集 {30, 40, 20, 50} 说明。首先对失效数据排序，得到 {20, 30, 40, 50}。

对于每个失效时间 $t_i$，使用累积二项式方程计算累积概率 $q_i$：

$0.5 = \sum_{k = j}^{N}{\binom{N}{k}q^{k}(1 - q)^{N - k}}$

每个 $t_j$ 对应一个累积概率 $q_j$。下一步是建立概率纸。对 X 轴和 Y 轴进行变换，使目标分布的 CDF 在概率纸上呈直线。

以下各节说明如何建立概率纸，并求解指数、威布尔、正态和对数正态分布的参数。

指数分布

指数分布的 CDF 为：

两边取自然对数：

$\ln\left( R(t) \right) = - \lambda t$

因此，$\ln(R(t))$ 对 $t$ 作图将形成一条斜率为 $-\lambda$ 的直线。

将坐标 $(t,\ln(R(t)))$ 绘制在指数概率纸上。由于 $R(t=0)=1$，拟合直线必须通过 $(t=0,\ln(R(0))=0)$。

为求斜率，可在拟合线上选取一个已知 Y 值的点。当 $Y=-2.303$ 时，对应时间 $t=75$。

指数分布参数为：

$\lambda = -\text{斜率} = - \frac{0 - \left( - \text{2.303~} \right)}{75} = 0.0307$

单位为每单位时间。

威布尔分布

威布尔分布的 CDF 可以写为：

$\Rightarrow \ ln\left( \ln\left( \frac{1}{1 - Q(t)} \right) \right) = \beta \bullet ln(t) - \beta \bullet ln(\eta)$

将方程线性化为 $Y=\beta X+C$，其中：

$Y = ln\left( \ln\left( \frac{1}{1 - Q(t)} \right) \right)$

$X = ln(t)$

$\text{斜率} = \beta$

因此，Y 对 X 作图将形成一条斜率为 $\beta$ 的直线。

将坐标 $(X,Y)$ 绘制在威布尔概率纸上。

由威布尔 CDF：

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

令 $t=\eta$，则：

$Q(\eta) = 1 - e^{- (1)^{\beta}} = 63.1\%$

因此，$Q=63.1\%$ 对应的时间就是 $\eta$（本例为 40）。

为求 $\beta$，可在拟合线上选择两个已知 Y 值的点并计算斜率。例如：

• 当 $Z=0.834$ 时，$Q(t)=0.9$。
• 当 $Z=-0.367$ 时，$Q(t)=0.5$。

$\beta = \frac{0.834 - ( - 0.367)}{\ln(55.4) - ln(34.6)} = 2.56$

正态分布

对于正态分布 $Nor(\mu,\sigma)$：

$Q(t) = \Phi\left( \frac{t - \mu}{\sigma} \right)$

其中，$\Phi^{-1}$ 为标准正态 CDF 的反函数。若以 $\Phi^{-1}(Q(t))$ 为 Y 轴，以未经变换的 $t$ 为 X 轴，则正态分布的 CDF 在变换后的概率纸上成为直线，其 Y 轴截距为 $-\mu/\sigma$，斜率为 $1/\sigma$。

将累积概率 $Q(t)$ 转换为标准正态分位数（Z 值）：

$Z = \Phi^{- 1}\left( Q(t) \right)$

可以使用标准正态表或 Excel 函数 NORM.S.INV(Q) 得到 Z 值。

将坐标 $(t,Z)$ 绘制在正态概率纸上。

数学上可由 Y 轴截距 $-\mu/\sigma$ 和斜率 $1/\sigma$ 求得均值和标准差。对于正态分布，$Q(t=\mu)=50\%$，从图中可得 $\mu=35$。

为求 $\sigma$，在拟合线上选取两个已知 Z 值的点（例如 $Z=1$ 和 $Z=-1$），再取斜率的倒数。

• 当 $Z=1$ 时，$Q(t)=0.841$。
• 当 $Z=-1$ 时，$Q(t)=0.159$。

可以使用标准正态表或 Excel 函数 NORM.S.DIST(Z) 得到累积概率。

标准差为：

对数正态分布

若变量 $t$ 服从对数正态分布，则 $\ln(t)$ 服从正态分布。因此，只需对失效时间取对数，再与标准正态分位数作图。

对数正态分布 $Log(\mu',\sigma')$ 的 CDF 为：

$Q(t) = \Phi\left( \frac{\ln(t) - \mu'}{\sigma'} \right)$

$\Rightarrow \ \Phi^{- 1}\left\lbrack Q(t) \right\rbrack = - \frac{\mu'}{\sigma'} + \frac{1}{\sigma'}\ln(t)$

其中：

• $\mu'$ 是对数变换后失效时间的均值。
• $\sigma'$ 是对数变换后失效时间的标准差。

若以 $\Phi^{-1}(Q(t))$ 为 Y 轴，以 $\ln(t)$ 为 X 轴，对数正态分布的 CDF 将形成一条斜率为 $1/\sigma'$ 的直线。

将 Y 轴的累积概率转换为标准正态分位数：

$Z=\Phi^{-1}(Q(t))$

同时将 X 轴转换为 $\ln(t)$。可使用标准正态表或 Excel 函数 NORM.S.INV(Q) 求取 Z 值。

将坐标 $(\ln(t),Z)$ 绘制在对数正态概率纸上。

对数均值 $\mu'$ 出现在 $Z=0$，即 $Q(t)=50\%$ 的位置，此时 $\ln(t)=\mu'$。本例中 $Q(t)=50\%$ 对应的 $t=33.1$，因此：

$\mu'=\ln(33.1)=3.5$

为求 $\sigma'$，在拟合线上选择两个已知 Z 值的点（例如 $Z=1$ 和 $Z=-1$），并计算斜率的倒数。本例中相应的 $t$ 值为 18.3 和 52.5。

对数标准差为：

秩回归

秩回归通过对观测数据拟合直线来估计分布参数，其目标是使观测值与拟合直线预测值之间的平方误差之和最小。

Y 方向秩回归

RRY（Rank Regression on Y）拟合线使各数据点与回归线之间 Y 方向的误差最小。当需要根据给定概率预测时间（例如在不可靠度为 90% 时预测失效时间）时，通常使用 RRY。

给定成对数据点 $(x_1,y_1),(x_2,y_2),\ldots,(x_n,y_n)$，并假定 X 值准确已知：

$y_{i} = a' + b'x_{i} + e_{i}$

其中：

$b' = \frac{\sum_{i = 1}^{n}{\left( y_{i} - \overline{y} \right)\left( x_{i} - \overline{x} \right)}}{\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}}$

$a' = \overline{y} - b^{'\overline{x}}$

$\overline{x} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}$

$\overline{y} = \frac{1}{n}\sum_{i = 1}^{n}y_{i}$

$e_i$ 为误差项。

X 方向秩回归

RRX（Rank Regression on X）拟合线使各数据点与回归线之间 X 方向的误差最小。当需要根据给定时间预测概率时，通常使用 RRX。

给定成对数据点，并假定 Y 值准确已知，为使 X 方向误差最小，对 X 关于 Y 进行回归：

$x_{i} = a' + b'y_{i} + e_{i}$

其中：

$b' = \frac{\sum_{i = 1}^{n}{\left( y_{i} - \overline{y} \right)\left( x_{i} - \overline{x} \right)}}{\sum_{i = 1}^{n}\left( y_{i} - \overline{y} \right)^{2}}$

$a' = \overline{x} - b'\overline{y}$

$\overline{x} = \frac{1}{n}\sum_{i = 1}^{n}x_{i}$

$\overline{y} = \frac{1}{n}\sum_{i = 1}^{n}y_{i}$

$e_i$ 为误差项。

相关系数

秩回归中的一个重要概念是相关系数，用来衡量所选分布对数据的拟合程度。相关系数 $\rho$ 定义为：

其中：

• $\overline{x}$ 和 $\overline{y}$ 为 X 和 Y 的样本均值。
• $Cov(X,Y)$ 为协方差。
• $\sigma_X$ 和 $\sigma_Y$ 为 X 和 Y 的标准差。

解释如下：

• 数值范围为 -1 至 +1。
• -1 表示负斜率的完全拟合。
• +1 表示正斜率的完全拟合。

考虑数据集 {20, 30, 40, 50}。下图显示使用 RRX 将该数据拟合至指数、威布尔、正态和对数正态模型时得到的 $\rho$ 值。

在概率图法（秩回归）中，$\rho$ 可作为所选分布拟合数据优劣的指标。

极大似然估计

极大似然估计（Maximum Likelihood Estimation，MLE）是一种估计概率分布参数的统计方法。它寻找一组参数，使在假定分布下观察到当前数据的可能性最大。其核心是建立一个似然函数，将观测数据的概率或概率密度表示为未知参数的函数，再求使该函数达到最大值的参数。

令：

• $f(t;\theta)$：分布的概率密度函数，例如威布尔或对数正态分布。
• $F(t;\theta)$：累积分布函数。
• $R(t;\theta)=1-F(t;\theta)$：可靠度函数（存活概率）。
• $\theta$：未知分布参数向量，例如形状参数和尺度参数。
• $n_1$：在时间 $t_i$ 发生的准确失效。
• $n_2$：在时间 $t_i$ 的右删失观测。
• $n_3$：在 $t_{L_i}$ 与 $t_{U_i}$ 之间的区间删失观测。

完整似然函数为：

$L(\theta) = \prod_{i = 1}^{n_{1}}{f\left( t_{i};\theta \right) \times \text{~}}\prod_{i = 1}^{n_{2}}{R\left( t_{i};\theta \right) \times \prod_{i = 1}^{n_{3}}{\left\lbrack F\left( t_{U_{i}};\theta \right) - F\left( t_{L_{i}};\theta \right) \right\rbrack\text{~}}}$

直接求乘积的最大值较为复杂，因此通常使用对数似然函数：

$\Lambda = ln\left( L(\theta) \right) = \sum_{i = 1}^{n_{1}}{\ln\left( f\left( t_{i};\theta \right) \right) \times \sum_{i = 1}^{n_{2}}{\ln\left( R\left( t_{i};\theta \right) \right) \times \sum_{i = 1}^{n_{3}}{\ln\left\lbrack F\left( t_{U_{i}};\theta \right) - F\left( t_{L_{i}};\theta \right) \right\rbrack\text{~}}}}$

MLE 寻找使对数似然 $\Lambda$ 最大的参数 $\theta$。

以下示例将一个数据集采用 MLE 拟合至威布尔分布。

PDF：

$f(t;\beta,\eta) = \frac{\beta}{\eta}\left( \frac{t}{\eta} \right)^{\beta - 1}e^{{- \left( \frac{t}{\eta} \right)}^{\beta}}$

CDF：

$Q(t;\beta,\eta) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

可靠度：

$R(t;\beta,\eta) = e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

对数似然函数：

$\small\Lambda = ln\left( f(20;\beta,\eta) \right) + ln\left( f(30;\beta,\eta) \right) + 2 \times ln\left( R(40;\beta,\eta) \right) + ln\left\lbrack F(30;\beta,\eta) - F(50;\beta,\eta) \right\rbrack$

代入威布尔分布函数：

$\small\Lambda = \left\lbrack \ln\left( \frac{\beta}{\eta} \right) + (\beta - 1)\ln(20) - \left( \frac{20}{\eta} \right)^{\beta} \right\rbrack + \left\lbrack \ln\left( \frac{\beta}{\eta} \right) + (\beta - 1)\ln(30) - \left( \frac{30}{\eta} \right)^{\beta} \right\rbrack$

$\small- 2\left( \frac{40}{\eta} \right)^{\beta} + ln\left( e^{- \left( \frac{30}{\eta} \right)^{\beta}} - e^{- \left( \frac{50}{\eta} \right)^{\beta}} \right)$

上述对数似然是 $\beta$ 和 $\eta$ 的函数。MLE 求取使 $\Lambda$ 最大的 $\beta$ 和 $\eta$。由于该函数通常为非线性且较复杂，需要使用数值优化方法。Weibull Toolbox 使用 Nelder–Mead 优化方法求解。

本例得到 $\beta=3.37$、$\eta=41.9$。将其代回对数似然方程：

$\small\Lambda = \left\lbrack \ln\left( \frac{\text{3.37}}{\text{41.9}} \right) + \left( \text{3.37} - 1 \right)\ln(20) - \left( \frac{20}{\text{41.9}} \right)^{\text{3.37}} \right\rbrack + \left\lbrack \ln\left( \frac{\text{3.37}}{\text{41.9}} \right) + \left( \text{3.37} - 1 \right)\ln(30) - \left( \frac{30}{\text{41.9}} \right)^{\text{3.37}} \right\rbrack$

在 Weibull Toolbox 中，对数似然值以 LK Value 显示在分析结果中。

若将同一数据集分别拟合至正态、对数正态、威布尔和指数分布，可比较它们的 LK Value：

正态分布得到最高（即最接近零）的 LK Value，因此在这些候选分布中拟合最佳。指数分布的 LK Value 明显较低，表示其对数据的拟合较差。

置信限

方法选择。 没有一种置信限方法在所有情况下都必然优于其他方法。Fisher 矩阵法、似然比法和模拟法采用不同的近似方法，特别是在数据量有限或存在大量删失数据时，可能得到不同的结果。所选方法应与数据、拟合模型以及需要解决的工程问题相适应。

Fisher 矩阵置信限

将威布尔模型拟合至失效数据后，可以利用 Fisher 矩阵置信限（FMCB）量化参数估计值 $\widehat{\beta}$ 和 $\widehat{\eta}$ 的不确定性，并由此求得 $\beta$ 和 $\eta$ 的方差。

参数的 Fisher 矩阵置信限

以下使用威布尔分布说明 Fisher 矩阵置信限；相同概念也适用于其他寿命分布模型。

二参数威布尔分布的似然函数为：

$L(\beta,\eta) = \prod_{i = 1}^{n}\left\lbrack \left( \frac{\beta}{\eta}\left( \frac{t_{i}}{\eta} \right)^{\beta - 1} \right)^{\delta_{i}} \cdot exp\left( - \left( \frac{t_{i}}{\eta} \right)^{\beta} \right) \right\rbrack$

其中，$t_i$ 是第 $i$ 个样品的观测时间；失效数据的 $\delta_i=1$，删失数据的 $\delta_i=0$。

对数似然函数为：

$\ln\left( L(\beta,\eta) \right) = \sum_{i = 1}^{n}\delta_{i}\left\lbrack \ln(\beta) - \beta ln(\eta) + (\beta - 1)\ln\left( t_{i} \right) \right\rbrack - \sum_{i = 1}^{n}\left( \frac{t_{i}}{\eta} \right)^{\beta}$

威布尔分布的 Fisher 信息矩阵

Fisher 信息矩阵 $I(\beta,\eta)$ 描述对数似然函数在 MLE 点附近的曲率，定义为二阶导数负期望：

其逆矩阵给出估计值的协方差矩阵：

假设 $\ln(\beta)$ 和 $\ln(\eta)$ 近似服从正态分布，以确保 $\beta$ 和 $\eta$ 的置信限始终为正。

形状参数 $\beta$ 的置信限

$\ln\left( \beta_{U,L} \right) = ln\left( \widehat{\beta} \right) \pm z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \ln\left( \widehat{\beta} \right) \right)}$

利用 Delta 方法：

$Var(ln(\widehat{\beta}) \approx \frac{Var(\widehat{\beta})}{\beta^{2}}$

因此：

$\ln\left( \beta_{U,L} \right) = \ln\left( \widehat{\beta} \right) \pm z_{\frac{\alpha}{2}} \cdot \frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}$

对于 90% 置信水平，$\alpha=0.1$，$z_{\alpha/2}=1.64$。

置信水平 $\alpha$ 下 $\beta$ 的双侧置信限为：

$\beta_{U,L} = \widehat{\beta} \bullet e^{\pm z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}}$

尺度参数 $\eta$ 的置信限

$\ln\left( \eta_{U,L} \right) = ln\left( \widehat{\eta} \right) \pm z_{\frac{\alpha}{2}} \cdot \frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}$

$\eta_{U,L} = \widehat{\eta} \bullet e^{\pm z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\eta} \right)}}{\widehat{\eta}}}$

示例

将数据集 {20, 30, 40, 50} 使用 RRX 拟合至二参数威布尔分布，求参数的 90% 置信限。

RRX 解为：$\beta=2.55$，$\eta=39.9$。

Fisher 信息矩阵：

协方差矩阵：

因此：

$Var\left( \widehat{\beta} \right)=1.215$

$Var\left( \widehat{\eta} \right)=76.42$

$\beta_{U} = \widehat{\beta} \bullet e^{z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}}$

$\beta_{L} = \widehat{\beta} \bullet e^{- z_{\frac{\alpha}{2}}\frac{\sqrt{Var\left( \widehat{\beta} \right)}}{\widehat{\beta}}}$

$\eta_{U} = 39.9 \times e^{1.64\left( \frac{\sqrt{76.42}}{39.9} \right)} = 57.23$

$\eta_{L} = 39.9 \times e^{-1.64\left( \frac{\sqrt{76.42}}{39.9} \right)} = 27.84$

Fisher 矩阵置信限的概念可概括为：

1. 通过统计分析（本例为 MLE 或 RRX）得到估计值，作为均值。
2. 由 Fisher 信息矩阵求得标准差。
3. 使用由标准正态分位数定义置信水平的正态分布，对估计值建立置信限。

时间的 Fisher 矩阵置信限

将数据集 {20, 30, 40, 50} 拟合至二参数威布尔分布，求总体中 10% 失效（不可靠度 $Q=10\%$）所对应时间的 90% 置信限。

RRX 解为：$\beta=2.55$，$\eta=39.9$。

威布尔 CDF：

$Q(t) = 1 - e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

当 $Q=10\%$ 时，$t=16.5$。以下在均值为 16.5 的时间上建立正态置信限：

由于 X 轴采用自然对数尺度 $\ln(t)$，实际均值是 $\ln(16.5)$。还需要方差 $Var(\ln(\widehat{t}))$ 才能完整定义正态置信限。

令 $u=\ln(t)$，则：

$u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

$u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

已知 $Var(\widehat{\beta})$ 和 $Var(\widehat{\eta})$，但不知道 $Var(\widehat{u})$，因此使用 Delta 方法求取。

由威布尔 CDF 线性化：

$lnln\left( \frac{1}{1 - Q(t)} \right) = \beta\left\lbrack ln(t) - ln(\eta) \right\rbrack$

$u=ln(t)=\frac{1}{\beta}lnln\left( \frac{1}{1 - Q(t)} \right)+ln(\eta)$

利用 Delta 方法：

$Var(u) = \left( \frac{\partial u}{\partial\beta} \right)^{2}Var(\beta) + \left( \frac{\partial u}{\partial\eta} \right)^{2}Var(\eta) + 2\left( \frac{\partial u}{\partial\beta} \right)\left( \frac{\partial u}{\partial\eta} \right)Cov(\beta,\eta)$

$= \frac{1}{\beta^{4}}\left( lnln\left( \frac{1}{(1 - Q)} \right) \right)^{2}Var(\beta) + \left( \frac{1}{\eta} \right)^{2}Var(\eta) - \frac{2}{{\eta\beta}^{2}}\left( lnln\left( \frac{1}{(1 - Q)} \right) \right)Cov(\beta,\eta)$

代入：

$Var(\widehat{\beta})=1.215,\quad Var(\widehat{\eta})=76.42,\quad Q=0.1,\quad \widehat{\beta}=2.55,\quad \widehat{\eta}=39.9$

得到：

$Var(\widehat{u})=0.1854$

$\widehat{u}=\ln(16.5)=2.804$

因此：

$u_U=3.513,\qquad u_L=2.096$

由于 $u=\ln(t)$，所以 $t=\exp(u)$：

$t_U=\exp(u_U)=33.5$

$t_L=\exp(u_L)=8.13$

在 $Q=10\%$ 处得到时间 $t$ 的 90% 置信区间。对不同的 Q 值重复计算：

将概率图上的各点连接起来，即得到置信限曲线：

可靠度的 Fisher 矩阵置信限

将数据集 {20, 30, 40, 50} 拟合至二参数威布尔分布，求时间 $t=20$ 时不可靠度的 90% 置信限。

RRX 解为：$\beta=2.55$，$\eta=39.9$。由威布尔 CDF，当 $t=20$ 时，$Q=15.8\%$。

以下在均值为 15.8% 的不可靠度上建立正态置信限：

由于 Y 轴采用双对数倒数尺度 $lnln(1/(1-Q(t)))$，实际均值为 $lnln(1/(1-0.158))$。

令：

$u = lnln\left( \frac{1}{1 - Q(t)} \right)$

则：

$u_{U} = \widehat{u} + z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

$u_{L} = \widehat{u} - z_{\frac{\alpha}{2}} \cdot \sqrt{Var\left( \widehat{u} \right)}$

由威布尔 CDF 线性化：

$u=\beta\left\lbrack ln(t)-ln(\eta)\right\rbrack$

应用 Delta 方法：

$Var(u) = \left( \ln(t) - ln(\eta) \right)^{2}Var(\beta) + \left( \frac{\beta}{\eta} \right)^{2}Var(\eta) - \frac{2\beta}{\eta}\left( \ln(t) - ln(\eta) \right)Cov(\beta,\eta)$

代入 $Var(\widehat{\beta})=1.215$、$Var(\widehat{\eta})=76.42$、$t=20$、$\widehat{\beta}=2.55$ 和 $\widehat{\eta}=39.9$：

$Var(\widehat{u})=0.8513$

$\widehat{u}=-1.762$

得到：

$u_U=-0.2447,\qquad u_L=-3.280$

再将变换尺度转换回不可靠度，得到：

$q_U=0.543,\qquad q_L=0.0369$

对不同时间 $t$ 重复计算，可以把置信限叠加在概率图上。

似然比置信限

似然比原理

似然比置信限以似然比检验为基础。以下以估计数据所拟合的威布尔分布形状参数 $\beta$ 的置信限为例。

似然比统计量为：

$\small\Lambda(\beta) = - 2\ln\left( \frac{L\left( \beta,\widehat{\eta} \right)}{L\left( \widehat{\beta},\widehat{\eta} \right)} \right) = - 2\left\lbrack \ln\left( L\left( \beta,\widehat{\eta} \right) \right) - ln\left( L\left( \widehat{\beta},\widehat{\eta} \right) \right) \right\rbrack$

其中：

• $L(\beta,\widehat{\eta})$：将 $\eta$ 固定为 $\widehat{\eta}$ 时的似然函数。
• $L(\widehat{\beta},\widehat{\eta})$：在估计参数处求得的似然值。

二参数威布尔分布的似然函数为：

$\small L(\beta,\eta) = \prod_{i = 1}^{n}\left\lbrack \left( \frac{\beta}{\eta}\left( \frac{t_{i}}{\eta} \right)^{\beta - 1} \right)^{\delta_{i}} \cdot exp\left( - \left( \frac{t_{i}}{\eta} \right)^{\beta} \right) \right\rbrack$

对于小样本，$\Lambda(\beta)$ 的准确分布未知；根据 Wilks 定理，当样本量较大时，它收敛于自由度等于被检验参数数目的卡方分布。本例只检验一个参数，因此自由度为 1。

在 $100(1-\alpha)\%$ 置信区间和一个自由度下：

$\Lambda(\beta) \leq \chi_{1,1 - \alpha}^{2}$

因此，$\beta$ 的上下限可通过求解满足以下等式的最大和最小 $\beta$ 得到：

$- 2\left\lbrack \ln\left( L\left( \beta,\widehat{\eta} \right) \right) - ln\left( L\left( \widehat{\beta},\widehat{\eta} \right) \right) \right\rbrack = \chi_{1,1 - \alpha}^{2}$

示例

将数据集 {20, 30, 40, 50} 拟合至二参数威布尔分布，采用似然比法求 $\beta$ 的 90% 置信限。

MLE 解为：$\beta=3.57$，$\eta=39$。

求解满足以下等式的 $\beta$ 上下限：

$\small- 2\left\lbrack \ln\left( L(\beta,39) \right) - ln\left( L(3.57,\ 39) \right) \right\rbrack = \chi_{1,0.9}^{2}$

由此得到 $\beta$ 的上下置信限。$\eta$ 的置信限可按相同方法估计。

时间的置信限

考虑相同的数据集 {20, 30, 40, 50}，拟合至二参数威布尔分布。其 B(10) 寿命为 $t=20.8$ 单位时间，即 $Q(t=20.8)=10\%$。现在求 90% 双侧置信限。

以时间 $t$ 表示似然比统计量：

$\small\Lambda(t) = - 2\ln\left( \frac{L\left( t,\widehat{\beta} \right)}{L\left( \widehat{t},\widehat{\beta} \right)} \right) = - 2\left\lbrack \ln\left( L\left( t,\widehat{\beta} \right) \right) - ln\left( L\left( \widehat{t},\widehat{\beta} \right) \right) \right\rbrack$

其中：

• $L(t,\widehat{\beta})$：将 $\beta$ 固定为 $\widehat{\beta}$ 时的似然函数。
• $L(\widehat{t},\widehat{\beta})$：在估计时间和估计 $\beta$ 处求得的似然值；它与 $L(\widehat{\beta},\widehat{\eta})$ 相同。

$L(t,\beta)$ 可通过在 $L(\beta,\eta)$ 中代入下式得到：

$\eta = \frac{\widehat{t}}{\left\lbrack \ln\left( \frac{1}{1 - Q} \right) \right\rbrack^{\frac{1}{\beta}}}$

对于完整数据：

$L(\beta,t) = \prod_{i = i}^{N}\left( \frac{\beta}{\eta} \right)\left( \frac{t}{\eta} \right)^{\beta - 1}e^{- \left( \frac{t}{\eta} \right)^{\beta}}$

其中 $Q=0.1$。求满足以下等式的最大和最小 $t$：

$- 2\left\lbrack \ln\left( L(t,3.57) \right) - ln\left( L(20.8,\ 3.57) \right) \right\rbrack = \chi_{1,0.9}^{2}$