You can find MLE equations for this distribution in [Mahdvai and Kundu (2017)](https://www.tandfonline.com/doi/abs/10.1080/03610926.2015.1130839) (accessible version [here](http://home.iitk.ac.in/~kundu/abbas-kundu-rev-1.pdf)).  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.  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.  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, we 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