# Standard error with maximum likelihood estimation

ASTM E2238 describes the prediction of the largest inclusion in a polished steel surface of 150,000 mm² based on a Gumbel distribution of the largest inclusions measured in each of 24 polished microsections (150 mm² each). The location and scale parameters $$\lambda$$ and $$\delta$$ of the distribution are determined using the maximum likelihood method. Then, section 6.9.5 of the standard states:

The standard error, SE, for any inclusion of length $$x$$ based upon the ML method is $$SE(x)=\delta_{ML}\cdot \sqrt{(1.109 + 0.514\cdot y+0.608\cdot y^{2})/n}$$

where $$y = (x-\lambda )/ \delta$$ and $$n$$ is the number of measurements (24 in this standard).

(I can't upload the complete ASTM standard for copyright reasons, but there are a number of openly accessible sources that use this method, e.g. here, p. 22 here or eq. (6) here. They all use, without further comment, this formula).

My question is, where do these coefficients $$1.109$$, $$0.514$$ and $$0.608$$ come from?

It seems to me that the answer should be obvious to the professional as everyone uses that formula and no-one wonders about its orgin, so I'd like apologize in advance if this is a dumb question but I'm not a statistician. A reference to any textbook or other source would be highly welcomed.

Additional information as requested in comment: The SE formula applies to all values (including the largest value determined by the procedure) and is used to draw the confidence bands $$95\%CI=\pm 2\cdot SE(x)$$ as in fig 1. Yes, the cdf used is $$F(x) = e^{-e^{-(x-\lambda)/\delta}}$$. The largest inclusion is calculated for a return period of $$T = 1000$$ (i.e. $$A_{ref}=1000 A_{0}$$ where $$A_{0}=150 mm^{2}$$) by $$-\delta_{ML}\cdot ln(-ln((T-1)/T))+\lambda_{ML}$$.

• Can you provide some more details? Does this formula only apply for this specific situation where the entire surface is 1000 times larger in area than the area of the microsections? Does the reference define the Gumbel distribution as $F(y)=e^{-e^{-(y-\lambda)/\delta}}$? Is this supposed to be the standard error of the estimated largest inclusion in the entire surface? If so, does it provide a formula for the estimate of the largest inclusion in the entire surface? Commented Apr 19, 2021 at 14:52
• I see 24 points in Fig. 1. The x-axis is the length of the largest inclusion in each of those 24 microsections. Can you explain what the reduced $y_i$ means on the y-axis- is it $y=(x-\lambda)/\delta$? Then, the dashed curves are drawn by adding and subtracting 1.96 times the SE(x) formula in your question? Commented Apr 19, 2021 at 15:44
• @JohnL the ordinate is $-ln(-ln(F(y))) = -ln(-ln(P)))$ where P is the cumulative probabilty plotting position $P_{i}=i/(n+1)$ and yes, the dashed line is the solid line $\pm$ 2 SE (I guess they just rounded 1.96 to 2 in the standard).
– Stef
Commented Apr 19, 2021 at 16:03
• Are the coefficients just entries in the covariance matrix obtained from the Hessian of the maximum likelihood estimate? Is this "SE(x)" (should be SE(y)) the SE of a confidence interval for a predicted value obtained using the $\delta$-method? Commented Apr 19, 2021 at 18:05
• @Adam0 It is SE(x). x is a function of y. it's inverted from the way you would expect. the confidence limits are symmetric vertically, not horizontally in the figure. Commented Apr 19, 2021 at 18:45

Let $$\hat{\delta}_{ML}$$ and $$\hat{\lambda}_{ML}$$ be the maximum likelihood estimates of $$\delta$$ and $$\lambda$$.
Also, let $$\gamma$$ denote the Euler-Mascheroni constant (approximately 0.5772).

The asymptotic covariance matrix of $$\left(\hat{\lambda}_{ML},\hat{\delta}_{ML} \right)'$$ is:

$$\frac{\delta^{2}}{n \pi^2}\begin{bmatrix}6(\gamma-1)^2+\pi^2 & 6(1-\gamma)\\6(1-\gamma) & 6\end{bmatrix}$$

That is found from the inverse of the Fisher Information matrix. The Fisher information for a single observation $$X$$ is found by finding the expected value of the second derivatives $$\frac{d^2}{d\lambda^2},\frac{d^2}{d\lambda d\delta},\frac{d^2}{d\delta^2}$$ of $$\log f(x)$$; then finding the expected value of those where you change $$x$$ to a random variable $$X$$ having that Gumbel distribution. Then, take the negative of that matrix of expected values of second derivatives. Those are not trivial to find, but I did it in Mathematica. The Fisher information adds from independent samples, so you multiply by $$n$$ to get the Fisher information for the entire sample.

fx = D[E^(-E^(-(x - lam)/del)), x];
InputForm[
Inverse[-n Assuming[Element[lam, Reals] && del > 0,
Integrate[
{{D[Log[fx], lam, lam], D[Log[fx], lam, del]},
{D[Log[fx], lam, del], D[Log[fx], del, del]}} fx,
{x, -Infinity, Infinity}]]]]


returns

{{(del^2*(6*(-1 + EulerGamma)^2 + Pi^2))/(n*Pi^2),
(6*del^2*(1 - EulerGamma))/(n*Pi^2)},
{(6*del^2*(1 - EulerGamma))/(n*Pi^2), (6*del^2)/(n*Pi^2)}}


Now, to find the approximate variance of $$x=\hat{\delta}y+\hat{\lambda}$$, you use $$(1,y)\frac{\hat{\delta}_{ML}^{2}}{n \pi^2}\begin{bmatrix}6(\gamma-1)^2+\pi^2 & 6(1-\gamma)\\6(1-\gamma) & 6\end{bmatrix}(1,y)'$$

This is equal to $$\frac{\hat{\delta}_{ML}^{2}}{n}\left(\frac{6-12\gamma+6\gamma^2+\pi^2}{\pi^2}+\frac{12-12\gamma}{\pi^2}y+\frac{6}{\pi^2}y^2 \right)\approx\frac{\hat{\delta}_{ML}^{2}}{n}(1.10866+0.514044y+0.607927y^2)$$
The estimated standard error is the square root of the estimated variance:
$$\hat{\delta}_{ML}\sqrt{\frac{1.10866+0.514044y+0.607927y^2}{n}}$$

• At what point did you compute the inverse of the matrix? (This might explain the differences in signs.) It would help to explain why the factor $\delta^2/n$ is missing from your last equation.
– whuber
Commented Apr 19, 2021 at 18:10
• @whuber corrected now and I added Mathematica code used to find the Cramer-Rao lower bound Commented Apr 19, 2021 at 18:46
• @Stef it's fixed now. Mathematica uses a different definition of the Gumbel distribution. It uses $1-e^{-e^{(x-\lambda)/\delta}}$ I re-did the calculations with the definition $e^{-e^{-(x-\lambda)/\delta}}$ and the formula is correct. Commented Apr 20, 2021 at 14:46