Optimization using the optim function in R with a two parameter exponential distribution

Question

I'm having trouble trying to optimize a two-parameter exponential distribution, by finding the maximum likelihood function and then using the function optim() in R

log.lik.exp <- function(x, parametros){ 
n<-length(x)
y <- -1 *(n*log(parametros[1])-parametros[1]*(sum(x)-n*parametros[2]))
}
lambda<-2
gamma<-1
mydata <- rand.exp(100, gamma, 1/lambda) 
optim(par=c(1,1),fn=log.lik.exp,x=mydata)$par

I'm also trying to use the function persp() to build a 3d plot to get a better look at the maximum values but I keep getting errors.

Thank you in advance, G.

Answer 1

As has been pointed out already, the issue is that you are not taking the bounds into consideration. You can fix this by adding an if-statement to your log likelihood function.

log.lik.exp <- function(x, parametros){ 
   n<-length(x)
   if(min(x) < parametros[2]){
      return(-Inf)
   }
   y <- n*log(parametros[1]) - parametros[1]*sum(x-parametros[2])
}
lambda <- 2
gamma <- 1
mydata <- rand.exp(100, gamma, 1/lambda) 
optim(par=c(1,1), fn=log.lik.exp, x=mydata, control=list(fnscale=-1))$par

I have also added an argument to optim() so that you are maximizing the log-likelihood rather than minimizing the negative log-likelihood. These are of course equivalent, but intuitively I prefer the first.

While we're at it, technically you could get an error if optim proposes a negative value for $\lambda$. This can be fixed in many ways. The simplest of which is to optimize over $\lambda_\star = log(\lambda)$ which can take any real value. Then just back transform, i.e. $\lambda = e^{\lambda_\star}$ at the end.