How can I find MLE estimates for the parameters in the function below?

Question

I have a new function given which I would like to estimate its parameters; a,b,alpha,vartheta using MLE. How do I make an estimation?

EMHL<-function(a,b,alpha,vartheta) {
  (2*a*b*alpha*vartheta*
    (x^(vartheta-1))* exp(-x^vartheta) *
    ((1-exp(-x^vartheta))^(a-1)) * 
    ((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) * 
    (1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}

for a given dataset

x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 
      1.7, 2.7, 4.1, 1.8, 1.5, 1.2, 1.4, 3, 
      1.7, 2.3, 1.6, 2.0)

Do you have the values of your function observed for the given values of x? — qdread, Jan 18 '22 at 14:24
The distribution is identical to the one you posted in a previous question. https://stackoverflow.com/a/70688756/9463489 When `x > 1`, the `(1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))` term quickly gets large and dominates the other terms as `vartheta` gets large. What are the range of values each of the parameters can take? — jblood94, Jan 18 '22 at 14:39
Another data I can use for `x<1` is `x<- c(0.067, 0.068, 0.076, 0.081, 0.084, 0.085, 0.085, 0.086, 0.089, 0.098, 0.098, 0.114, 0.114, 0.115, 0.121, 0.125, 0.131, 0.149, 0.160, 0.485) ` — Bakang, Jan 18 '22 at 15:09

G. Grothendieck · Accepted Answer · 2022-01-19T12:39:54.700

There are syntax errors in the function shown in the question so we used the one in the Note at the end.

If that function is a density and you want to minimize the corresponding negative log likelihood then trying a few different starting values these seem to result in convergence.

(For the x given in a comment below the question list(a = 1, b = 1, alpha = 575, vartheta = 0.01) seems to work as starting values.)

NLL <- function(par) -sum(log(EMHL(par[1], par[2], par[3], par[4])))
st <- list(a = 1, b = 1, alpha = 525, vartheta = .2)
res <- optim(st, NLL); res

giving:

$par
           a            b        alpha     vartheta 
8.845296e-01 1.211526e+00 5.315759e+02 1.326975e-03 

$value
[1] -10904.36

$counts
function gradient 
     327       NA 

$convergence
[1] 0

$message
NULL

L-BFGS-B

Using L-BFGS-B with the function above results in problems but an answer can be obtained if we constrain the parameters. For example, this converges with all the constraints being active, i.e. the solution is on the boundary of the feasible region. Note that tests and p values derived from estimates on a boundary may not be valid.

optim(c(1, 1, 1, 1), NLL, method = "L-BFGS-B", lower = 0.01, upper = 2)

Other densities

Another possibility is to use a different distribution. The Cullen & Frey diagram produced by descdist suggests that the gamma distribution may be close

library(fitdistrplus)
descdist(x)

(continued after graph)

or we could try the generalized gamma (dgengamma) in the flexsurv package.

library(bbmle)
library(flexsurv)
NLLgeng <- function(mu, sigma, Q) -sum(dgengamma(x, mu, sigma, Q, log = TRUE))
m <- mle2(NLLgeng, list(mu = 1, sigma = 1, Q = 1), optimizer = "nlminb")
summary(m)

library(fitdistrplus)
fit <- fitdist(x, "gengamma", start = list(mu = 1, sigma = 1, Q = 1))
summary(fit)
plot(fit)

Note

EMHL<-function(a,b,alpha,vartheta) {
  2 * a * b * alpha * vartheta * 
  (x^(vartheta-1))* exp(-x^vartheta) *((1-exp(-x^vartheta))^(a-1)) * 
  ((1 - (1 -( 1 - exp(-x^vartheta))^a)^b)^(alpha-1)) * 
  (1- (1 - exp(-x^vartheta))^a )^(-b*(alpha+1))
}
x<- c(1.1, 1.4, 1.3, 1.7,1.9, 1.8, 1.6, 2.2, 1.7, 2.7, 4.1, 1.8,
   1.5, 1.2, 1.4, 3, 1.7, 2.3, 1.6, 2.0)

Thank you a lot. But when I use the above estimates as initial values in MLE below and try to get `-2 Log L` and `p-values` I get NaNs.... The code I used for `MLE estimates, P-values and -2Log L ` is `EHML.result<-mle2(EMHL,hessian = NULL,start=list(a=1.8845296,b=1.21152,alpha=531.5,vartheta=0.013),optimizer="nlminb",lower=0)` — Bakang, Jan 18 '22 at 16:32
The first argument to mle2 is the neg log lik, not the density. If EMHL and st are as in the answer then: `bbmle::mle2(function(a, b, alpha, vartheta) -sum(log(EMHL(a,b,alpha,vartheta))), st, method = "Nelder")` — G. Grothendieck, Jan 18 '22 at 17:51

How can I find MLE estimates for the parameters in the function below?

1 Answers1

L-BFGS-B

Other densities

Note