4

I just simulated 100 randoms observations from a gamma density with alpha(shape parameter)=5 and lambda(rate parameter)=5 :

x=rgamma(100,shape=5,rate=5)

Now, I want to fin the maximum likelihood estimations of alpha and lambda with a function that would return both of parameters and that use these observations.

Any hints would be appreciate. Thank you.

Mercier
  • 91
  • 3
  • 10

1 Answers1

6

You can use fitdistr(...) for this in the MASS package.

set.seed(1)   # for reproducible example
x <- rgamma(100,shape=5,rate=5)

library(MASS)
fitdistr(x, "gamma", start=list(shape=1, rate=1))$estimate
#    shape     rate 
# 6.603328 6.697338

Notice that with a small sample like this you don't get great estimates.

x <- rgamma(10000,shape=5,rate=5)
library(MASS)    # may be loaded by default
fitdistr(x, "gamma", start=list(shape=1, rate=1))$estimate
#    shape     rate 
# 4.984220 4.971021

fitdistr(...) also returns the standard error of the estimates and the log-likelihood.

jlhoward
  • 58,004
  • 7
  • 97
  • 140
  • Thank you sir! Very well explained. – Mercier Sep 14 '15 at 13:12
  • For the standard error of the estimates, is it the square root of the asymptotic variance? If not, is there an interesting way to find it? I am looking forwar the function optim in R to do that. – Mercier Sep 14 '15 at 13:30
  • 1
    One way to get at this is to type `fitdistr` at the commend line (no `?` and no `(...)`, just the function name). This will list the code for that function. From this you can see that `fitdistr(...)` calculates the se as the sqrt of the diagonal of a variance-covariance matrix, which in turn is the inverse of the hessian returned by `optim(...)`. So I would say, yes, it is the square root of the asymptotic variance. – jlhoward Sep 14 '15 at 13:38