Weibull Life Analysis: Method

Weibull Life Analysis

The Weibull distribution and the shape parameter

The two-parameter Weibull distribution gives the fraction of a population failed by age t as F(t) = 1 − exp[−(t/η)^β], so the reliability is R(t) = exp[−(t/η)^β], the density f(t) = (β/η)(t/η)^β−1 exp[−(t/η)^β] and the hazard rate (instantaneous failure rate of the survivors) h(t) = f/R = (β/η)(t/η)^β−1. The scale η is the characteristic life, the age at which 63.2 % have failed whatever the shape. The shape β sets how the hazard changes with age, which is the physical content of the fit: β < 1 means a falling hazard (infant mortality from defects and damage, the reason for burn-in), β = 1 a constant hazard (random failures, the exponential distribution, where age-based replacement gains nothing), β > 1 a rising hazard (wear-out: fatigue, corrosion, wear), and β above about 4 a narrow life distribution that often signals a mixture of modes, a non-zero time origin or a deterministic process. Weibull's 1951 paper proposed the form for its flexibility; it is also the limit distribution of the weakest link of many, which is why it fits fatigue and fracture so often.

Median ranks and Bernard's approximation

To plot a sample, each failure needs a plotting position: the fraction of the population it represents. For the i-th of n ordered failures that fraction is itself a random variable with a beta distribution, and its median, the median rank, is the F for which the regularized incomplete beta function I_F(i, n − i + 1) equals ½. Bernard and Bos-Levenbach (1953) gave the approximation F ≈ (i − 0.3)/(n + 0.4), accurate to about 0.001 for n of five or more, and it is what the handbooks tabulate; the tool shows both and uses the one chosen. (ReliaWiki's worked example uses the exact ranks and is reproduced to its printed digits with that choice.)

Rank adjustment for suspensions

A suspension is a unit whose life is known only to exceed its stated time: still running when the data were gathered, removed from test, or failed by another mode. It cannot be plotted but it must not be dropped, because leaving survivors out makes every later failure look earlier than it is. Johnson's method assigns each failure a mean order number among the possible orderings of the suspensions: in the increment form of Kececioglu and ReliaSoft, I = (n + 1 − previous order number)/(1 + number of units at or beyond the present failure), order number = previous + I; a run of failures with no suspension between them gets increments of exactly one, so complete data reduce to the ordinary ranks. The median rank of the (generally fractional) order number is then the plotting position. The method uses only the positions of the suspensions, not their exact times, which is the main reason maximum likelihood is preferred when suspensions are many.

Rank regression (least squares on the probability plot)

Taking logarithms twice, ln(−ln(1 − F)) = β ln t − β ln η, so the Weibull plots as a straight line on axes x = ln t, y = ln(−ln(1 − F)); the fitted slope is β and the intercept gives η = exp(−a/b). Rank regression on Y minimizes the vertical deviations of the failure points from the line, the choice of Abernethy's handbook because the time is measured and the rank is the uncertain coordinate; rank regression on X minimizes the horizontal deviations (ReliaSoft's default), giving β = 1/b′ and η = exp(a′) from the line x = a′ + b′ y. The two differ by a factor equal to the squared correlation and coincide for collinear points; the correlation coefficient r of the points with the line is reported as a measure of the fit. With Bernard's ranks the six-point example gives β = 1.427 by RRY and 1.440 by RRX; with the exact ranks 1.430 and 1.443.

Maximum likelihood

The likelihood of the sample is the product of the density at each failure and the reliability at each suspension, L = Π_failures f(t_i) · Π_suspensions R(t_j), and the estimates maximize ln L. Setting ∂ln L/∂η = 0 gives η = (Σ t_i^β/r)^1/β, the sum over all n units and r the number of failures, and substituting it leaves one equation in β, Σ t_i^β ln t_i/Σ t_i^β − 1/β − (1/r) Σ_failures ln t_i = 0, whose left side rises monotonically with β, so bisection finds the single root and a few Newton steps polish it to machine precision. The MLE uses the exact suspension times and works with few failures among many survivors, which rank regression handles poorly; for small complete samples (under about twenty failures) its β is biased high by roughly 10 % at n = 6 and the handbook prefers rank regression there. The tool reports both and lets the primary method be chosen; with suspensions present maximum likelihood is the better default.

Fisher-matrix confidence bounds

The observed Fisher information is the negative matrix of second derivatives of ln L at the estimate; its inverse is the asymptotic variance-covariance matrix of (β, η). Because both parameters are positive, Nelson treats ln β and ln η as normal: the two-sided bounds at confidence C are β exp(∓K √Var(β)/β) and η exp(∓K √Var(η)/η) with K the standard normal quantile at (1 + C)/2. The bounds on a B-life follow from Var(ln t_p) with ln t_p = ln η + ln(−ln(1 − p))/β, and those on R(t) from Var(u) with u = β ln(t/η) and R = exp(−e^u), both by the delta method (first-order propagation of the variances and the covariance). These are large-sample bounds; with few failures they are optimistic, and the likelihood-ratio bounds of Meeker and Escobar are better but not offered here.

B-lives, MTTF and the other life quantities

The B-life B_p is the age by which the fraction p has failed, t_p = η(−ln(1 − p))^1/β; B10 is the bearing industry's L10 and the common warranty figure, B1 the usual aerospace design figure. The mean life is MTTF = η Γ(1 + 1/β) and the median η (ln 2)^1/β; for β near 1 the mean exceeds the median because the distribution is skewed right. The reliable life for a target R is the same formula with p = 1 − R. The expected number of failures by t in a fleet of N new units is N F(t); for a fleet of mixed ages, or for renewal after replacement, it is not, and a separate calculation is needed.

Lognormal and exponential alternatives and AIC

The lognormal, ln t normal with mean μ and standard deviation σ, describes lives governed by a multiplicative process (crack growth, degradation) and its hazard rises and then falls; it is fitted by maximum likelihood with the suspensions entering through their survival terms, 1 − Φ((ln t − μ)/σ), by Newton's method on (μ, ln σ). The exponential is the Weibull with β = 1 and one parameter, λ = r/Σt_i (failures over the total time on test of all units). The three fits are compared by Akaike's information criterion AIC = 2k − 2 ln L with k the number of parameters; the lowest AIC is the best supported model, and a difference under about 2 means the data do not distinguish them, in which case the physics of the failure mode should decide. With a handful of failures the comparison is weak; it becomes informative at twenty or more.

Kaplan-Meier estimate

The product-limit estimate of Kaplan and Meier (1958) assumes no distribution: at each distinct failure time, S is multiplied by (1 − d/n) where n is the number of units still at risk just before that time (neither failed nor suspended earlier) and d the number failing there; between failure times S is constant, and a suspension only reduces the number at risk. Greenwood's formula gives its standard error, se(S) = S √Σ d/[n(n − d)] over the failure times so far. The steps are the data's own survival curve; a fitted model should run through them, and a systematic departure is the plainest sign of a wrong distribution or a mixture of modes.

Probability plotting

The probability plot draws the failures at their median ranks on the linearizing axes, time on a logarithmic scale against ln(−ln(1 − F)), with the fitted lines. The chart helper cannot label the vertical axis in percent failed, so the transformed value is shown; the reference points are F = 1 %, 10 %, 50 %, 63.2 % and 90 % at −4.60, −2.25, −0.37, 0 and 0.83. Curvature of the points indicates a time origin other than zero (a three-parameter Weibull) or a mixture of failure modes (a dog-leg); a straight line through the points with r near one supports the Weibull. The rank-regression line passes through the plotted points by construction; the maximum-likelihood line need not, since it weights the suspensions and the exact times rather than the plotting positions.

Limits

Two-parameter Weibull only: no location parameter (three-parameter Weibull), no interval or left censoring, no grouped inspection data beyond counts at a time, no mixtures or competing modes (analyze each mode separately with the other modes' failures as suspensions), no Bayesian or Weibayes fit with an assumed β, no likelihood-ratio bounds, and no warranty forecasting beyond N F(t) for new units. At least two failures are required; with fewer than five, any estimate of β is weak and the bounds say so.