CODEWITHSUNDEEP
rmle
Travis
Travis

Reputation: 1241

get MLE of gamma distribution parameters(especially location parameter) in R

Gamma distribution p.d.f

Hi, I want to estimate gamma distribution parameters hand by hand! I know a lot of R functions to estimate shape and scale parameters, but it seems hard to find code about estimating location parameter.

x <- c(108,91,62,59,84,60,71,105,70,69,66,65,78,83,82,68,107,68,68,69,80,
       75,89,68,64,68,70,57,62,87,51,55,56,57,75,98,60,68,81,47,76,48,63,
       58,40,62,61,58,38,40,45,68,56,64,49,53,50,39,54,47,37,50,54,70,49,
       57,52,47,43,52,57,46,63,56,50,51,50,42,46,56,52,59,45,50,59,44,52,
       54,53,63,45,56,55,53,56,46,45,49,63,50,41,42,53,50,58,50,37,53,58,
       49,53,51,64,44,53,53,55,43,50,60,51,55,56,52,51,45,49,51,63,48,51,
       60,45,40,50,66,62,69,53,54,49,47,63,55,62,57,58,51,50,57,62,45,47,
       52,35,41,53,48,59,45,41,52,36,84,62,31,41,48,47,50,50,57,53,37,46,
       41,56,51,39,59,53,51,49,45,42,32,55,34,43,35,48,33,41,38,57,37,40,
       34,44,43,62,36,41,51,48,31,28,33,35,48,31)
# estimate shape and scale parameter
gamma_likelihood <- function(para){
  sum (  (para[2] -1)*log(x) - para[2]*log(para[1]) - log(gamma(para[2])) - x/para[1] + 1/para[1])
}

MLE = optim(c(10,10), 
            fn = gamma_likelihood, 
            method = "L-BFGS-B", 
            lower = 0.00001, 
            control = list(fnscale = -1), 
            hessian = T 
)
MLE$par


# estimate location, shape and scale parameter
gamma_likelihood <- function(para){
  x = x[x > para[1]]
  sum (  (para[3] -1)*log(x - para[1]) - para[3]*log(para[2]) - 
           log(gamma(para[3])) - x/para[2] + para[1]/para[2] )
}

MLE = optim(c(23,6,7), 
            fn = gamma_likelihood,
            method = 'L-BFGS-B',
            lower = 0.00000001,
            control = list(fnscale = -1)
)
MLE$par

This is my code, I can estimate shape and scale parameters.

However, when it comes to add location parameters into log likelihood. The result seems incorrect.The TRUE parameters are c(21.4, 5.47, 6.0).

Upvotes: 1

Views: 2790

Answers (1)

VFreguglia
VFreguglia

Reputation: 2311

If you have any observed value less or equal than your location parameter, your whole likelihood for that value of lambda must be 0 (remember it's a function of parameters, not observations).

x = x[x > para[1]] is cutting observations that don't make sense for a specific location parameter, making your function return a valid number, when it should return -Inf if any of the x is "invalid", since you'd have 0 likelihood.

Here's a corrected version of your log-likelihood function:

# estimate location, shape and scale parameter
gamma_likelihood <- function(para){
  if(min(x) < para[1]) return(-Inf)
  sum (  (para[3] -1)*log(x - para[1]) - para[3]*log(para[2]) - 
           log(gamma(para[3])) - x/para[2] + para[1]/para[2] )
}

MLE = optim(c(23,6,7), 
            fn = gamma_likelihood,
            method = 'L-BFGS-B',
            lower = 0.00000001,
            control = list(fnscale = -1)
)
MLE$par

retults in: [1] 21.161109 5.394343 6.136862

Upvotes: 1

Related Questions