1
votes
nll <- function(lambda, kappa){
  logit=function(x) {log(x/(1-x))}
  a=c(1-exp(-(15/lambda)^kappa), 1-exp(-(25/lambda)^kappa), 1-exp(-(35/lambda)^kappa))
  a=logit(a)
  mu = c(0.1, 0.2, 0.3)
  mu = logit(mu)
  cov = matrix(c(0.18830690, 0.00235681, 0.00071954, 0.00235681, 0.00736811, 0.00110457, 0.00071954, 0.00110457, 0.00423955), nrow =3)
  L1 = dmvnorm(a, mu, cov)

  a=c(1-exp(-(25/lambda)^kappa), 1-exp(-(35/lambda)^kappa), 1-exp(-(45/lambda)^kappa))
  a=logit(a)
  mu = c(0.4, 0.1, 0.9)
  mu = logit(mu)
  cov = matrix(c(2.7595442, 0.0045178, 0.0010505, 0.0045178, 0.00972309, 0.0015120, 0.0010505, 0.0015120, 0.0088425), nrow =3)
  L2 = dmvnorm(a, mu, cov)
  -sum(log(L1*L2))
}


> mle(nll, start = list(lambda = 1, kappa = 1))
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  initial value in 'vmmin' is not finite

I'm trying to find the lambda and kappa values that maximize the function above.

My original likelihood function returns L1*L2, but because the mle function requires the negative log-likelihood to be passed in, I modified the function to return -sum(log(L1*L2)) instead.

However, I ran into the above error. And I've also tried specifying dmvnorm(... ,log = TRUE) but that didn't solve the problem.

1
your variance-covariance matrices look fishy. If they're all constant then they will be non-positive-definite (i.e. not legal variance-covariance matrices) and dmvnorm will return a non-finite value ...Ben Bolker
@BenBolker I changed my covariance matrices. But the error still persists?Adrian

1 Answers

0
votes
  • L1 and L2 are both scalars. Assuming we're going to pass log=TRUE to dmvnorm so they are each log-likelihoods, do you mean just -(L1+L2) in the final line?
  • by specifying debug(nll) and nll(lambda=1,kappa=1) , then stepping through the code, waiting til we find an infinite value, and then backtracking, we see that 1-exp(-(45/lambda)^kappa) is exactly 1 (exp(-45) is less than 1e-16, the smallest value for which 1+x is > 1, so that the final element of logit(a) is infinite ...

So if I make dmvnorm(...,log=TRUE) in both places, change the last line to return(-(L1+L2)), and change the initial value of lambda to 10, I get a finite value for nll(10,1) (4474), and stats4::mle(nll,start=list(lambda=10,kappa=1)) gives:

Coefficients:
   lambda     kappa 
40.622673  4.883857