# Numerically computing the MLEs using Newton's method and the invariance proprty

I'm currently addressing maximum likelihood estimation for a three-parameter probability distribution with all parameters being positive real-valued. I'm using Newton's method to calculate the MLEs using R. Sometimes, in this case of small sample sizes, I may obtain: (1) negative estimators, (2) ill-conditioned Hessian matrix, or (3) execution error. I have implemented other "constrained optimization" algorithms in R, and I was able to get the right solution in all most all cases. However, the problem that I'm currently investigating is forcing me to focus only on Newton's method. Therefore, I thought about using the invariance property of the MLEs. To illustrate the problem that I have and the solution that I'm thinking about, consider the following example. Using R, I simulated

x <- c(2.325710e-07, 6.327766e-04, 1.877846e-03, 1.357671e-30, 2.550096e-159, 1.141821e-29, 1.180393e-37, 7.076208e-52, 4.134248e-20, 1.410155e-17)


from the gamma distribution with shape parameter $\alpha = 0.01$ and rate parameter $\beta = 1.0$. I'm interested in finding the MLEs using any optimization method, such as the optim function, for example. I simply define the log-likelihood function as follows:

gloglik <- function(par, dat){
-sum(dgamma(dat, shape = par[1], scale = par[2], log = TRUE))
}


Afterwards, I implement:

> optim(c(0.2, 0.9), gloglik, dat = x)$par > warnings()  to get: > optim(c(0.2, 0.9), gloglik, dat = x)$par
[1] 0.01297942 0.01947207
There were 13 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
2: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
3: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
4: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
5: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
6: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
7: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
8: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
9: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
10: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
11: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
12: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced
13: In dgamma(dat, shape = par[1], scale = par[2], log = TRUE) : NaNs produced


Clearly I have some warnings due to negative values (which would be errors in my case). Nevertheless, I can use the invariance property of the MLEs to bypass these warnings. In fact, by assuming that $\alpha = \exp(a)$ and $\beta = \exp(b)$, with $a = \log(0.01)$ and $b = 0.0$, and accordingly defining the log-likelihood function as follows:

gloglik2 <- function(par, dat){
-sum(dgamma(dat, shape = exp(par[1]), scale = exp(par[2]), log = TRUE))
}


Then by implementing:

exp(optim(c(-2.0, 0.01), gloglik2, dat = x)\$par)


I get:

[1] 0.01298458 0.01948171


which is a solution similar to the previous one but without any warnings.

My questions here, however, are:

1. Is what I did both computationally and statistically acceptable?
2. What is the impact on the efficiency of estimation (bias, variance, etc.)?
3. Are there any references (paper, textbook, etc.) that address the above issue?

Any advice would be appreciated. Thanks!