# Can we fit extreme value distribution by build-in package?

I try to find a package in R to fit Gumbel distribution by Block Maxima Approach using maximal likelihood function (see here) $$G(x; \mu , \sigma)=\exp[-e^{-\frac{x-\mu}{\sigma}}].$$

The block maxima approach is that we have $$n$$ iid random samples $$X_1, X_2, \dots, X_n$$ from $$exp(1)$$ (the limiting distribution of extreme value is Gumbel), and choose $$m$$ blocks with same size for $$X_i$$ and choose the maximal value of each block denoted by $$Z_1, Z_2,\dots, Z_m$$. Since $$Z_i$$ is the maximal value in $$i-th$$ block, we can approximate its distribution by Gumbel distribution.

In this website https://www.stat.purdue.edu/~huang251/1018.html, he used the package fevd to fit the generalized extreme value distribution with three parameters $$(\mu, \sigma, \xi)$$. Can we just fit two parameters? Is there any such package to fit Gumbel distribution?

• Looking at the documentation for the fevd() function it seems like you can just add the argument type = "Gumbel". Is there a reason you haven't tried this, or that it won't work? Commented Nov 6, 2022 at 6:13
• @awhug I am not familiar with this package. Where do I add this argument? Commented Nov 7, 2022 at 17:45
• @awhug Do you mean fevd(data, type="Gumbel")? Commented Nov 7, 2022 at 17:46
• Yes. E.g. with simulated data x <- rexp(1000, 1.3); x_agg <- aggregate(x, list(gl(100, 10)), max); fit <- fevd(x_agg$x, type = "Gumbel") Commented Nov 8, 2022 at 1:24 • Are you sure?$-log(X)$where$X$is exponential does give you that Gumbel distribution, but here you're fitting a model to the maxima across blocks of$X$, not every observation. Commented Nov 8, 2022 at 6:52 ## 1 Answer Your simulation appears to be working as expected. Here's a verification using optim() to directly maximise the log-likelihood function you've written yourself in the question (noting that the number of blocks, $$m$$, should actually be 1000 and not 100). The results are equivalent to what fevd() is estimating (within a margin of error). # Data generating process n <- 100 # observations within each block m <- 1000 # number of blocks x <- rexp(n*m, 1) x_agg <- aggregate(x, list(gl(m, n)), max) # Optional: Rescale to force mu=0 and sigma=1 # x_agg$$x <- x_agg$$x - log(n) # Log-likelihood function: Eq. 3.9 Coles (2001) pg. 55 ll_fn <- function(par, z, m){ mu <- par[1] sigma <- par[2] -m * log(sigma) - sum((z - mu)/sigma) - sum(exp(-((z - mu)/sigma))) } # Fit using optim ---- # Note fnscale set to maximise ll_fn # hessian = TRUE for Fisher Info if needed fit1 <- optim(par = c(1, 1), fn = ll_fn, z = x_agg$$x, m = m, hessian = TRUE, control = list(fnscale = -1)) #$$par # [1] 4.661471 1.006802 # #$value
# [1] -1588.832

# Fit using extRemes ----
fit2 <- fevd(x_agg$x, type = "Gumbel") # Negative Log-Likelihood Value: 1588.832 # # Estimated parameters: # location scale # 4.661528 1.006741  In fact, without fitting a model at all you can verify that if the block maxima in this simulation are Gumbel distributed, this distribution can't be centred on $$\mu = 0$$. Here's what the raw x_agg$x looks like just plotting the data

# Simple histogram and density overlay
hist(x_agg$x, probability = TRUE, breaks = 38, xlab = "Block Maxima - rexp(1)") lines(density(x_agg$x), lty = 2)


The proof you've reproduced above leading you to expect the parameter estimates to approach the standard Gumbel ($$\mu = 0$$, $$\sigma = 1$$) is actually based on rescaling the block maxima by the log number of observations within each block, i.e., $$\text{max}\{X_1, \dots, X_m\} - \log(n)$$. This rescale does indeed result in the standard Gumbel, which you can verify using the code above, but note that when working with the raw block maxima, the parameter estimates will differ.

A more comprehensive discussion on this is provided in Chapter 3 of Coles, S. (2001). An introduction to statistical modeling of extreme values. London: Springer.

• I can't look at this again until tonight, but just quickly - I'd guess that $\mu$ here may be centred on the inverse CDF of the exponential at the quantile $1 - (1/n)$ where $n$ is the number of observations in each block (e.g. qexp(1 - (1/100)) above). Can't prove this though. It looks like you've lifted your proof above directly from Example 3.1 in Coles (2001, pg. 52), but it's difficult to follow how it relates to the rest of the question. Commented Nov 9, 2022 at 0:35
• I've just revised the answer and modified the code. You should be able to rescale the block maxima to get your standard Gumbel and check the MSE that way. Commented Nov 9, 2022 at 16:52
• See point #13 here for the same proof with cleaner notation. Commented Nov 9, 2022 at 16:56