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')
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=8)
plt.tight_layout(); plt.savefig('reliability.png', dpi=120, bbox_inches='tight'); plt.show()
No description has been provided for this image

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.