[R] find parameters for a gamma distribution

[Prof Brian Ripley]

First, you want to fit the data not the histogram counts (binned data).

Second, glm does not do a very principled fit of a gamma (it is a moment 
estimator).  Something like

> d <- rgamma(1000, 20, scale = 2)
> summary(glm(d ~ 1, Gamma))
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0253072  0.0001822   138.9   <2e-16

(Dispersion parameter for Gamma family taken to be 0.05185392)

where the shape = 1/0.05185392 (no s.e.) and 0.0253072 is the mean
= 1/(scale*shape)

Try instead (on 100000 obs here)

library(MASS)
fitdistr(d, "Gamma")
       shape           rate
   19.812444027    0.495157589
  ( 0.087858938) ( 0.002223778)

which works for me.  You can make it use scale =1/rate by
> fitdistr(d, dgamma, start=list(shape=20,scale=2))
       shape          scale
   19.812702659    2.019419771
  ( 0.087774433) ( 0.009060263)


On Wed, 5 Jan 2005, Andrew Collier wrote:

> hello,
>
> i have just started exploring R as an alternative to matlab for data analysis. so
> +far everything is _very_ promising. i have a question though regarding parameter
> +estimation. i have some data which, from a histogram plot, appears to arise from
> +a gamma distribution. i gather that you can fit the data to the distribution
> +using glm(). i am just not quite sure how this is done in practice... so here is
> +a simple example with artificial data:
>
> d <- rgamma(100000, 20, scale = 2)
> h <- hist(d, breaks = c(seq(10, 80, 2), 100))
>
> H <- data.frame(x = h$mids, y = h$density)
>
> g <- glm(y ~ x, data = H, family = Gamma)
> summary(g)
>
> Call:
> glm(formula = y ~ x, family = Gamma, data = H)
>
> Deviance Residuals:
>    Min       1Q   Median       3Q      Max
> -3.8654  -2.0887  -0.7685   0.7147   1.4508
>
> Coefficients:
>            Estimate Std. Error t value Pr(>|t|)
> (Intercept)  30.4758    26.7258   1.140    0.262
> x             1.0394     0.6825   1.523    0.137
>
> (Dispersion parameter for Gamma family taken to be 1.343021)
>
>    Null deviance: 119.51  on 35  degrees of freedom
> Residual deviance: 116.28  on 34  degrees of freedom
> AIC: -260.49
>
> Number of Fisher Scoring iterations: 7
>
> now i suppose that the estimates parameters are:
>
>        shape = 30.4758
>        scale = 1.0394
>
> am i interpreting the output correctly? and, if so, why are these estimates so
> +poor? i would, perhaps naively, expected the parameters from an artificial
> +sample like this to be pretty good.
>
> my apologies if i am doing something stupid here but my statistics capabilties
> +are rather limited!
>
> best regards,
> andrew collier.
> -- 
> Andrew B. Collier
>
> Antarctic Research Fellow                                   tel: +27 31 2601157
> Space Physics Research Institute                            fax: +27 31 2616550
> University of KwaZulu-Natal, Durban, 4041, South Africa
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595