How do I find the MLE of the APE distribution in R?

Question

The random variable $Y$ is said to have a two-parameter APE distribution, denoted by $\text{APE}(\alpha, \lambda)$, with the shape parameter $\alpha>0$ and scale parameter $\lambda>0$ if the density function is:

$$f_Y(y) = \begin{cases} \log (\frac{\log \alpha}{\alpha-1}) \cdot \lambda e^{-\lambda y} \cdot \alpha^{1-e^{-\lambda y}} & & & \text{for } \alpha \neq 1 \\[8pt] \lambda e^{-\lambda y} & & & \text{for } \alpha = 1 \\[6pt] \end{cases}$$

Let $Y_1, Y_2,..., Y_n \sim \text{IID APE}(\alpha,\lambda)$ be a random sample from the APE distribution. Then the log-likelihood function is:

$$\ell_\mathbf{y}(\alpha,\lambda) = n \log \alpha + n \log \bigg( \frac{\log \alpha}{\alpha-1} \bigg) + n \log \lambda - \lambda \sum_i y_i - (\log \alpha) \sum_i e^{-\lambda y_i}.$$

How do I find MLE of this distribution in R? My data is given below:

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

Answer 1

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. 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