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

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

• Questions only about R is prone to be closed, consider reformulate to have emphasis on the statistical question about estimation. Jul 12, 2020 at 15:27

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

• How should i do that? @Ben - Reinstate Monica Jul 12, 2020 at 21:52
• The value of log likelihood should be in negative(preferably -701.2132) , and the values of MLE of alpha and lambda is 0.2807 and 0.0030 according to the same paper i'm working on. Jul 13, 2020 at 4:05
• The value I have obtained appears to give a higher log-likelihood, so please check your numbers.
– Ben
Jul 13, 2020 at 4:27
• I'm not sure why the values are different when we are using the same data!! Jul 13, 2020 at 4:55
• @matrikashukla, please avoid asking too many or new questions in comments. It is OK to ask the OP to get greater clarity about an aspect that is unclear to you, but there is a limit to this. If you have a new or follow-up question, start a new thread. You can link back to this one for context, if appropriate. Jul 20, 2020 at 1:06