[R] find parameters for a gamma distribution
Prof Brian Ripley
ripley at stats.ox.ac.uk
Wed Jan 5 23:25:53 CET 2005
On Wed, 5 Jan 2005, Prof Brian Ripley wrote:
> 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)
[In case it confuses you, that got garbled in a cut-and-paste:
0.0253072 is the rate = 1/mean = 1/(scale*shape) as the default link for
the Gamma glm is reciprocal.
]
>
> 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
>
--
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
More information about the R-help
mailing list