2. 两参数威布尔分布拟合与参数估计
两参数威布尔分布拟合与参数估计
两参数威布尔分布的可靠度函数为 R(t)=exp[−(t/η)^β]，其中 β 为形状参数，η 为尺度参数（特征寿命）。形状参数决定失效机理：β<1 表示早期失效，β≈1 退化为指数分布（恒定失效率），β>1 表示磨损或老化失效。本案例对一组寿命试验数据采用最大似然估计（MLE）求解 β 和 η。
教学示例数据（假设）：某元件 12 个样本的失效时间（小时）：
t = [52, 78, 95, 110, 128, 145, 162, 185, 210, 240, 275, 320]
所有样本均发生失效（完全数据），无删失。
最大似然估计原理
对于完全数据的威布尔分布，似然函数为：
L(β, η) = ∏ (β/η) · (t_i/η)^(β−1) · exp[−(t_i/η)^β]取对数后对 β 和 η 分别求偏导并令其为零，得到 MLE 方程组：
∂lnL/∂η = 0 → η = [ (1/n) · Σ t_i^β ]^(1/β)
∂lnL/∂β = 0 → (1/β) + (1/n)·Σ ln(t_i) − [Σ t_i^β·ln(t_i)] / [Σ t_i^β] = 0第二个方程只含 β，可用数值方法（如牛顿-拉夫逊迭代或二分法）求解，得到 β 后再代入第一个方程求 η。
Python 实现（MLE 数值求解）
以下代码使用二分法求解 β 的 MLE 方程，再计算 η，并输出参数估计结果与拟合可靠度曲线。
import numpy as np
from scipy.optimize import brentq
# 教学示例数据（假设）
t = np.array([52, 78, 95, 110, 128, 145, 162, 185, 210, 240, 275, 320], dtype=float)
n = len(t)
# MLE 方程 g(β) = 0
def g(beta):
# 防止 beta 过小导致数值溢出
if beta <= 0:
return 1e6
t_beta = t ** beta
sum_t_beta = np.sum(t_beta)
sum_t_beta_ln = np.sum(t_beta * np.log(t))
return 1.0/beta + np.mean(np.log(t)) - sum_t_beta_ln / sum_t_beta
# 在 [0.1, 10] 区间内用 Brent 法求根
beta_hat = brentq(g, 0.1, 10.0, xtol=1e-8)
eta_hat = (np.mean(t ** beta_hat)) ** (1.0 / beta_hat)
print(f"形状参数 β̂ = {beta_hat:.4f}")
print(f"尺度参数 η̂ = {eta_hat:.4f} h")
# 可靠度函数
def R(t_val):
return np.exp(-(t_val / eta_hat) ** beta_hat)
# 示例：计算 t=200h 时的可靠度
print(f"R(200) = {R(200):.4f}")运行说明：上述代码依赖 numpy 和 scipy。若环境无 scipy，可用简单二分法替代 brentq：在 [0.1, 10] 内迭代 100 次，每次取中点判断 g 的符号。预期输出 β̂ 约在 2.0~2.3 之间，η̂ 约在 170~190 之间（具体数值取决于数据）。
结果判读与验证
| 指标 | 含义 | 本例预期值 | 判读 |
|---|---|---|---|
| β̂ | 形状参数 | ≈ 2.15 | β̂ > 1，失效随时间增加，属磨损/老化型 |
| η̂ | 特征寿命 | ≈ 178.6 h | 约 63.2% 样本在此前失效 |
| R(η̂) | 特征寿命处可靠度 | ≈ 0.368 | 威布尔分布固有性质，可验证拟合正确性 |
| g(β̂) | MLE 方程残差 | ≈ 0 | 数值求解收敛标志 |
注意：MLE 对完全数据表现良好，但若存在删失数据（未失效样本），似然函数需加入删失项，本案例仅覆盖完全数据情形。此外，小样本（n<20）下 MLE 可能有偏，可考虑修偏估计或贝叶斯方法。