Maximum Likelihood Estimation of Gumbel Distribution
February 6, 2020 · 3 min read. Edited on May 14, 2020.
This write up assumes basic knowledge of the estimation theory.
Let’s start with the probability density function of the gumbel distribution, which is
$$p(x)=\frac{1}{\beta}e^{\frac{x-\alpha}{\beta}}e^{-e^{\frac{x-\alpha}{\beta}}}$$
$$\text{where }\alpha \text{ and } \beta \text{ are parameters and }\alpha \in \mathbb{R}, \beta > 0$$
Maximum liklihood solution can be written as,
$$\mathrm{f}_{M L}=\underset{\alpha \in \mathbb{R}, \beta>0}{\operatorname{argmax}} \in \mathbb{R}, \beta>0(P(D | \alpha, \beta))$$
where D represents the set of datapoints.
$$\mathrm{f}{M L}=\underset{\alpha \in \mathbb{R}, \beta>0}{\operatorname{argmax}} \in \mathbb{R}, \beta>0(\prod{i=1}^n P(x_i | \alpha, \beta))$$
where xi represents each datapoint.
Now onto calculating the log-likelihood,
$$=\ln \left(\prod_{i=0}^{n} P\left(x_{i} | \alpha, \beta\right)\right)$$
$$=\sum_{i=0}^{n} \ln P\left(\left(x_{i} | \alpha, \beta\right)\right)$$
$$=\sum_{i=0}^{n} \ln \left(\frac{1}{\beta} e^{-\frac{y_{i}-\alpha}{\beta}} e^{-e^{-\frac{j_{j}-a}{\beta}}}\right)$$
We end up with an equation with 2 variables in it, we’ll be using the Newton-Raphson method to approximate the roots of the equation. I recommend reading through the link if you’re familiar with the approximation method.
Finding derivatives, with respect to both α and β.
$$\frac{\partial f}{\partial \beta}=\sum_{i=0}^{n} \frac{x_{i}-n}{\beta^{2}}-\frac{n}{\beta}-\sum_{i=0}^{n} \frac{x_{i}-\alpha}{\beta} e^{-\frac{z_{i}-\alpha}{\beta}}$$
$$\frac{\partial f}{\partial \alpha}=\frac{n}{\beta}-\frac{1}{\beta} \sum_{i=0}^{n} e^{-\frac{z_{i}-\alpha}{\beta}}$$
For the Newron-Raphson approximation method, we require double derivatives, which we use to calculate the Hessian matrix.
$$\frac{\partial^{2} f}{\partial \beta^{2}}=\frac{n}{\beta^{2}}-\frac{2}{\beta^{2}} \sum_{i=0}^{n}\left(x_{i}-\alpha\right)+\frac{2}{\beta^{3}} \sum_{i=0}^{n}\left(x_{i}-\alpha\right) e^{\frac{-\left(x_{i}-\alpha\right)}{\beta}}+\frac{2}{\beta^{4}} \sum_{i=0}^{n}\left(x_{i}-\alpha\right)^{2} e^{\frac{-\left(x_{j}-\alpha\right)}{\beta}}$$
$$\frac{\partial^{2} f}{\partial \alpha^{2}}=\frac{-i}{\beta^{2}} \sum_{i=0}^{n} e^{-\frac{x_{i}-a}{\beta}}$$
$$\frac{\partial^{2} f}{\partial \alpha \beta}=-\frac{n}{\beta^{2}}+\frac{1}{\beta^{2}} \sum_{i=0}^{n} e^{\frac{-\left(x_{i}-\alpha\right)}{\beta}}-\frac{1}{\beta^{3}} \sum_{i=0}^{n}\left(x_{i}-\alpha\right) e^{\frac{-\left(x_{i}-a\right)}{\beta}}$$
The Hessian matrix can be calculated by,
$$H=\left[\begin{array}{ll}\frac{\partial^{2} f}{\partial \alpha^{2}} & \frac{\partial^{2} f}{\partial \alpha \beta} \ \frac{\partial^{2} f}{\partial \alpha \beta} & \frac{\partial^{2} f}{\partial \beta^{2}}\end{array}\right]$$
We also need a f matrix with the equations we’re solving for,
$$f=\left[\begin{array}{l}\frac{\partial f}{\partial \alpha} \ \frac{\partial f}{\partial \beta}\end{array}\right]$$
Now we can follow the following algorithm to estimate the parameters.
Step 1 : We’ll choose starting values for α and β, α(0) and β(0) using the method of moment estimators.
β = 0.7977∗ standard deviation of the dataset α = mean of the dataset - 0.4501 * standard deviation
Step 2 :
Choose a tolerance t for the change, 10^(−10) in our case. If the level of changes are less than the tolerance the iterations will break.
Step 3 :
Obtain inverse of the Hessian Matrix inv(H(α(0))(β(0))) and f(α(0))(β(0))
Step 4 :
Obtain new values of α and β, α(new) and β(new), from the Newton Raphson algorithm,
$$\left[\begin{array}{l}\alpha^{(n e w)} \ \beta^{(n e w)}\end{array}\right]=\left[\begin{array}{l}\alpha^{(o l d)} \ \beta^{(o l d)}\end{array}\right]-H^{-1}\left(\alpha^{(o l d)}, \beta^{(o l d)}\right) f\left(\alpha^{(o l d)}, \beta^{(o l d)}\right)$$
Step 5 :
Check to see if the differences between new and the old values are small enough and compare it by the set tolerance,
$$\left[\left(\alpha^{(n e w)}-\alpha^{(o l d)}\right)^{2}+\left(\beta^{(n e w)}-\beta^{(o l d)}\right)^{2}<t\right]$$
- if No, then return to Step(3), and calculate H−1(α(new))(β(new)) and f(α(new))(β(new)) and keep iterating.
- if Yes, then stop.
Following are the results of the above algorithm for a gumbel destibution generated using α = 2.3 and β = 4.0 ran multiple times and averaged.
