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

Question

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!

Answer 1

It is very common practice to reparameterize a constrained parameter, such as the shape and scale of a gamma distribution, to a parameterization in which the domain is the real line. Again using the gamma distribution, this means using log shape/scale.

In 999 out 1000 times, this is a very good idea computationally; unconstrained methods are easier than constrained methods. An exception to that rule is when solution contains parameter values on the boundaries (i.e. estimated probabilities equal to 0 for example). In such cases, this often means that solutions in the unconstrained space have non-finite solutions. This is bad, computationally.

In terms of the statistical ramifications, it's not a bad idea either. As you point out, the MLE is invariant, so you can always transform the unconstrained solution back to the constrained solution, and it will be the correct answer. This does not affect the bias in any way; the invariance of the MLE tells us it's the exact same solution.

We do have to be slightly careful about standard errors; we have to remember that when using the expected Fisher's Information to find the covariance, this is reference the current parameterization. So if we used log shape and log scale in our gamma distribution, our standard errors will be in regards to log shape, not shape. Not really a problem; if we want something like a confidence interval for shape, we can first make a confidence interval for log shape, and then merely transform the confidence interval.

Closing remark; transforming the variables for simplifying the computations is not only not a problem, it's actually common practice.