You can find MLE equations for this distribution in Mahdvai and Kundu (2017) (accessible version here). As you can see from the paper, computing the MLE requires you to solve a critical point equation for $\lambda$ and you can then compute the MLE for $\alpha$ from this. This The paper also contains further information on the asymptotic distribution of the MLE, etc.
Implementation in R: This can be done in R using nonlinear optimisation with the nlm function, or by solving the critical point equation with the uniroot function. The paper also contains further information on the asymptotic distribution Using one of the critical point equations, Mahdvai and Kundu (2017) give the MLE of the first parameter as the function:
$$\hat{\alpha}(\mathbf{y},\lambda) = \exp \Bigg( \frac{\sum_i y_i - n/\lambda}{\sum_i y_i e^{-\lambda y_i}} \Bigg).$$
You can substitute this function into the log-likelihood function or the remaining critical point equation. In the code below, etcwe will substitute into the log-likelihod function and then maximise using the nlm function. (As the starting point for the iterative optimisation procedure, we will use the MLE for the exponential distribution.)
#Set the MLE function for alpha
LOG_ALPHA_HAT <- function(y, lambda) {
n <- length(y);
NUM <- sum(y) - n/lambda;
DEN <- sum(y*exp(-lambda*y));
NUM/DEN; }
#Set the log-likelihood function
LOGLIKE <- function(y, lambda) {
la <- LOG_ALPHA_HAT(y, lambda);
if (la == 0) {
LL <- n*log(lambda) - lambda*sum(y); } else {
LL <- n*la + n*log(la/expm1(la)) + n*log(lambda) -
lambda*sum(y) - la*sum(exp(-lambda*y)); }
LL; }
#Input the data
DATA <- c(1, 4, 4, 7, 11, 13, 15, 15, 17, 18, 19, 19, 20, 20, 22, 23, 28,
29, 31, 32, 36, 37, 47, 48, 49, 50, 54, 54, 55, 59, 59, 61, 61,
66, 72, 72, 75, 78, 78, 81, 93, 96, 99, 108, 113, 114, 120, 120,
120, 123, 124, 129, 131, 137, 145, 151, 156, 171, 176, 182, 188,
189, 195, 203, 208, 215, 217, 217, 217, 224, 228, 233, 255, 271,
275, 275, 275, 286, 291, 312, 312, 312, 315, 326, 326, 329, 330,
336, 338, 345, 348, 354, 361, 364, 369, 378, 390, 457, 467, 498,
517, 566, 644, 745, 871, 1312, 1357, 1613, 1630);
#Maximise the log-likelihood function
OBJECTIVE <- function(lambda) { - LOGLIKE(y = DATA, lambda) }
START <- c(1/mean(DATA))
NLM <- nlm(OBJECTIVE, p = START);
LLMAX <- NLM$minimum;
MLE_LAMBDA <- NLM$estimate;
MLE_ALPHA <- exp(LOG_ALPHA_HAT(y, MLE_LAMBDA));
MLE <- data.frame(alpha = MLE_ALPHA, lambda = MLE_LAMBDA, loglike = LLMAX);
rownames(MLE) <- 'MLE';
We can now display the MLE computed using this optimisation:
#Show the MLE
MLE;
alpha lambda loglike
MLE 0.00366583 0.0009550325 700.6492
