# [R] fit a gamma pdf using Residual Sum-of-Squares

Peter Ehlers ehlers at ucalgary.ca
Mon Mar 8 18:53:18 CET 2010

```On 2010-03-08 8:24, vincent laperriere wrote:
> Hi all,
>
> I would like to fit a gamma pdf to my data using the method of RSS (Residual Sum-of-Squares). Here are the data:
>
>   x<- c(86,  90,  94,  98, 102, 106, 110, 114, 118, 122, 126, 130, 134, 138, 142, 146, 150, 154, 158, 162, 166, 170, 174)
>   y<- c(2, 5, 10, 17, 26, 60, 94, 128, 137, 128, 77, 68, 65, 60, 51, 26, 17, 9, 5, 2, 3, 7, 3)
>
> I have typed the following code, using nls method:
>
> fit<- nls(y ~ (1/((s^a)*gamma(a))*x^(a-1)*exp(-x/s)), start = c(s=3, a=75, x=86))
>

There are a couple of problems:
1) don't include a start value for x; it's not a 'parameter';
2) you're trying to fit a *density* function to data that's
clearly not normalized.

A quick check shows that your empirical curve integrates to

y[-1] * diff(x)
# 3992

This almost works:

fit<- nls(y ~ 4000*(1/((s^a)*gamma(a))*x^(a-1)*exp(-x/s)),
start = c(s=3, a=75))

but not quite; I still get an error. So let's do the right
thing and plot the data and some test fits:

plot(y ~ x)
curve(4000 * dgamma(x, shape=75, scale=3), add=TRUE)
# no good
curve(4000 * dgamma(x, shape=75, scale=1), add=TRUE)
# no good
curve(4000 * dgamma(x, shape=75, scale=1.6), add=TRUE)
# pretty good!

fit<- nls(y ~ 4000*(1/((s^a)*gamma(a))*x^(a-1)*exp(-x/s)),
start = c(s=1.6, a=75))

coef(fit)
#        s         a
# 1.399638 86.395409

xx <- seq(86, 174, length=100)
yy <- predict(fit, data.frame(x=xx))
plot(y ~ x)
lines(yy ~ xx, col='red')

-Peter Ehlers

> But I have the following message error (sorry, this is in German):
>
>
> Fehler in qr(.swts * attr(rhs, "gradient")) :
>    Dimensionen [Produkt 3] passen nicht zur L�nge des Objektes 
> Zus�tzlich: Warnmeldung:
> In .swts * attr(rhs, "gradient") : L�nge des l�ngeren Objektes
>    ist kein Vielfaches der L�nge des k�rzeren Objektes
>
> Could anyone help me with the code?
> I would greatly appreciate it.
> Sincerely yours,
> Vincent Laperri�re.
>
>
>
> 	[[alternative HTML version deleted]]
>
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help