Coding likelihood and log-likelihood function in r to perform optimization

I am working on a paper that requires me to find the MLE of Gumbel’s type I bivariate exponential distribution. I have proved the likelihood and log-likelihood functions likelihood and log-likelihood but I am struggling to implement it in r to perform optimization with Optim function. My code generates NA values. Below are my codes.

# likelihood function of x
likelihood.x = function(params, data) {
  lambda1 = params[1]
  lambda2 = params[2]
  theta = params[3]
  A = (1 - theta) * (lambda1 * lambda2)
  B = theta * (lambda1 ^ 2) * lambda2 * data$X1
  C = theta * lambda1 * (lambda2 ^ 2) * data$X2
  D = (theta ^ 2) * (lambda1 ^ 2) * (lambda2 ^ 2) * data$X1 * data$X2
  E = (lambda1 * data$X1) + (lambda2 * data$X2) + (theta * lambda1 * lambda2 * data$X1 * data$X2)
  f = sum(log(A + B + C + D)) - sum(E)
  return(exp(f))
}

# Log-likelihood function of x
log.likelihood.x = function(params, data){
  lambda1 = params[1]
  lambda2 = params[2]
  theta = params[3]
  A = (1 - theta) * (lambda1 * lambda2)
  B = theta * (lambda1 ^ 2) * lambda2 * data$X1
  C = theta * lambda1 * (lambda2 ^ 2) * data$X2
  D = (theta ^ 2) * (lambda1 ^ 2) * (lambda2 ^ 2) * data$X1 * data$X2
  E = (lambda1 * data$X1) + (lambda2 * data$X2) + (theta * lambda1 * lambda2 * data$X1 * data$X2)
  f = sum(log(A + B + C + D)) - sum(E)
  return(-f)
}

Here's the function for generating the data

# Simulating data
rGBVE = function(n, lambda1, lambda2, theta) {
  x1 = rexp(n, lambda1)
  lambda12 = lambda1 * lambda2
  pprod = lambda12 * theta
  C = exp(lambda1 * x1)
  A = (lambda12 - pprod + pprod * lambda1 * x1) / C
  B = (pprod * lambda2 + pprod ^ 2 * x1) / C
  D = lambda2 + pprod * x1
  wExp = A / D
  wGamma = B / D ^ 2
  data.frame(x1, x2 = rgamma(n, (runif(n) > wExp / (wExp + wGamma)) + 1, D))
}

data = rGBVE(n=100, lambda1 = 1.2, lambda2 = 1.4, theta = 0.5)
colnames(data) = c("X1", "X2")

My goal is to find MLE for lambda1, lambda2 and theta using Optim() in r.

Kindly assist me to implement my likelihood and log-likelihood function in r. Thank you.

Solution

Your concern appears to be about the warning message

In log(A+B+C+D): NaNs produced

Such warnings are usually harmless — it just means that the optimization algorithm tried a set of parameters somewhere along the way that violated the condition A+B+C+D ≥ 0. Since these are reasonably complex expressions it would take a little bit of effort to figure out how one might constrain the parameters (or reparameterize the function, e.g. fitting some of the parameters on the log scale) to avoid the warning, but taking a guess that keeping the parameters non-negative will help, we can try using the L-BFGS-B algorithm (which is the only algorithm available in optim() that allows multidimensional bounded optimization).

r1 <- optim(par = c(1,2,1),
      fn = log.likelihood.x,
      dat  = data)

r2 <- optim(par = c(1,2,1),
      fn = log.likelihood.x,
      lower = rep(0,3),
      method = "L-BFGS-B",
      dat  = data)

The second does not generate warnings, and the results are close (if not identical):

all.equal(r1$par, r2$par)
## "Mean relative difference: 0.0001451953"

You might want to use bbmle, which has some additional features for likelihood modeling:

library(bbmle)
fwrap <- function(x) log.likelihood.x(x, dat = data)
parnames(fwrap) <- c("lambda1", "lambda2", "theta")
m1 <- mle2(fwrap, start = c(lambda1 = 1, lambda2 = 2, theta = 1), vecpar = TRUE,
           method = "L-BFGS-B", lower = c(0, 0, -0.5))
pp <- profile(m1)
plot(pp)
confint(pp)
confint(m1, method = "quad")