Forum: help

RE: maxLik-estimate parametrs of gamma distribution [ Reply ] By: Ott Toomet on 2019-11-08 17:58	[forum:47132]
Hi Vaclav -- it seems to me that your log-likelihood function is wrong. If you use something like dgamma(X, k, theta, log=TRUE) instead of the formula you use in the function it works as expected (although the results are rather noisy with only 42 observations). Cheers, Ott

maxLik-estimate parametrs of gamma distribution [ Reply ]
By: Vaclav Kalous on 2019-11-05 14:09

[forum:47125]

Hello,
I generated a dataset from exponential distribution:
a=rexp(42,0.3)
len_a=length(a)
I would like to fit it with gamma distribution, so I defined its log likelihood function:
ll_gamma_a=function(param){

k=param[1]
theta=param[2]
if (k<=0){
return(NA)
}
if (theta<=0){
return(NA)
}

return((k-1)*(sum(log(a))-sum(a/theta)-len_a*k*log(theta)-len_a*log(factorial(k-1))))
}

Now I want to maximalize it using maxLik:

mle_gamma_a= maxLik(ll_gamma_a, start = c(shape=1,rate=1))
summary(mle_gamma_a)

The output I get is this:

--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 5 iterations
Return code 3: Last step could not find a value above the current.
Boundary of parameter space?
Consider switching to a more robust optimisation method temporarily.
Log-Likelihood: 8644.971
2 free parameters
Estimates:
Estimate Std. error t value Pr(> t)
shape 0.16217947 0.01600630 10.13 <0.0000000000000002 ***
rate 0.01491913 0.00008685 171.78 <0.0000000000000002 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Which is not values I would expect i.e. something like shape=1 and rate=0.3 approximately. Can anyone please suggest what I am doing wrong? I think one of the problems is the "start" argument- I am not really sure how to define it.
Thank you very much for all suggestions!