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@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
______________________________________________
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.
