Distribution Fitting

Distribution Fitting

Maximum Likelihood, Goodness of Fit, Tolerance Limits and Basis Values

What the tool does

A sample of measurements of one quantity (strengths, lives, dimensions, times) is fitted by several candidate probability distributions, each fit is tested against the data, and the fits are used for the quantities a design or a specification needs: percentiles, tolerance limits that cover a stated fraction of the population with a stated confidence, and the A-basis and B-basis allowables of MIL-HDBK-17 / CMH-17. The sample is treated as one unstructured population: the handbook's screens for outliers (maximum normed residual, Section 8.3.3) and for batch-to-batch differences (k-sample Anderson–Darling, Section 8.3.2) must be applied before pooled data are used for an allowable.

The candidate distributions

• Normal, parameters μ, σ; fitted by x̄ and s (n − 1 denominator, as the handbook does). The log-likelihood is evaluated at the maximum-likelihood σ (n denominator) so that AIC compares like with like.
• Lognormal: ln x is normal with parameters μ_L, σ_L; everything is the normal case on the logarithms, transformed back by exp.
• Two-parameter Weibull, F(x) = 1 − exp[−(x/α)^β], shape β and scale α; the handbook's model of choice for the strength of brittle fibres and composites (chain-of-bundles argument). The maximum-likelihood shape solves G(β) = Σ x^β ln x / Σ x^β − 1/β − mean(ln x) = 0, a monotone equation solved by bisection (Section 8.3.4.2.1), and α = (Σ x^β/n)^1/β. Standard errors come from the observed information matrix.
• Three-parameter Weibull (optional), x − γ Weibull with a threshold γ below the smallest value. The threshold is found by maximising the profile likelihood; the likelihood is unbounded as γ approaches x₍₁₎ when the local shape falls below one, so the search stops short of the minimum and the page flags a fit that ran to that edge. Its extra parameter is charged by AIC.
• Exponential, F(x) = 1 − e^−λx, λ = 1/x̄: the constant-hazard model, usually a poor fit to strength data but the reference for life data.
• Gamma, shape k and scale θ: the shape solves ln k − ψ(k) = ln x̄ − mean(ln x) by Newton's method on the digamma function, θ = x̄/k, and the standard errors come from the Fisher information n[ψ′(k), 1/θ; 1/θ, k/θ²].

Models with positive support are skipped when the sample contains a value ≤ 0.

Comparing the fits

AIC = 2k − 2 ln L penalises each parameter by two units of log-likelihood; AICc adds the small-sample correction 2k(k + 1)/(n − k − 1) and BIC uses k ln n. The model with the smallest AIC ranks first; a difference under about 2 between two models is not decisive, and above 10 the poorer model has essentially no support. AIC says which model describes the sample best; it does not say whether any of them describes it adequately, which is the job of the goodness-of-fit tests.

Goodness-of-fit tests

Anderson–Darling. A² = −n − (1/n) Σ (2i − 1)[ln F(x_(i)) + ln(1 − F(x_(n+1−i)))] weights the tails of the distribution more than the Kolmogorov–Smirnov statistic and is the test the handbook prescribes. Because the parameters are estimated from the same sample, the critical values of the fully specified case do not apply; D'Agostino and Stephens (1986) give, for each distribution, a small-sample modification of A² and the critical values of the estimated-parameter case: normal and lognormal A* = A²(1 + 0.75/n + 2.25/n²) with the 15 / 10 / 5 / 2.5 / 1 % points 0.561, 0.631, 0.752, 0.873, 1.035 and the p-value formulas of their Table 4.9; Weibull (through the extreme-value distribution of ln x) A²(1 + 0.2/√n) with 0.637, 0.757, 0.877, 1.038 at 10 / 5 / 2.5 / 1 %; exponential A²(1 + 0.6/n) with 0.922, 1.078, 1.341, 1.606, 1.957; gamma by the fitted shape, Table 4.21. Between tabulated points the p-value is interpolated on a log scale and beyond them it is reported as a bound. The handbook (MIL-HDBK-17-1F 8.3.4.2.2 and 8.3.4.3.2) uses its own modifications, AD* = A²(1 + 0.2/√n) for the Weibull and A²(1 + 4/n − 25/n²) for the normal and lognormal, with the logistic fits OSL = 1/{1 + exp[−0.10 + 1.24 ln AD* + 4.48 AD*]} and 1/{1 + exp[−0.48 + 0.78 ln AD* + 4.58 AD*]} for the observed significance level, and accepts a model when OSL exceeds 0.05. Both are shown; they agree closely in the region that matters.

Kolmogorov–Smirnov. D is the largest vertical distance between the empirical and the fitted cumulative distribution. Its classical p-value (Kolmogorov's limit with Stephens' finite-n correction) assumes the distribution was specified in advance; with parameters estimated from the sample the true p-value is smaller (Lilliefors 1967), so the classical value is shown as a lower bound on the evidence and, for the normal and lognormal, the Lilliefors p-value is computed by the Dallal–Wilkinson (1986) approximation, the same as in R's nortest package.

Chi-square. The range is cut into k equiprobable classes under the fitted model (k = 2n^0.4, between 4 and 20, reduced so that at least five values are expected per class), χ² = Σ (O − E)²/E with k − 1 − (parameters) degrees of freedom. It is the least powerful of the three and depends on the class count; it is reported because many specifications still ask for it.

Shapiro–Wilk. W = (Σ a_i x_(i))² / Σ (x − x̄)² tests normality alone; the coefficients and the p-value follow Royston's AS R94 (1995), valid for 3 ≤ n ≤ 5000, the algorithm behind R and scipy. It is generally the most powerful omnibus test of normality for small samples.

Percentiles

x_p is the value below which a fraction p of the fitted population lies: the Weibull quantile α(−ln(1 − p))^1/β, the normal μ + z_pσ, the lognormal exp(μ_L + z_pσ_L), the exponential −ln(1 − p)/λ, and the gamma by inverting its cdf. The empirical quantile (R's type 7) is shown beside them. A fitted 0.1 % or 99.9 % percentile is an extrapolation whose reliability depends entirely on the model's tail being right, which is why the handbook insists on its order of preference rather than on the best-fitting model.

Tolerance limits

A tolerance limit is a bound that covers a stated fraction p of the population with a stated confidence γ, allowing for the sampling error of x̄ and s. For the normal model the one-sided limit is x̄ − k₁s with k₁ = t′_γ(n − 1; δ = z_p√n)/√n, the γ-quantile of the noncentral t distribution: exact, and the definition behind the handbook's k_B (p = 0.90) and k_A (p = 0.99) tables. The two-sided interval x̄ ± k₂s has no closed form; the tool solves the exact condition P[Φ(X̄ + k₂S) − Φ(X̄ − k₂S) ≥ p] = γ by integrating over the distribution of the sample mean (Wald and Wolfowitz 1946; Odeh and Owen 1980), and prints two approximations beside it: Howe's (1969) k₂ = z_(1+p)/2√[ν(1 + 1/n)/χ²_1−γ,ν], the form in the NIST e-Handbook, and the Wald–Wolfowitz approximation that produced the classical tables (NBS Handbook 91). For n = 10, p = 0.99, γ = 0.95 they give 4.437, 4.445 and 4.433. The lognormal limits are the same factors applied to ln x.

Without a distribution, order statistics give distribution-free limits. The rth smallest value is a lower limit covering a fraction p with confidence P[Bin(n, 1 − p) ≥ r]: the tool reports the largest r that reaches the chosen confidence, or, when even r = 1 does not, the smallest sample size that would (29 for 95 / 90, 299 for 95 / 99, the handbook's thresholds). The pair (x₍₁₎, x_(n)) covers a fraction p with confidence 1 − n pⁿ⁻¹ + (n − 1) pⁿ.

A-basis and B-basis values

An A-basis value is a 95 % lower confidence bound on the first percentile of the population, a B-basis value the same bound on the tenth percentile (MIL-HDBK-17-1F 8.1.4; the structural "A" and "B" allowables of MIL-HDBK-5 and MMPDS have the same definition). Section 8.3.4 of the handbook computes them for unstructured data in a fixed order:

1. Fit the two-parameter Weibull and test it; if OSL > 0.05 use it, even when another model fits better. B = Q exp[−V_B/(β√n)] with Q = α(0.10536)^1/β, the fitted tenth percentile, and V_B from Table 8.5.8 (the 95 % point of the conditional distribution of √n[ln Q̂ − ln Q]β̂, Lawless 1982); for A, Q = α(0.01005)^1/β and V_A from Table 8.5.9. The tables start at n = 10 and are interpolated linearly between the tabulated sample sizes (V is convex, so the interpolation is conservative); the handbook's approximation formulas 8.5.8(h) and 8.5.9 are used beyond their last entries. Two misprinted cells of the 1F tables (n = 59 and 425 of Table 8.5.8, n = 196 of Table 8.5.9) are replaced by the mean of their neighbours.
2. Otherwise test the normal: B = x̄ − k_Bs, A = x̄ − k_As with the one-sided factors above, computed exactly from the noncentral t (the handbook's Tables 8.5.10 / 8.5.11 used the rounded noncentralities 1.282√n and 2.326√n, which puts their entries up to 0.001 above the exact values).
3. Otherwise test the lognormal: the normal procedure on ln x, exponentiated.
4. When all three are rejected, the nonparametric method. For n ≥ 29 (B) or n ≥ 299 (A) the value is the rth smallest observation with r from Tables 8.5.12 / 8.5.13, which the tool reproduces from their defining criterion, the largest r with P[Bin(n, 1 − p) ≥ r] ≥ 0.95. For smaller samples the Hanson–Koopmans (1964) bound B = x_(r)(x₍₁₎/x_(r))^k with r and k from Table 8.5.14 (n = 2 … 28), and A = x_(n)(x₍₁₎/x_(n))^k_A with k_A from Table 8.5.15 (n = 2 … 299; between tabulated sizes the next lower entry is used, which is conservative). The method assumes a population whose log-cdf is concave, which strength data satisfy in practice.

The page prints every method's A and B values for comparison and marks the one the rule selects. The handbook reports a basis value rounded to the precision of the data (104.41 becomes 104) and only after the batch and outlier screens.

Worked example on the page

The default sample is Problem 1 of MIL-HDBK-17-1F Section 8.3.7: thirty compressive strengths from ten batches (the recorded 1444.5 corrected to 144.4 as the handbook's outlier step does). The handbook prints β̂ = 15.35, α̂ = 128.39, AD = 0.699, AD* = 0.7245, OSL = 0.0576, V_B = 5.057, Q = 110.88 and B = 104.41; the page gives β̂ = 15.37 (its bisection is tighter than the handbook's), α̂ = 128.39, OSL = 0.061 and B = 104.44, and the test script checks the library against these values, the handbook's normal, lognormal, nonparametric and Hanson–Koopmans problems, its k and V tables, the NIST e-Handbook tolerance-factor examples, the classical two-sided tables, the Shapiro–Wilk (1965) example and an independent scipy run.

Limits of the tool

• Complete samples only: censored or interval data (suspensions) belong in the Weibull Life Analysis tool.
• No batch structure: pooled data must first pass the handbook's k-sample Anderson–Darling test; structured data need the ANOVA or regression basis values of Section 8.3.5.
• The outlier screen (maximum normed residual) is not applied; inspect the probability plots.
• The Weibull V tables of the 1F handbook are used as printed (apart from the three misprints); CMH-17-1G revised some entries by a few thousandths.
• The three-parameter Weibull standard errors are from a numerical Hessian and are omitted when the threshold sits at the edge of its search range.