手里攥着一批部件的运行时长记录四十个样本一大半已经失效剩下十几个到测试结束还在正常跑这时候新人最常干的事就是把这些记录一平均报一个平均寿命上去。我见过太多次这样的报表被推翻原因不是算错了加减乘除而是这批数据里混着两种完全不同性质的信息一种是确切的死亡时间另一种是我只知道它活到了这个点之后的事我不清楚。把后者当成前者参与算术平均等于主动往低估方向拉了一把后面的所有结论都建在流沙上。威布尔分布就是用来收拾这类局面的主力工具它把生存时间数据的非负性、长尾特征和删失结构一次性装进一个只有两个参数的框架里还能顺手把故障是早期磨合问题还是累积磨损问题这个机理信号从数据里读出来。这篇内容写给已经在跑寿命试验、做可靠性评估或者分析用户流失时间的同行也写给刚接触生存分析、被指数分布拟合出一条直线的残差图搞得一头雾水的初学者我会把参数化约定、似然函数怎么搭、B10寿命怎么算、哪些坑反复有人踩一条一条摊开讲。1. 生存时间数据的三处硬骨头删失、非负、长尾1.1 右删失把算术平均这条路堵死了先把这个概念钉死。右删失指的是观测在某个时间点结束时对象还没有发生我们关心的事件。部件还在转用户还在用受试者还健在。这时候你手里唯一的确定信息是它的生存时间大于等于这个观测时长而不是它的事件时间等于这个观测时长。这两种解读的差别有多大举个极端例子一百个样本九十个在第 50 小时失效十个跑到第 500 小时还没坏。真实平均寿命一定高于 50 小时如果把这十个样本按 500 小时算事件时间正常情况下会高估可如果干脆把它们从分母里剔掉、只算那九十个你得到的 50 小时就是个明显的低估因为存活到 500 小时的那些大概率还能继续跑。两种粗暴处理都会把结论带偏方向还不一样。正确的做法是从生存函数 S(t) P(T t) 出发让每个样本按它贡献的信息量进入似然。失效样本贡献的是概率密度 f(t)删失样本贡献的是生存概率 S(t)。这一步转换看着只是一行公式实际是整个生存分析的地基删失样本不再是缺失数据而是部分信息它带着约束进来把参数估计往正确的方向推。我一般会让新人先手推一遍两参数威布尔的似然函数推到能自己写出δ标记怎么进指数为止之后就很少犯删失当事件这类低级错误了。还有一类更麻烦的区间删失你只知道事件发生在两次巡检之间确切时间未知。工控设备的定期检查、体检随访都属于这种。两参数威布尔配合区间似然同样能处理只是似然项变成两点 CDF 之差代价是需要写一个自定义的负对数似然函数现成的库未必直接支持。这块如果样本量不大我倾向于退回离散化处理或者干脆用参数化 bootstrap 跳过解析式。1.2 指数分布为什么总在关键时刻不够用指数分布是威布尔在形状参数等于 1 时的特例它的核心特征是风险率恒定。也就是说一个用了几千小时的部件下一小时失效的概率和一个刚上机的新件完全一样。这个假设在某些场景下勉强站得住比如随机过载击穿、偶发外力冲击但凡是涉及材料疲劳、润滑劣化、触点磨损、电池衰减的场合它几乎必然错。错在哪里错在数据里明明有趋势你硬套一个无趋势模型后果不是精度差一点而是残差里留下系统性结构。我在早期做过一次轴承寿命拟合用指数分布跑完概率图上点子在低分位端向上弯、高分位端向下弯形成一个明显的弧形。当时没在意直接报了参数后来拿外部验证集一跑实际早期失效率远高于模型预期。换成威布尔形状参数落在 2.3 附近风险率随时间单调上升概率图立刻拉直。那次之后我形成了一个习惯任何生存数据先用威布尔跑一遍看形状参数它落在 [0.9, 1.1] 区间内我再考虑退回指数模型简化落在这个区间外指数模型直接出局。威布尔之所以能顶上来还因为它和极值理论有血缘关系一个系统由多个串联环节组成系统失效时间取各环节失效时间的最小值这类最弱环结构在弱条件下会趋近威布尔分布。这个背景解释了一件事——威布尔在工程可靠性里适用面这么宽不是巧合背后有结构性的道理。1.3 两个参数各自在管什么两参数威布尔的生存函数写出来是 S(t) exp(-(t/λ)^β)累积分布函数 F(t) 1 - exp(-(t/λ)^β)概率密度 f(t) (β/λ)(t/λ)^(β-1) exp(-(t/λ)^β)风险函数 h(t) (β/λ)(t/λ)^(β-1)。四个式子看着差不多但真正做工程决策时天天用的是风险函数因为它直接对应现在这个时点下一单位时间出问题的强度。尺度参数 λ 管时间轴的横向拉伸量纲和时间一致含义是特征寿命t λ 时S(λ) exp(-1) ≈ 0.368也就是大约 63.2% 的样本已经失效。注意它不等于平均寿命这一点被搞混的频率高得吓人。形状参数 β 无量纲管风险率的走向是纯粹的机理参数β 1 风险递减β 1 恒定β 1 递增β 2 时风险率是一条过原点的直线恰好对应瑞利分布。这里有一个必须先说清楚的问题参数化约定不统一。有的教材和软件用 rate 参数化的形式 f(t) βλt^(β-1)exp(-λt^β)这里的 λ 量纲是时间的 -β 次方数值上等于尺度参数的 -β 次方另一些用 scale 参数化λ 就是特征寿命。这两种写法的 λ 完全不是一回事混用会得到一个荒谬的结果比如算出特征寿命 0.0003 小时之类的。套任何现成代码之前去文档里确认它用的是哪一种或者干脆拿一个已知答案的小样本喂进去验一下这一步花三分钟能省掉后面三小时。2. 形状参数 β 读出来的失效机理信号2.1 β 小于 1早期失效段的典型脸孔β 落在 0.5 到 0.9 之间风险率随时间单调下降意味着越新的东西越容易坏。听起来反直觉但制造业里这太常见了焊接虚接、装配应力、材料内部缺陷、工艺参数漂移造成的批次性问题都会在投入使用初期集中暴露。跑过一段时间的存活件反而越来越稳因为它们已经筛掉了。这类数据的处置思路和磨损失效完全不同。你不需要延长保修期或者加大备件库存你要做的是把早期失效段的缺陷件筛出来做老化筛选、加严出厂检验、追查工艺环节。我见过一个团队拿着 β 0.7 的数据去论证应该缩短更换周期方向整个反了按那个模型缩短周期反而把设备推进了高风险的早期区间维护成本上升故障率却没降。这里还有一个坑β 1 时密度函数在 t 0 处发散。有些数值优化器在这个区域会抖初值给不好直接发散到边界。我一般的处理是给 β 加一个下界约束比如 0.1并且在似然里对 t 0 的样本做单独处理或者干脆把它们当作左删失。2.2 β 接近 1别急着下结论说随机失效β 估计值在 1 附近的 10% 以内很容易被解释成失效是随机的没有老化趋势。这个结论有时候成立有时候是样本量不够导致区间太宽——β 的 95% 置信区间可能是 [0.7, 1.4]你把点估计 1.02 当成真值在讲其实是统计精度不足在骗自己。所以看到 β ≈ 1第一件事是去看它的置信区间宽度而不是直接下判断。第二件事是看数据结构。如果数据本身是异质总体混在一起的比如一批设备由两个供应商供货一个供应商的件有磨损趋势另一个是随机失效混在一起做单峰拟合形状参数会被拉向中间。这种混合分布伪装成随机失效的情况用单参数威布尔永远看不出来得上混合模型或者做分层分析。真正意义上的随机失效在工程里对应的是过载、误操作、外部冲击这些与累积使用时间无关的机理。如果你确认了机理确实是这一类那么指数模型的简化就是合理的而且能省掉一个待估参数在小样本下是实打实的收益。2.3 β 大于 1累积损伤的量化指标β 超过 1风险率单调上升这是磨损失效、疲劳、腐蚀、绝缘老化的标准特征。β 越大风险率上升越陡失效越集中在特征寿命附近。β 在 2 到 4 之间是最常见的磨损区间β 大于 5 的情况说明失效高度集中往往是单一主导机理在起作用。这里有个工程上很好用的推论β 大于 1 时预防性更换是有意义的而且最优更换间隔可以直接从模型里算出来。β 越陡这个窗口越窄意味着更换时间点必须卡得比较准早换浪费寿命晚换故障率飙升。相反β 小于 1 时定期更换策略在数学上就是低效的这时候该做的是筛选和 burn-in不是维护计划。β 区间风险率走势常见机理该做的动作0.3 – 0.8快速下降装配缺陷、材料夹杂、工艺波动出厂筛选、老化跑合、追工艺0.8 – 1.2基本恒定随机过载、外部冲击、误操作简化成指数模型关注备件1.2 – 2.5缓慢上升轻度磨损、润滑劣化定期巡检可做预防更换2.5 – 5明显上升疲劳、腐蚀、绝缘老化严格按周期更换监控趋势 5陡峭上升单一主导机理接近定时失效卡点更换容错设计这张表是我自己项目里常用的速查表但必须强调它是经验区间不是定律机理判断永远要回到失效分析报告上数据只能提供线索不能替代拆解。3. 参数估计的三条路线各有各的适用场景3.1 极大似然含删失样本下的标准解法两参数威布尔在含右删失样本时的对数似然是这个形状log L(β, λ) Σ_i δ_i · [ln β - ln λ (β-1)(ln t_i - ln λ)] - Σ_i (t_i/λ)^β其中 δ_i 1 表示观测到失效δ_i 0 表示右删失。注意删失样本只进入第二项它对参数的约束是这个时间点之前别失效体现在惩罚那些让 S(t_i) 太小的参数组合上。这个式子对 λ 有解析解对 β 没有实际求解用二维数值优化或者先对 λ 做剖面化再对 β 一维搜索。初值怎么给β 给 1.0 或者用概率图法的斜率λ 给样本中位数或者所有观测时间的几何均值。这两组初值在绝大多数情况下都能收敛。如果优化器报奇异 Hessian 或者参数跑到边界八成是数据里有 t 0 的记录或者删失标记全为 0那就是纯粹的参数拟合没有生存结构了。import numpy as np from lifelines import WeibullFitter # durations: 观测时长event_observed: 1失效, 0右删失 durations np.array([...]) event_observed np.array([...]) wf WeibullFitter() wf.fit(durations, event_observedevent_observed) print(shape (rho):, wf.rho_) print(scale (lambda):, wf.lambda_) print(log-likelihood:, wf.log_likelihood_)用 lifelines 的时候要留意它的参数名rho_是形状参数lambda_是尺度参数生存函数定义成 S(t) exp(-(t/λ)^ρ)和前面讲的 scale 参数化一致。不同版本对属性名的处理有过调整代码跑不通先去查你装的版本别急着怀疑数据。如果不想依赖高层封装用 scipy 自己写负对数似然也很直接好处是完全掌控参数化也方便加约束和自定义删失类型from scipy.optimize import minimize from scipy.stats import weibull_min def neg_log_lik(params, t, d): beta, lam params if beta 0 or lam 0: return 1e12 ll np.sum(d * (np.log(beta) - np.log(lam) (beta - 1) * (np.log(t) - np.log(lam))) - (t / lam) ** beta) return -ll res minimize(neg_log_lik, x0[1.0, np.median(durations)], args(durations, event_observed), methodNelder-Mead, options{xatol: 1e-8, fatol: 1e-8})我通常两种都写互为交叉验证。两个实现给出的 β 差超过 1%说明中间某一环的参数化或者删失标记出了问题这时候停下来查比继续往下跑划算得多。3.2 概率图法现场没有计算环境时的快刀概率图法的思路是把分布函数做两次对数变换把曲线拉成直线。对 F(t) 1 - exp(-(t/λ)^β) 取变换y ln(-ln(1 - F(t))) β · ln(t) - β · ln(λ)于是 y 对 x ln(t) 就是一条直线斜率就是 β截距是 -β ln λ。实操时横轴用 ln(t_i)纵轴用经验分布函数的估计值代入 y。经验分布怎么估样本量小的时候用中位秩公式 F_i ≈ (i - 0.3) / (n 0.4)比 i/n 偏差小样本量大用哪个都无所谓。这里有个细节容易忽略排序和删失的关系。如果存在右删失不能简单地按时间排序取秩次必须用 Kaplan-Meier 估计或者专门的秩回归方法处理否则删失点会被错误地分配到一个偏大的累积概率上整条直线被拉歪。import numpy as np t_sorted np.sort(durations) n len(t_sorted) i np.arange(1, n 1) F_hat (i - 0.3) / (n 0.4) x np.log(t_sorted) y np.log(-np.log(1 - F_hat)) slope, intercept np.polyfit(x, y, 1) beta_hat slope lambda_hat np.exp(-intercept / slope)这个方法我到现在还在用原因是它极其直观把点画上去肉眼一看线性度就知道威布尔适不适合弯曲的方向还能告诉你该往哪个分布上试。左端翘起、右端压平往往提示数据里有早期失效的子总体右端翘起则可能是磨损失效的尾部在作怪。这种判断速度是数值拟合给不了的。当然它的精度不如极大似然尤其在删失比例高的时候因为删失点没有被充分利用。我的做法是拿它出初值和做直观判断最终参数报极大似然的结果。3.3 小样本和先验信息贝叶斯这条路的性价比样本量掉到个位数或者十几个的时候极大似然估计会变得很不稳定β 的点估计可能今天 1.8 明天 3.2置信区间宽得没有实际意义。这时候引入先验是有实打实收益的尤其当你手上有同系列产品、同类工艺的历史数据。具体做法是在 β 和 λ 上放对数正态先验或者 Gamma 先验β 的先验中心根据历史机理判断给比如同类磨损件历史 β 在 2 附近就给一个以 2 为中心、对数标准差 0.4 的弱信息先验。这个先验强度不能太大否则数据说了不算等于把历史结论硬塞进新数据。我的经验是样本量低于 20 时先验带来的收缩是帮助样本量超过 50先验影响基本被数据淹没放不放差别不大。代价是计算复杂度和沟通成本。贝叶斯给的是后验分布报给业务方的时候要解释清楚这个区间是可信区间不是置信区间很多人到这一步就绕晕了。所以我的选择是常规项目走极大似然小样本或者需要融合历史信息的场景走贝叶斯两条路的结果碰一下差异大就说明数据信息量不足结论要谨慎表述。4. 拟合质量怎么判光看相关系数会翻车4.1 对数似然和 AIC 的正确用法判断威布尔拟合得好不好最常见的错误是拿概率图上的 R² 说事。R² 高只说明点子在 log-log 坐标下接近直线但这条直线的位置对不对、尾部拟合准不准它一点都给不了。而且 R² 会随样本量增加自动变高跟模型优劣没有必然关系。我的判定顺序是先看对数似然值再做同维度模型的 AIC 对比最后回到概率图看残差结构。AIC 2k - 2lnLk 是参数个数两参数威布尔 k 2。拿它和对数正态、Gamma、指数模型比差值小于 2 说明区分度不够别硬说谁更好差值超过 10 才有明确的倾向性。from lifelines import WeibullFitter, ExponentialFitter, LogNormalFitter models { Weibull: WeibullFitter(), Exponential: ExponentialFitter(), LogNormal: LogNormalFitter(), } for name, m in models.items(): m.fit(durations, event_observedevent_observed) print(name, AIC , m.AIC_partial_)要提醒的是这几个模型的参数个数不同指数模型只有一个参数AIC 对比时它天然占便宜一点所以当威布尔和指数的 AIC 接近时不能直接选指数得回到机理上看风险率到底是不是恒定的。4.2 尾部偏差为什么会被放大生存分析里最有价值的信息往往在尾部B1 寿命、B0.1 寿命这类极端分位点恰恰是工程上最关心的。但尾部也是拟合误差被放大最厉害的地方因为那里样本稀少一两个点的位置变化就能让曲线大幅摆动。具体表现是概率图中间段拟合得漂漂亮亮两端的点却系统性偏离。左端偏离说明早期失效段建模不准通常是混合总体没拆开右端偏离说明磨损尾部的形状和威布尔假设不符可能是失效机理在后期发生了切换。这两种情况都不是换个优化算法能解决的必须回到数据本身去找子总体。我做尾部评估时会专门算几个低分位点的拟合值和 Kaplan-Meier 经验值的差比如 B1、B5、B10看偏差是否在置信带内。这一步比看整体 AIC 更能暴露问题。4.3 和竞争分布的取舍逻辑威布尔不是唯一选项对数正态和 Gamma 在很多场景下表现相当甚至更好。三者的风险率行为差异明显威布尔的风险率是幂函数单调且没有上界对数正态的风险率先升后降存在一个峰值Gamma 的风险率从某个非零值起步再上升。如果你的数据里风险率明显在某个时间点之后开始回落威布尔是拟合不出来的硬拟合会在高分区产生系统偏差。我的经验判断是失效由累积损伤主导、风险率单调上升优先威布尔失效由多因素叠加、早期有个爬升过程试试对数正态需要更灵活的形状控制用 Gamma 或者广义 Gamma。实际项目里我一般三个都跑AIC 差距不明显就选机理上说得通的别为了小数点后的似然值去选一个解释不了的模型。5. 从参数到工程指标B10 寿命和置信区间的实操5.1 MTTF 和 B 寿命的换算拿到 β 和 λ 之后最常被问到的是两个数平均寿命和 B10 寿命。平均寿命 MTTF λ · Γ(1 1/β)其中 Γ 是伽马函数。这一步很多人直接用 λ 当平均寿命报上去在 β 2 的时候误差接近 11%在 β 3 的时候误差超过 10%在小 β 时误差更夸张。所以务必走一遍 Γ 函数代码里就是scipy.special.gamma(1 1/beta)。B10 寿命指累积失效概率达到 10% 对应的时间解 F(t) 0.1 得 t_B10 λ · (-ln 0.9)^(1/β) ≈ λ · 0.10536^(1/β)。同理 B1 λ · 0.01005^(1/β)。这类低分位寿命是可靠性指标里最常被引用的因为它对应的是早期批次风险比平均值更能反映用户实际感受。举个具体例子某批液压密封件n 40测试结束时 28 个失效、12 个右删失极大似然给出 β 2.1、λ 1800 小时。那么MTTF 1800 × Γ(1 1/2.1) 1800 × 0.8855 ≈ 1594 小时B10 1800 × 0.10536^(1/2.1) 1800 × 0.3425 ≈ 617 小时B1 1800 × 0.01005^(1/2.1)0.01005^(0.4762) ≈ exp(-4.6002 × 0.4762) exp(-2.1906) ≈ 0.1119得约 201 小时B1 和 MTTF 差了将近八倍这个差距在跟业务方沟通时非常关键。你要是只报平均寿命对方会以为产品能用一千五百小时实际上有 1% 的件在两百小时就挂了售后投诉全集中在那批。5.2 置信区间三种获取方式的选择点估计单独报出去是不负责任的尤其样本量不大的时候。三个常用路径Fisher 信息矩阵法用负对数似然在最优点的 Hessian 求逆得到协方差再取对角元素开方。速度快但依赖大样本正态近似样本小或者 β 接近边界时区间会失真甚至出现负的下界。我会在对参数做对数变换后再求 Hessian保证区间恒正然后用 delta 方法换算回原尺度。似然剖面法固定一个参数扫另一个找对数似然下降 1.92对应 95% 单参数的位置。这个方法在小样本下表现好得多也能处理不对称的区间形状。代价是计算量大一点但现在这点算力根本不算事。我现在基本默认用剖面法。Bootstrap对样本做有放回重抽样每次重新拟合取分位数。优点是几乎不依赖任何分布假设实现简单还能顺带给出 B10 的置信区间——这一点前两种方法做不到因为 B10 是参数的非线性函数。缺点是对删失数据的重抽样需要注意简单的样本重抽样可能把删失结构打乱稳妥的做法是做条件 bootstrap 或者用参数 bootstrap从拟合好的威布尔里生成新样本保留原来的删失时间结构。rng np.random.default_rng(42) boot_b10 [] for _ in range(2000): idx rng.integers(0, len(durations), len(durations)) d_b durations[idx] e_b event_observed[idx] if e_b.sum() 2: continue m WeibullFitter().fit(d_b, event_observede_b) if not np.isfinite(m.rho_): continue boot_b10.append(m.lambda_ * 0.10536 ** (1 / m.rho_)) lo, hi np.percentile(boot_b10, [2.5, 97.5])这段代码我基本是复制粘贴改改就用的2000 次重抽样在普通机器上几秒钟跑完。注意两个跳过条件失效样本太少、拟合参数非有限值这两种情况在重抽样里出现得比想象中频繁不处理会让分位数被 NaN 污染。5.3 把结论翻译成维护动作参数算完只是中间产物真正交付的是决策。β 2.1 这个结果对应的动作是失效由累积磨损主导定期更换有效B10 617 小时意味着如果要控制 10% 的早期失效率更换窗口应该卡在 600 小时之前MTTF 1594 小时说明平均表现还行但尾部风险集中在早期所以策略重心应该放在首检而不是后期维护。如果算出来是 β 0.7结论就完全反过来更换周期再密也压不住早期故障该做的是加严出厂筛选、做 50 小时的跑合筛选把缺陷件提前淘汰。同一套数据、同一个模型形状参数不同投入的方向就不同这也是我一直强调 β 是机理参数而不是普通拟合系数的原因。6. 反复踩到的几个坑写下来省点时间6.1 删失标记写反是最贵的错误这个错误我见过不止一次而且极难自己发现因为模型照常收敛、参数照常输出、图照常漂亮只是结论完全反了。原因在于不同的库对 event 标记的约定不一致有的用 1 表示失效有的用 True 表示观测到事件还有的用 0 表示删失。一旦反了模型会把大部分样本当成活到了很久把删失的少数当成很快就坏了形状参数和尺度参数都会跑偏到一个看起来很合理的假值上。防守办法很简单跑完先做一次合理性检查拿 Kaplan-Meier 估计的经验中位生存时间和你拟合出来的模型中位生存时间比一比差得离谱就说明有环节错了。另外把event_observed.sum()打出来跟你原始记录里的失效数量对一下这个数字对不上别的都不用看。6.2 位置参数该不该上三参数威布尔多了一个位置参数 γ形式是 S(t) exp(-((t-γ)/λ)^β)物理含义是在 γ 之前绝对不会失效也叫保证寿命或门槛值。理论上很吸引人实际上是个大坑。原因在于 γ 的估计极不稳定它和 β、λ 高度相关似然曲面在 γ 方向常常很平坦样本量不够时估计值会贴着最小值乱跑。更麻烦的是如果真实分布其实没有门槛很多失效机理确实没有硬拟合三参数会让另外两个参数的方差大幅膨胀得不偿失。我的原则是除非有明确的物理证据支持存在一个无失效区间比如某种涂层在厚度耗尽之前完全不失效否则坚持两参数。真要用也先固定 γ 做敏感性分析看结论对 γ 的取值敏不敏感敏感就别用。6.3 混合总体把单峰模型骗了一批产品来自两个工艺不同的批次或者同一个系统有多个失效模式数据里就会出现两个不同尺度的子总体。单峰威布尔拟合这种数据卡在中间形状参数被拉平尾部严重失真概率图上表现为一段一段的偏离而不是整体弯曲。识别方法是画经验风险率的非参数估计或者直接看概率图的残差模式。看到分段特征就该考虑混合威布尔或者竞争风险模型了。做这一步之前先确认数据来源很多混合总体其实是数据记录口径不一致造成的比如某个批次的计时单位不同、或者起始时间点定义不同先排除这个再上复杂模型。6.4 外推这件事要克制威布尔在观测区间内拟合得好不意味着往观测范围外推还准。特别是往远尾外推形状参数的微小不确定性会被指数放大。我用 1000 小时的数据去预测 10000 小时的失效率β 的 95% 区间是 [1.9, 2.3]那么外推出的高分量预测区间会宽到没有实用价值。诚实的做法是把外推的不确定性显式画出来用参数 bootstrap 生成预测带让对方看到这个外推值可能在很大一个范围内。比起到处给一个看起来精确的单点预测我宁愿报一个宽区间至少不会误导后续的备件和保修决策。最后分享一个我自己在项目里养成的小习惯每次做完威布尔拟合我都会在报告首页画一张概率图把原始点、拟合直线和置信带叠在一起然后在旁边用一行字写清楚形状参数的机理判读。目的是让看报告的人第一眼就能判断这个拟合靠不靠谱、结论能不能用而不是被一堆参数表淹没。这张图被追问的次数比所有数字加起来都多。