Reputation: 19
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?
Upvotes: 2
Views: 1193
Reputation: 226881
A couple of issues here:
optim()
minimizes, so you either need control=list(fnscale=-1))
or to redefine your function as the negative log-likelihoodoptim(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))
Upvotes: 2
Reputation: 2368
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)
Upvotes: 2