Reliability — Weibull life-data analysis (right-censored)¶
Reuses survival_mcmc.py; the engineering framing of the same Weibull¶
Same Weibull machinery as survival_python.ipynb (see it for the censored-likelihood derivation), but the reliability / failure-time framing: a single sample (no covariates) of times-to-failure, most units suspended (still working when the test ended = right-censored). Data: Meeker & Escobar's classic shock-absorber study — 38 vehicles, distance (km) to failure, 23 of 38 suspended.
The questions are distributional, not about a covariate effect:
- Shape $k$ — the failure mode: $k<1$ infant mortality, $k=1$ random, $k>1$ wear-out (the Weibull "bathtub").
- Characteristic life $\eta=\lambda^{-1/k}$ (the 63.2% quantile), $B_{10}$ (distance by which 10% fail), MTTF $=\eta\,\Gamma(1+1/k)$, and reliability $R(t)=S(t)$ at a mission distance.
A note on priors (numerical scaling)¶
With distances in the tens of thousands of km, the Weibull intercept $\beta_0=-k\log\eta\approx-30$ — far from 0 — so a tight $N(0,10^2)$ prior would distort the fit. We use a flat intercept prior (prior_sd=100) and a weakly-informative shape prior (logk_sd=5) so the data drive the estimate; the Bayes fit then matches the MLE. (Equivalently one could rescale time to 1000s of km.)
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import gamma as Gamma
from survival_mcmc import weibull_mle, weibull_rwm, km_estimator
d = pd.read_csv('shock_absorber.csv'); t = d['km'].to_numpy(float); delta = d['failed'].to_numpy(float)
X = np.ones((len(d), 1)) # single sample, intercept only
print('n=%d failures=%d suspended=%d (%.0f%% censored)' % (len(d), int(delta.sum()), int((1-delta).sum()), 100*(1-delta.mean())))
mle, se, _ = weibull_mle(X, t, delta); k_mle = np.exp(mle[1]); eta_mle = np.exp(mle[0])**(-1/k_mle)
o = weibull_rwm(X, t, delta, R=60000, burn=10000, seed=1, prior_sd=100.0, logk_sd=5.0, s=0.6)
dr = o['draws']; k = np.exp(dr[:,1]); lam = np.exp(dr[:,0]); eta = lam**(-1/k)
b10 = (-np.log(0.9)/lam)**(1/k); mttf = eta*Gamma(1+1/k); R10k = np.exp(-lam*10000.0**k)
def ci(a): return f'{a.mean():.0f} [{np.percentile(a,2.5):.0f}, {np.percentile(a,97.5):.0f}]'
print('\nMLE: shape k=%.2f characteristic life eta=%.0f km (acc %.2f)' % (k_mle, eta_mle, o['accept']))
print('Bayes posterior:')
print(' shape k = %.2f [%.2f, %.2f] (k>1 => WEAR-OUT)' % (k.mean(), np.percentile(k,2.5), np.percentile(k,97.5)))
print(' characteristic life = %s km (63.2%% failed)' % ci(eta))
print(' B10 life = %s km (10%% failed)' % ci(b10))
print(' MTTF = %s km' % ci(mttf))
print(' reliability @ 10000 km R(10k) = %.3f [%.3f, %.3f]' % (R10k.mean(), np.percentile(R10k,2.5), np.percentile(R10k,97.5)))
n=38 failures=15 suspended=23 (61% censored)
MLE: shape k=3.27 characteristic life eta=25012 km (acc 0.48) Bayes posterior: shape k = 3.32 [2.22, 4.97] (k>1 => WEAR-OUT) characteristic life = 25552 [21832, 31012] km (63.2% failed) B10 life = 12695 [9171, 16138] km (10% failed) MTTF = 22929 [19584, 27617] km reliability @ 10000 km R(10k) = 0.948 [0.880, 0.989]
In [2]:
fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
# (1) Weibull probability plot: ln(t) vs ln(-ln R); straight line <=> Weibull, slope = k
tk, Sk = km_estimator(t, delta)
ev = (Sk > 0) & (Sk < 1) & (tk > 0)
ax[0].plot(np.log(tk[ev]), np.log(-np.log(Sk[ev])), 'o', color='steelblue', label='failures (KM)')
gx = np.linspace(np.log(t.min()), np.log(t.max()), 50)
ax[0].plot(gx, mle[0] + k_mle*gx, 'r--', lw=1.3, label=f'Weibull fit (k={k_mle:.2f})')
ax[0].set_xlabel('ln(distance)'); ax[0].set_ylabel('ln(-ln R)'); ax[0].set_title('Weibull probability plot'); ax[0].legend(fontsize=8)
# (2) reliability function with 95% band + B10 / characteristic life
grid = np.linspace(1, 40000, 300)
Rg = np.exp(-(lam[:,None]) * grid[None,:]**(k[:,None])) # R(t) per draw
Rm = Rg.mean(0); Rlo = np.percentile(Rg,2.5,axis=0); Rhi = np.percentile(Rg,97.5,axis=0)
ax[1].plot(grid, Rm, color='seagreen', lw=2, label='R(t) posterior mean')
ax[1].fill_between(grid, Rlo, Rhi, color='seagreen', alpha=.2, label='95% band')
ax[1].axhline(0.9, color='grey', ls=':', lw=.8); ax[1].axvline(b10.mean(), color='goldenrod', ls='--', lw=1, label=f'B10 ~ {b10.mean():.0f} km')
ax[1].axhline(np.exp(-1), color='grey', ls=':', lw=.8); ax[1].axvline(eta.mean(), color='firebrick', ls='--', lw=1, label=f'char. life ~ {eta.mean():.0f} km')
# MTTF is the remaining headline quantity and was the only one not on the plot
ax[1].axvline(mttf.mean(), color='rebeccapurple', ls='-.', lw=1, label=f'MTTF ~ {mttf.mean():.0f} km')
# the warranty question the whole fit exists to answer, read off the curve
ax[1].plot([10000, 10000], [0, R10k.mean()], color='k', lw=.8, ls=':')
ax[1].plot([0, 10000], [R10k.mean()]*2, color='k', lw=.8, ls=':')
ax[1].plot(10000, R10k.mean(), 'o', color='k', ms=5,
label=f'R(10,000 km) = {R10k.mean():.3f}')
# where the data actually stop -- everything right of here is extrapolation
ax[1].axvspan(t.max(), 40000, color='grey', alpha=.10, lw=0)
ax[1].text(t.max()+700, .93, 'beyond the\nlast observation', fontsize=7.5, color='#64748b')
ax[1].set_xlabel('distance (km)'); ax[1].set_ylabel('reliability R(t)'); ax[1].set_ylim(0,1.02)
ax[1].set_title('Reliability function'); ax[1].legend(fontsize=7.5, loc='lower left')
print(f'largest distance actually observed : {t.max():,.0f} km')
print(f'of the three life metrics, {(np.array([b10.mean(), mttf.mean(), eta.mean()]) > t.max()).sum()} lie beyond it')
print(f'model-implied reliability at that point: R({t.max():,.0f}) = {np.exp(-lam*t.max()**k).mean():.3f}')
print(f' and well past it: R(40,000) = {np.exp(-lam*40000.0**k).mean():.4f}')
print('\nThat first count is worth pausing on, because the intuition runs the other way. With 61% of')
print('the sample still unfailed you might expect the headline life numbers to be extrapolations --')
print('but all three sit comfortably inside the observed range, below the 28,100 km the longest-')
print('surviving unit reached. What censoring costs here is precision, not reach: every one of these')
print('quantities has a wide interval (B10 runs 9,200-16,100 km) because only 15 units actually')
print('failed, yet each is still interpolation rather than invention.')
print('\nThe shaded region is where that stops being true. Past the last observation the curve is')
print('carried entirely by the fitted shape k = %.2f, and a Kaplan-Meier estimate would simply stop' % k.mean())
print('at the edge of the shading and decline to answer. That is the trade the parametric model makes,')
print('and it is why k -- not the fit to the plotted points -- is the assumption worth arguing about:')
print('R(40,000) above is a statement about a distance no vehicle in this study travelled.')
plt.tight_layout(); plt.savefig('reliability.png', dpi=120, bbox_inches='tight'); plt.show()
largest distance actually observed : 28,100 km
of the three life metrics, 0 lie beyond it
model-implied reliability at that point: R(28,100) = 0.249
and well past it: R(40,000) = 0.0282
That first count is worth pausing on, because the intuition runs the other way. With 61% of
the sample still unfailed you might expect the headline life numbers to be extrapolations --
but all three sit comfortably inside the observed range, below the 28,100 km the longest-
surviving unit reached. What censoring costs here is precision, not reach: every one of these
quantities has a wide interval (B10 runs 9,200-16,100 km) because only 15 units actually
failed, yet each is still interpolation rather than invention.
The shaded region is where that stops being true. Past the last observation the curve is
carried entirely by the fitted shape k = 3.32, and a Kaplan-Meier estimate would simply stop
at the edge of the shading and decline to answer. That is the trade the parametric model makes,
and it is why k -- not the fit to the plotted points -- is the assumption worth arguing about:
R(40,000) above is a statement about a distance no vehicle in this study travelled.
Results¶
- Shape $k\approx3.3$ (CI ~2.2–5.0), well above 1 — a clear wear-out failure mode (hazard rises with distance), exactly what you'd expect for a mechanical part fatiguing over use.
- Characteristic life $\eta\approx25{,}500$ km (63% of units failed by here), $B_{10}\approx12{,}700$ km (the distance by which 10% fail — the usual warranty-style metric), MTTF $\approx22{,}900$ km. Reliability at a 10,000 km mission is high (~0.95).
- The Weibull probability plot (failures fall close to the fitted straight line) confirms the Weibull is an adequate model; the slope is the shape $k$.
- Censoring is essential here: with 61% suspended, ignoring the still-working units would massively understate the life — the censored likelihood uses each suspension as "lasted at least this far."
Takeaway¶
Methodologically this is the same single-sample Weibull as a survival fit with no covariates; the reliability value-add is the interpretation — shape as failure-mode, and the engineering life metrics ($\eta$, $B_{10}$, MTTF, mission reliability) read straight off the posterior. Cross-checks: PyMC and R (survreg) follow.