estimating lambda for a exponential distribution using method of MLE

2
votes

I am attempting to estimate lambda using the method of maximum likelihood estimation.

Here is the code I am using:

set.seed(0)
x=rexp(100,10)
plot(x)
fn <- function(lambda){
   length(x)*log(lambda)-lambda*sum(x)
 }
plot(fn)
optim(lambda, fn)

Error in optim(lambda, fn) : object 'lambda' not found nlm(fn, lambda) Error in nlm(fn, lambda) : object 'lambda' not found

Clearly my issue is with the error I am not sure how to fix it. can anyone tell me how to fix this so that i can get the estimation or perhaps recommend a better method?

2
rmle

2 Answers

2
votes

Sure thing. So this is a case for the optimize function

# Build Data
set.seed(0)
x=rexp(100,10)

Modifying the equation for log-likelihood slightly we have (still numerically equivalent):

log_like_lambda <- function(lambda, x){
  length(x) * log(lambda) - lambda*length(x)*mean(x)
}

We can then use our optimize function to find the maximum.

optimize(f = log_like_lambda, x, interval = c(1,15), maximum = TRUE)

Which outputs our maximum as expected from the fake data

$maximum
[1] 9.688527

$objective
[1] 127.0944

Aside:

Note that it is also helpful to plot the log-likelihood to make sure you are optimizing what you think you are optimizing:

log_like <- log_like_lambda(1:50)
plot(log_like)

log-like-exp-func

2
votes

A couple of issues here:

  1. The first argument should be a numeric vector (of length 1 in this case)
  2. optim() minimizes, so you either need control=list(fnscale=-1)) or to redefine your function as the negative log-likelihood
optim(1,fn,control=list(fnscale=-1))

works, although it gives a warning suggesting that you should use method="Brent".

You might want to consider the fitdistr() function in the MASS package (for MLE fits to a variety of distributions), or the mle2() function in the bbmle package (for general MLE, including this case, e.g. mle2(x ~ dpois(lambda), data=data.frame(x), start=list(lambda=1))