A technical guide to Maximum Likelihood Estimation, Weibull Probability Plotting, Method of Moments, and Python implementation for reliability engineering and survival analysis.
Analyzing when a jet engine turbine blade will crack, how long a customer subscription will last before churning, or how fast the wind blows at an offshore turbine site frequently leads to the Weibull distribution.
Named after Swedish mathematician Waloddi Weibull, who popularized it in 1951 to describe the breaking strength of materials, the Weibull distribution is the foundational probability distribution of survival analysis and reliability engineering.
By adjusting just two parameters—the shape parameter (k) and the scale parameter (λ)—the Weibull distribution can represent an exponential distribution, approximate a normal distribution, or model heavily skewed failure rates across every stage of a system’s lifecycle.
However, given a raw array of time-to-failure observations t1, t2, …, tn, those two parameters cannot be calculated by simply plugging the sample mean and standard deviation into a basic closed-form formula.
This article breaks down the physical intuition behind the Weibull parameters, derives the three primary mathematical methods used to estimate them (including adjustments for right-censored data), and provides a complete Python implementation.
1. The Anatomy of the Weibull Distribution
For a continuous random variable T > 0 (typically representing time-to-failure), the two-parameter Weibull Probability Density Function (PDF) and Cumulative Distribution Function (CDF) are defined as:
f(t; k, λ) = (k / λ) ⋅ (t / λ)k−1 ⋅ e−(t/λ)k
F(t; k, λ) = 1 − e−(t/λ)k
The Shape Parameter (k > 0) — The Weibull Slope (β)
The shape parameter dictates the failure mechanism via the hazard function h(t) = (k / λ) ⋅ (t / λ)k−1, modeling the three phases of the classic Bathtub Curve of reliability:
- k < 1 (Infant Mortality): The failure rate decreases over time. Defective units fail early, while surviving units prove more reliable (e.g., manufacturing defects in electronics or early subscription churn).
- k = 1 (Random Failures): The Weibull distribution simplifies to the Exponential distribution. The hazard rate remains constant (1 / λ), meaning the component’s failure probability is independent of its age.
- k > 1 (Wear-out / Aging): The failure rate increases over time as components degrade from mechanical fatigue, corrosion, or friction. At k = 2, it forms the Rayleigh distribution; around k ≈ 3.6, it closely approximates a Normal distribution.
The Scale Parameter (λ > 0) — The Characteristic Life (η)
The scale parameter stretches or compresses the distribution along the time axis. Evaluating the CDF at t = λ yields:
F(λ; k, λ) = 1 − e−(λ/λ)k = 1 − e−1 ≈ 0.6321
Regardless of the shape parameter k, λ represents the exact time at which 63.2% of the population is expected to fail.
2. Three Methods for Estimating k and λ
Selecting an estimation method depends on sample size, computational requirements, and the presence of censored observations or extreme outliers.
| Estimation Method | Statistical Efficiency | Computational Effort | Primary Use Case |
|---|---|---|---|
| 1. Maximum Likelihood (MLE) | Highest (Asymptotically unbiased & minimum variance) | Iterative numerical root-finding (Newton-Raphson) | Medium-to-large samples (n > 30) & censored survival data |
| 2. Linearized Rank Regression (OLS) | Moderate (Sensitive to extreme tail points) | Closed-form linear regression | Small samples (n < 20), visual diagnostics (Weibull plots) |
| 3. Method of Moments (MoM) | Low-to-Moderate (High variance when k < 1) | Numerical inversion of Gamma function ratios | Quick initialization values for MLE solvers |
Method 1: Maximum Likelihood Estimation (MLE)
Given n independent failure times t1, t2, …, tn, the likelihood function L(k, λ) is the product of the individual PDFs:
L(k, λ) = ∏i=1n (k / λ) ⋅ (ti / λ)k−1 ⋅ e−(ti/λ)k
Taking the natural logarithm yields the Log-Likelihood function ℓ(k, λ):
ℓ(k, λ) = n ln(k) − nk ln(λ) + (k − 1) ∑i=1n ln(ti) − ∑i=1n (ti / λ)k
To find the estimators k̂ and λ̂ that maximize ℓ(k, λ), the partial derivatives with respect to λ and k are set to zero.
Step A: Differentiate with respect to λ
∂ℓ / ∂λ = −(nk / λ) + (k / λk+1) ∑i=1n tik = 0
Solving for λ produces an explicit closed-form expression for λ̂ in terms of k̂:
λ̂ = ( (1 / n) ∑i=1n tik̂ )1/k̂
Mathematically, λ̂ is the generalized power mean of order k̂.
Step B: Differentiate with respect to k and Substitute λ̂
Setting ∂ℓ / ∂k = 0 and substituting the expression for λ̂k eliminates λ, reducing the system to a single nonlinear equation in k:
g(k) = [ ∑i=1n tik ln(ti) / ∑i=1n tik ] − (1 / k) − (1 / n) ∑i=1n ln(ti) = 0
Because k appears inside both a power summation and a reciprocal, no algebraic closed-form solution exists. Standard statistical packages solve g(k) = 0 numerically using 1D Newton-Raphson iteration (km+1 = km − g(km) / g′(km)), then substitute the converged k̂ into Step A to obtain λ̂.
Method 2: Weibull Probability Plotting & Rank Regression (OLS)
Before iterative numerical solvers were ubiquitous, reliability engineers linearized the Weibull CDF to plot failure data on logarithmic probability paper. This technique—Rank Regression—remains widely used for small sample sizes (n < 20) and visual goodness-of-fit diagnostics.
Starting from the survival function 1 − F(t) = e−(t/λ)k, a double-logarithmic transformation linearizes the equation:
- Take the reciprocal and natural logarithm once: ln(1 / (1 − F(t))) = (t / λ)k
- Take the natural logarithm a second time:
ln(−ln(1 − F(t))) = k ln(t) − k ln(λ)
This expression maps directly to a linear equation Y = mX + c:
- Y = ln(−ln(1 − F(t)))
- X = ln(t)
- Slope (m) = k
- Intercept (c) = −k ln(λ) ⇒ λ̂ = exp(−c / k̂)
Estimating F(ti) via Median Ranks
When failure times are sorted in ascending order t(1) ≤ t(2) ≤ … ≤ t(n), the naive empirical CDF i / n cannot be used because at i = n, ln(1 − 1) = ln(0) = −∞.
Instead, Bernard’s Median Rank Approximation is used to estimate the cumulative probability of the i-th ordered failure:
F̂(t(i)) = (i − 0.3) / (n + 0.4)
Once Xi = ln(t(i)) and Yi = ln(−ln(1 − F̂(t(i)))) are computed, Ordinary Least Squares (OLS) linear regression provides k̂ and λ̂. Furthermore, a high coefficient of determination (R2) on this plot visually confirms that the data follows a two-parameter Weibull distribution.
Method 3: Method of Moments (MoM)
The Method of Moments equates theoretical population moments to sample moments. For a Weibull distribution, the r-th raw moment about zero is:
피[Tr] = λr Γ(1 + r / k)
where Γ(z) = ∫0∞ uz−1e−udu is the Gamma function. Consequently, the theoretical mean (μ) and variance (σ2) are:
μ = λΓ(1 + 1 / k), σ2 = λ2 [ Γ(1 + 2 / k) − Γ2(1 + 1 / k) ]
Dividing the standard deviation σ by the mean μ cancels out λ, leaving the Coefficient of Variation (CV) as a function of k alone:
CV = σ / μ = √[ Γ(1 + 2 / k) − Γ2(1 + 1 / k) ] / Γ(1 + 1 / k)
Replacing μ and σ with the sample mean t̄ and sample standard deviation s allows solving for k̂ numerically (or via analytical approximations such as k̂ ≈ (s / t̄)−1.086 for k > 1), after which the scale parameter is recovered via λ̂ = t̄ / Γ(1 + 1 / k̂).
3. Handling Right-Censored Data
In industrial life testing and clinical survival studies, tests often terminate before 100% of units fail. Units that remain operational at the end of the observation window are right-censored.
Discarding right-censored observations biases λ̂ downward, while treating censoring times as actual failures distorts k̂.
Maximum Likelihood Estimation accommodates right-censoring directly. Let δi = 1 if unit i failed at time ti, and δi = 0 if unit i was right-censored at time ti. With r = ∑i=1n δi representing the total number of observed failures, the censored MLE equations become:
λ̂ = ( (1 / r) ∑i=1n tik̂ )1/k̂
[ ∑i=1n tik ln(ti) / ∑i=1n tik ] − (1 / k) − (1 / r) ∑i=1n δi ln(ti) = 0
The power sums of tik run across all n units (accounting for total accumulated exposure time), whereas the divisor r and the unweighted logarithmic sum count only observed failures (δi = 1).
4. Python Implementation
The following script simulates a Weibull dataset with true parameters k = 2.2 and λ = 100.0, comparing Newton-Raphson MLE, Rank Regression (OLS), and scipy.stats.weibull_min.
Implementation Note on
scipy: By default,scipy.stats.weibull_min.fit(data)estimates a 3-parameter Weibull distribution (including a location/threshold shift parameterloc). For standard 2-parameter Weibull estimation, passingfloc=0is required to fix the threshold at zero.
import numpy as np
import scipy.stats as stats
from scipy.optimize import newton
# 1. Generate synthetic Weibull data (True k = 2.2, True lambda = 100)
np.random.seed(42)
true_k, true_lam = 2.2, 100.0
n = 200
data = stats.weibull_min.rvs(c=true_k, scale=true_lam, size=n)
# ---------------------------------------------------------
# Method A: Custom Maximum Likelihood (Newton-Raphson)
# ---------------------------------------------------------
def weibull_mle(t):
ln_t = np.log(t)
mean_ln_t = np.mean(ln_t)
def g(k):
t_k = t ** k
return np.sum(t_k * ln_t) / np.sum(t_k) - (1.0 / k) - mean_ln_t
k_hat = newton(g, x0=1.5)
lam_hat = (np.mean(t ** k_hat)) ** (1.0 / k_hat)
return k_hat, lam_hat
# ---------------------------------------------------------
# Method B: Weibull Probability Plot (Rank Regression / OLS)
# ---------------------------------------------------------
def weibull_rank_regression(t):
t_sorted = np.sort(t)
n_obs = len(t_sorted)
i = np.arange(1, n_obs + 1)
F_hat = (i - 0.3) / (n_obs + 0.4)
X = np.log(t_sorted)
Y = np.log(-np.log(1.0 - F_hat))
slope, intercept, r_value, _, _ = stats.linregress(X, Y)
k_hat = slope
lam_hat = np.exp(-intercept / k_hat)
return k_hat, lam_hat, r_value**2
# ---------------------------------------------------------
# Method C: SciPy Built-in MLE (Fixing location floc=0)
# ---------------------------------------------------------
k_scipy, loc_scipy, lam_scipy = stats.weibull_min.fit(data, floc=0)
k_mle, lam_mle = weibull_mle(data)
k_ols, lam_ols, r2 = weibull_rank_regression(data)
print(f"True Parameters -> k: {true_k:.4f}, lambda: {true_lam:.4f}")
print(f"Custom MLE -> k: {k_mle:.4f}, lambda: {lam_mle:.4f}")
print(f"SciPy MLE (floc=0) -> k: {k_scipy:.4f}, lambda: {lam_scipy:.4f}")
print(f"Rank Regression (OLS)-> k: {k_ols:.4f}, lambda: {lam_ols:.4f} (R² = {r2:.4f})")5. Key Modeling Considerations
- Diagnostic Linearized Plotting: Even when MLE is selected for final parameter estimation, plotting ln(−ln(1 − F̂(t))) against ln(t) reveals structural anomalies such as dog-leg bends or curvature—which indicate competing failure modes (a mixture of Weibull distributions) or the presence of a non-zero failure-free threshold requiring a 3-parameter Weibull model (t0 > 0).
- Small-Sample Bias in MLE: Maximum Likelihood estimators for k̂ exhibit an upward bias in small samples (n < 25). For small sample sizes, applying Ross’s unbiasing correction factor k̂unbiased ≈ k̂MLE ⋅ (n − 2) / (n − 0.68) or using Median Rank Regression (regressing X on Y) mitigates this distortion.
- Fixing the Location Parameter: Allowing a 3-parameter optimization routine to float the location parameter unrestricted on 2-parameter Weibull data frequently leads to boundary convergence issues and distorted shape estimates.
No comments:
Post a Comment