I'm trying implement the Maximum Likelihood Estimation in R to Gumbel distribution, but the algorithm doesn't converge.
I'm using this parametrization to gumbel:
$${\frac {1}{\sigma}{{\rm e}^{{\frac {-(x-\mu)}{\sigma}}-{{\rm e}^{{\frac {-(x-\mu)}{\sigma}}}}}}}$$
Then the log-likelihood is:
$$-n\ln \left( \sigma \right) +\sum _{i=1}^{n}\ln \left( {{\rm e}^{Z_{ {i}}-{{\rm e}^{Z_{{i}}}}}} \right)$$ where $Z_i = \dfrac{-(x_i + \mu)}{\sigma}$
The Newton-Raphson technique to find the $\hat{\mu}$ and $\hat{\sigma}$ is:
$$\hat{\theta} = \theta_0 - H^{-1}(\theta_0)\,U(\theta_0)$$
where $\theta = (\mu, \sigma)^{T}$, $H$ is the hessian matrix and $U$ is the score vector.
For Gumbel distribution we have to:
\begin{equation} U(\mu,\sigma\mid x) = \begin{cases} {\frac {n}{\sigma}}+\sum _{i=1}^{n}-{\frac {{{\rm e}^{Z_{{i}}}}}{ \sigma}}\\ {\frac {\sum _{i=1}^{n} \left( -x_{{i}}+\mu \right) {{\rm e}^{Z_{{i}}} }+x_{{i}}+ \left( -\mu-\sigma \right) n}{{\sigma}^{2}}} \end{cases} \end{equation}
and $$ H(\mu, \sigma \mid x) = \left[ \begin {array}{cc} -{\frac {\sum _{i=1}^{n}{{\rm e}^{Z_{{i}}}} }{{\sigma}^{2}}}&{\frac {-n\sigma+\sum _{i=1}^{n}{{\rm e}^{Z_{{i}}}} \left( \sigma-x_{{i}}+\mu \right) }{{\sigma}^{3}}} \\ {\it\textrm{Sym}}&{\frac {2\,n\mu\,\sigma+n{\sigma}^{2}- \sum _{i=1}^{n} \left( -x_{{i}}+\mu \right) \left( \mu+2\,\sigma-x_{{i}} \right) {{\rm e}^{Z_{{i}}}}+2\,x_{{i}}\sigma}{{\sigma}^{4}}}\end {array} \right]$$
Simulated random variate from Gumbel I use the Inverse-Transform Method, that is:
$$\mu -\sigma\ln \left( -\ln \left( U(0,1) \right) \right)$$
My R code is:
log-Likelihood
logLH <- function(theta){
mu <- theta[1]
sigma <- theta[2]
Z <- -( (x - mu)/sigma)
t1 <- log(sigma)
t3 <- sum(log(exp(Z - exp(Z))))
l <- -t1 * n + t3
return(l)
}
Score Vector
U <- function(n, x, theta){
mu <- theta[1]
sigma <- theta[2]
Z <- -((x - mu)/sigma)
Es <- matrix(nrow=2)
Es[1] <- n / sigma - sum(1 / sigma * exp(Z))
Es[2] <- (sum((mu-x) * exp(Z) + x) - (mu + sigma) * n) / sigma ^ 2
return(Es)
}
Hessian Matrix
H <- function(n, x, theta){
mu <- theta[1]
sigma <- theta[2]
Z <- -((x - mu)/sigma)
Hs <- matrix(ncol=2,nrow=2)
Hs[1,1] <- -1 / sigma ^ 2 * sum(exp(Z))
Hs[1,2] <- (-n * sigma + sum(exp(Z) * (sigma - x + mu))) / sigma ^ 3
Hs[2,1] <- Hs[1,2]
Hs[2,2] <- (2 * n * mu * sigma + n * sigma ^ 2 - sum((-x + mu) *
(mu + 2 * sigma - x) * exp(Z) + 2 * x * sigma)) / sigma^4
return(Hs)
}
Random Variable Generation
n <- 100
mu <- .5
sigma <- 2
x <- mu - sigma*log(-log(runif(n)))
Newton-Raphson
tol <- 1e-6
iter <- 1
t.hat <- matrix(c(.3, 4))
while(max(abs(U(n,x,t.hat))) > tol & iter <= 20){
t.hat <- t.hat - solve(H(n,x,t.hat))%*%U(n,x,t.hat)
cat(iter, t.hat, "\n")
iter <- iter + 1
}
Using the library maxLik the method converge
library(maxLik)
maxLik(logLik = logLH, start = c(.3,4))
