I'm doing a project on Extreme Value Theory in R and I'm having to code my own maximum likelihood estimates for the parameters for each of the 3 distributions.
The mle's for extreme value distributions don't have analytic solutions, they must be approximated numerically. The Gumbel case was easy, as I could get one parameter by itself and use the result to optimize the other parameter numerically.
For the Fréchet and Weibull cases, you have to optimize more than one parameter simultaneously. I'm trying to use the optim() function in R but getting very weird results. I'm not sure if this is a problem with my code, a problem with my understanding of optim(), or with the theory itself (or some combination of the 3).
Anyway, here's my code:
library(evd) #using this for frechet random variables
#generate maxima from 2 distributions: 1st should converge to Frechet, 2nd to Weibull
Mfrechet<-function (m) {
d1<-matrix(rfrechet(1000000),nrow=m)
apply(d1,1,max)/(1000000/m)
}
Munif <- function (m) {
d1<-matrix(runif(10000000),nrow=m)
n<-ncol(d1)
Xn<-apply(d1,1,max)
n*Xn -n +1
}
x1<-Mfrechet(10000)
x2<-Munif(10000)
#negative log-likelihood for frechet
lf.fr <- function(y) {
x<-x1
n<-length(x1)
a<-y[1]
b<-y[2]
g<-y[3]
-(n*log(g)+n*g*log(b)-(g+1)*sum(log(x-a),na.rm=T)-sum((b/log(x-a))^g,na.rm=T))
}
#negative log likelihood function for weibull
lf.weibull <- function(y) {
x<-x2
n<-length(x2)
a<-y[1]
b<-y[2]
g<-y[3]
-(n*log(g)+n*g*log(b)+(g-1)*sum(log(a-x),na.rm=T)-sum(((a-x)/b)^g,na.rm=T))
}
optim(c(1,1,1),lf.fr)
optim(c(1,1,1),lf.weibull)
The results are way off, especially compared with fevd from extRemes package:
library(extRemes)
fevd(x1)
fevd(x2)
I tried changing the starting values for optim, but it gives me wildly different results but not close to correct ones. I've used getAnywhere to look at most gev fitting in different packages, they mostly use optim as well, but they are pretty difficult to sort through.
Can anyone show me where I am going wrong? And hopefully a way to fix it that isn't too complicated...
Your problem is two-fold. First, the GEV distribution has support which depends on the parameters. So if you know the parameters, it is a simple matter to calculate the support, but if you are searching for the parameters in an optimiizer is is easy to land on a combination which is inconsistent with your data. Your loglik functions need to deal with that.
Second, whether GEV corresponds to Weibull or Frechet depends on the shape parameter, xi (xi > 0 implies Frechet; xi < 0 implies Weibull). So you need to use an optimizer which allows you set lower and upper bounds on the parameters. If you use optim(...) your only choice is "L-BFGS-B". Unfortunately, this method is somewhat fragile (see this post). An alternative, identified in that same post, is nmkb(...) from the dfoptim package.
Here is a working example, based on the parameterization of the GEV distribution defined in Wikipedia.
library(evd)
library(dfoptim)
#
ll = \(par, x) {
mu <- par[1]; sigma <- par[2]; xi <- par[3]
s <- (x - mu)/sigma
t <- (1 + xi*s)^(-1/xi)
f <- t^(xi+1) * exp(-t) / sigma
if(any(is.na(f)) | any(f <= 0)) return(1e6)
return (-sum(log(f)))
}
par <- c(mu=0, sigma=.5, xi=-0.1)
result <- nmkb(par, ll, lower = c(-Inf, 0, -Inf), upper = c(Inf, Inf, 0), x=x2)
result$par
## [1] -0.001692221 0.998771474 -0.997121661
fevd(x2)$results$par
## location scale shape
## -0.001711532 0.998739426 -0.997070437
#
par <- c(mu=0, sigma=0.5, xi=+0.1)
result <- nmkb(par, ll, lower = c(-Inf, 0, 0), upper = c(Inf, Inf, Inf), x=x1)
result$par
## [1] 0.9997442 1.0003228 1.0004027
fevd(x1)$results$r
## location scale shape
## 0.9999759 1.0005464 1.0003256
Note that you did not set the seed for the random number generator in your post, so your x1 and x2 will be different. I used set.seed(1).