On 06-Aug-11 19:37:49, Alexander Engelhardt wrote:
> On 08/06/2011 09:23 PM, Jorge Ivan Velez wrote:
>> Hi Alex,
>>
>> Try
>>
>>> require(MASS)
>> Loading required package: MASS
>>> b<- c(2039L, 2088L, 5966L, 2353L, 1966L, 2312L, 3305L, 2013L, 3376L,
>> + 3363L, 3567L, 4798L, 2032L, 1699L, 3001L, 2329L, 3944L, 2568L,
>> + 1699L, 4545L)
>>> fitdistr(b, 'gamma')
>> shape rate
>> 6.4528939045 0.0021887943
>> (0.7722567310) (0.0001854559)
>>
>
> Hi,
>
> I'm getting a weird error here, pasted below. Do you know what
> could cause that?
>
> > fitdistr(b,"gamma")
> Error in optim(x = c(2039L, 957L, 2088L, 5966L, 2307L, 2044L,
> 2353L, 1966L, :
> non-finite finite-difference value [2]
> In addition: Warning message:
> In dgamma(x, shape, scale, log) : NaNs produced
I'm surprised no=-one has suggested going straight down the
maximum-likelihood route. The solution for the 'rate' parameter
is immediate, while the solution for the 'shape' paramater
needs the solution of an equation in one variable imvolving
the digamma function psi(k).
Writing the density of the gamma distribution as
(1/G(k))*(x^(k-1))*exp(x/r)/(r^k)
where k is the shape parameter and r is the rate, and
G(k) denotes the gamma function Gamma(k), we have for
the estimates
[1] r = mean(x.i)
[2] psi(k) = mean(log(x.i))
where {x.1,x.2,...,x.n} is the sample, and psi(k) denotes
the digamma function (G'(k))/G(k) = (log(G(k))' (the "'"
denotes first derivative).
So you just need to solve [2] for k.
This is well summarised in the Wikipedia article:
http://en.wikipedia.org/wiki/
Gamma_distribution#Maximum_likelihood_estimation
in the Section "Parameter estimation".
Since the digamma function digamma() is in the base package,
all you need to do is submit it to a root-finder, here done
in the naivest possible way.
With your values of
x<- c(2039L, 2088L, 5966L, 2353L, 1966L, 2312L, 3305L,
2013L, 3376L, 3363L, 3567L, 4798L, 2032L, 1699L,
3001L, 2329L, 3944L, 2568L, 1699L, 4545L)
you get
mean(log(x))
# [1] 7.92335
which is quite a modest number (though in a somewhat tense
relationship with the digamma function).
A bit of experimentation with plot() will lead you home:
t<-(50:200); plot(t,digamma(t))
t<-(100:2000); plot(t,digamma(t))
t<-20*(100:200); plot(t,digamma(t))
from which the root is somewhere near 2750.
t<-(2700:2800); plot(t,digamma(t))
t<-(2745:2865); plot(t,digamma(t))
t<-2755+0.1*(0:100); plot(t,digamma(t))
So:
mean(log(x)) - digamma(2760)
# [1] 0.0005452383
mean(log(x)) - digamma(2769)
# [1] -0.002710915
mean(log(x)) - digamma(2765)
# [1] -0.001265045
mean(log(x)) - digamma(2763)
# [1] -0.0005413247
mean(log(x)) - digamma(2762)
# [1] -0.0001792682
mean(log(x)) - digamma(2761)
# [1] 0.0001829194
mean(log(x)) - digamma(2761.5)
# [1] 1.809230e-06
mean(log(x)) - digamma(2761.75)
# [1] -8.873357e-05
mean(log(x)) - digamma(2761.625)
# [1] -4.34632e-05
mean(log(x)) - digamma(2761.5625)
# [1] -2.082724e-05
mean(log(x)) - digamma(2761.53125)
# [1] -9.509068e-06
mean(log(x)) - digamma(2761.515625)
# [1] -3.849935e-06
mean(log(x)) - digamma(2761.5078125)
# [1] -1.020356e-06
mean(log(x)) - digamma(2761.50390625)
# [1] 3.944361e-07
mean(log(x)) - digamma(2761.505859375)
# [1] -3.129603e-07
mean(log(x)) - digamma(2761.5048828125)
# [1] 4.073787e-08
And so it goes ... Anyway, to 7 significant figures,
the answer is that k = 2761.505, and that is "by hand"
using a simple "interval-halving" procedure.
I don't know why fitdistr() threw an error. It may be
that the starting-value for the iterations was to far
off and threw it into some difficult space. The good
starting value above was obtained by a bit of graphical
trial and error. Maybe something similar can be emulated
in R?
Hoping this helps,
Ted.
PS
Oh, and while I'm at it, the estimate of the rate is
r = mean(x) = 2948.15
--------------------------------------------------------------------
E-Mail: (Ted Harding) <ted.hard...@wlandres.net>
Fax-to-email: +44 (0)870 094 0861
Date: 07-Aug-11 Time: 02:03:29
------------------------------ XFMail ------------------------------
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
