手搓 Python /JMP搞定可靠性中的 Weibull 数据分

[加可靠性交流群] 

[5万字搞透可靠性工程（上）] 

[往期高赞干货大合集] 

[新朋友必看，本号正确“食用”说明 & 阿超的自我介绍]

[付费精品]

import numpy as npfrom scipy.optimize import minimizefrom scipy.special import gamma# =====================================================================# 1. 录入数据 (包含 确切失效、区间失效、右删失 三种类型)# =====================================================================# 1.1 确切失效时间格式: (发生故障的具体时间 t, 数量 qty)exact_failures = [(450, 1), (520, 1), (680, 1), (750, 1), (810, 1)]# 1.2 区间失效时间格式: (起始时间 t1, 结束时间 t2, 数量 qty)interval_failures = [(400, 600, 2),  # 400到600小时巡检时，发现坏了2个(600, 800, 3),  # 600到800小时巡检时，发现坏了3个(800, 1000, 3)  # 800到1000小时巡检时，发现坏了3个]# 1.3 悬挂/右删失时间格式: (中止测试时产品仍完好的时间 t, 数量 qty)suspensions = [(1000, 3),   # 有 3 个产品跑了 1000 小时没坏，中止了(1500, 14)   # 有 14 个产品跑了 1500 小时没坏，结束测试]# =====================================================================# 2. 定义负对数似然函数 (Negative Log-Likelihood)# =====================================================================def neg_log_likelihood(params):beta, eta = params# 设定物理边界：形状参数和尺度参数都必须大于0if beta <= 0 or eta <= 0:return np.inflog_likelihood = 0.0# 2.1 计算【确切失效数据】的似然 (使用的是 Weibull 概率密度函数 PDF)for t, q in exact_failures:# PDF = (beta/eta) * (t/eta)**(beta-1) * exp(-(t/eta)**beta)pdf = (beta / eta) * ((t / eta)**(beta - 1)) * np.exp(-(t / eta)**beta)if pdf <= 1e-15: return np.inflog_likelihood += q * np.log(pdf)# 2.2 计算【区间失效数据】的似然 (使用的是 F(t2) - F(t1) 累积概率差)for t1, t2, q in interval_failures:# F(t) = 1 - exp(-(t/eta)**beta)p = np.exp(-(t1 / eta)**beta) - np.exp(-(t2 / eta)**beta)if p <= 1e-15: return np.inflog_likelihood += q * np.log(p)# 2.3 计算【悬挂/右删失数据】的似然 (使用的是 可靠度函数 R(t))for t, q in suspensions:# R(t) = exp(-(t/eta)**beta)log_likelihood -= q * (t / eta)**beta# 优化器默认求函数的最小值，所以我们要返回负的对数似然值return -log_likelihood# =====================================================================# 3. 运行优化器求解最佳参数# =====================================================================# 给出初始猜测值 [beta_guess, eta_guess]# 对于常规机械或电子产品，可以粗略猜测 beta=1.5, eta 取大致的试验平均时长initial_guess = [1.5, 2000.0] # 使用 Nelder-Mead 算法寻找极值点 (该算法对无需计算梯度的似然函数非常稳定)result = minimize(neg_log_likelihood, initial_guess, method='Nelder-Mead')beta, eta = result.x# =====================================================================# 4. 计算工程中常用的衍生指标# =====================================================================# B10 寿命：预期有 10% 产品发生失效的时间b10 = eta * (-np.log(0.9))**(1/beta)# MTTF 平均寿命：期望寿命mttf = eta * gamma(1 + 1/beta)# =====================================================================# 5. 打印结果对齐商业软件# =====================================================================print("=" * 45)print("📊 混合数据 Weibull 极大似然估计 (MLE) 结果")print("=" * 45)print(f"优化状态        : {'✅ 成功'if result.success else'❌ 失败'}")print(f"形状参数 (Beta) : {beta:.3f}")print(f"尺度参数 (Eta)  : {eta:.1f} (在 JMP 中也称为 Alpha)")print("-" * 45)print(f"B10 寿命        : {b10:.1f} 小时")print(f"MTTF (平均寿命) : {mttf:.1f} 小时")print("=" * 45)