Reliability — R cross-check (survival::survreg)¶

Companion to reliability_python.ipynb / reliability_pymc.ipynb. Single-sample Weibull MLE + the reliability metrics + a Weibull probability plot.

In [1]:
suppressMessages(library(survival))
d <- read.csv('shock_absorber.csv')
S <- Surv(d$km, d$failed)
wb <- survreg(S ~ 1, data = d, dist = 'weibull')
k <- 1/wb$scale; eta <- exp(coef(wb)[['(Intercept)']])
cat(sprintf('survreg: shape k = %.3f   characteristic life eta = %.0f km\n', k, eta))
cat(sprintf('  B10 = %.0f km   MTTF = %.0f km   R(10000km) = %.3f\n',
            eta*(-log(0.9))^(1/k), eta*gamma(1+1/k), exp(-(10000/eta)^k)))

survreg: shape k = 3.273   characteristic life eta = 25012 km

  B10 = 12576 km   MTTF = 22427 km   R(10000km) = 0.951
In [2]:
# Weibull probability plot (base R)
km <- survfit(S ~ 1, data = d); ev <- km$n.event > 0 & km$surv > 0 & km$surv < 1
x <- log(km$time[ev]); y <- log(-log(km$surv[ev]))
plot(x, y, pch=19, col='steelblue', xlab='ln(distance)', ylab='ln(-ln R)',
     main='Weibull probability plot (shock absorbers)')
abline(coef(wb)[['(Intercept)']]*(-k)*0 + log(1/eta^k), k, col='red', lwd=2, lty=2)
legend('topleft', legend=c('failures (KM)', sprintf('Weibull fit (k=%.2f)', k)),
       pch=c(19,NA), lty=c(NA,2), col=c('steelblue','red'), bty='n', cex=0.85)
No description has been provided for this image

Results¶

survreg matches the Bayesian fits: shape k ≈ 3.27 (wear-out), characteristic life ≈ 25,000 km, B10 ≈ 12,600 km, MTTF ≈ 22,400 km. The probability-plot points fall along the fitted line — Weibull is adequate. Three engines agree (from-scratch ≈ PyMC ≈ survreg); the engineering metrics follow directly from $(k,\eta)$.