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

[Peter Ehlers]

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
about 4000:

  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 [23]
> 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
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Ehlers
University of Calgary