MATLAB: How to fit a Gumbel Maximum distribution using both Maximum Likelihood Estimates and Method of Moments Estimates

Question

How can I fit a Gumbel Maximum distribution using both Maximum Likelihood Estimates and Method of Moments Estimates?

Answer 1

The usual Gumbel distribution models the *minimum *of a sample and it is captured in MATLAB by the "Extreme Value Distribution":

https://www.mathworks.com/help/stats/extreme-value-distribution.html

The maximum likelihood estimates in this distribution for a particular dataset are obtained by means of the "evfit" function:

https://www.mathworks.com/help/stats/evfit.html

However, if we would like to model the *maximum* instead, we need to resort to a different approach. Essentially, the Gumbel maximum distribution is the mirror image of the Gumbel minimum distribution and, therefore, we can still model it using the "Extreme Value Distribution". Transforming the dataset by taking negatives allows us to send maximums to minimums (and vice versa). As a consequence, fitting a "Extreme Value Distribution" will now model the maximum:

 

rng default;
xMaxima = max(randn(1000,500), [], 2);
pd = fitdist(-xMaxima,'ev')

The maximum likelihood estimates are given by "evfit":

 

paramEstsMaxima = evfit(-xMaxima)

Bear in mind that any subsequent computations are now based on the transformed dataset. For example, if we would like to calculate a maximum extreme value at p=0.99, we would need to obtain a minimum extreme value at p=0.01 in the counterpart dataset and transform it back to maximum by taking negatives:

 

value = - icdf(pd,0.01)

Yet if we would like to use method of moments estimates, we would need to resort to another approach. In this case, the starting point is the "Generalized Extreme Value" distribution:

https://www.mathworks.com/help/stats/generalized-extreme-value-distribution.html

If the shape parameter "k" of this distribution is set to 0, then this distribution is suitable for modelling the maximum as well. The method of moments estimates for Gumbel maximum are given by:

 

gamma = 0.5772;
sigmaHat = sqrt(6)*std(xMaxima)/pi;
muHat = mean(xMaxima) - gamma*sigmaHat;

where "gamma" is the Euler constant. Instantiating the "Generalized Extreme Value" distribution with the above parameters produces a Gumbel Maximum with method of moments estimates:

 

pd = makedist('GeneralizedExtremeValue','k',0,'sigma',sigmaHat,'mu',muHat);

In this case, the maximum extreme value at p=0.99 would be given by:

 

extremeValue = icdf(pd,0.99);