There's a number of issues with this:
(a) your data appear to be binned counts, not measurements along a curve.
(b) The function you are trying to fit is the Weibull _density_ This has
integral 1, by definition, whereas any curve anywhere near your y's would have
integral near sum(y)=127
(c) SSweibull is for growth curves which are proportional to the cumulative
Weibull distribution.
(d) SelfStart functions do *not* need starting values, that is the whole point
So you need to study things a bit more...
The expedient way would be this:
> MASS::fitdistr(rep(x,y), "Weibull")
shape scale
2.4207659 10.5078293
( 0.1530137) ( 0.4079979)
Warning message:
In densfun(x, parm[1], parm[2], ...) : NaNs produced
> plot(y~x, ylim=c(0,20), xlim=c(0,24))
> curve(127*dweibull(x,2.42,10.5), add=TRUE)
It doesn't actually fit very well, but there are quite a few observations out
in what was supposed to be the tail of the distribution.
If you want to play with SSweibull, you might do something like
> yy <- cumsum(y)
> nls(yy~SSweibull(x, Asym, Drop, lrc, pwr))
Nonlinear regression model
model: yy ~ SSweibull(x, Asym, Drop, lrc, pwr)
data: parent.frame()
Asym Drop lrc pwr
122.417 122.471 -6.944 3.129
residual sum-of-squares: 187
This gives a nonlinear least squares fit to the cumulative distribution (I am
_not_ advocating this, but you said that you were trying to figure out what
others had been up to...). If you plot it, you get
> plot(yy~x)
> curve(SSweibull(x, 122.42, 122.47, -6.94, 3.13), add=TRUE)
which _looks_ nicer, but beware! Everything looks nicer when cumulated and the
fit strongly underemphasizes that the data curve is clearly growing past x=15.
Notice also that there is a parametrization difference. SSweibull has Asym and
Drop which are F(inf) and F(inf)-F(0) respectively; in this case one could fix
both at 127. pwr is equal to a in the Weibull density, whereas lrc (log rate
constant) comes from writing exp(-(x/b)^a) as exp(-exp(lrc)*x^a), so
b = exp(-lrc)^(1/a) -- i.e. exp(6.94)^(1/3.13) = 9.18 which is in the same
range as the estimate from fitdistr().
You could also fit the weibull density directly using least squares
> nls(y~127*dweibull(x,shape,scale), start=c(shape=3,scale=10))
Nonlinear regression model
model: y ~ 127 * dweibull(x, shape, scale)
data: parent.frame()
shape scale
3.419 9.574
residual sum-of-squares: 230.6
Number of iterations to convergence: 6
Achieved convergence tolerance: 6.037e-06
> plot(y~x, ylim=c(0,20), xlim=c(0,24))
> curve(127*dweibull(x,2.42,10.5), add=TRUE)
> curve(127*dweibull(x,3.419,9.574), add=TRUE)
This fits the peak quite a bit better than the fitdistr() version, but notice
again that there are also more observations in regions where there shouldn't
really be any according to the fitted curve. This is a generic difference
between maximum likelihood and the curve fitting approaches.
-pd
On 13 Oct 2015, at 23:42 , Aditya Bhatia <aditya.bhati...@gmail.com> wrote:
> I am trying to fit this data to a weibull distribution:
>
> My y variable is:1 1 1 4 7 20 7 14 19 15 18 3 4 1 3 1 1 1
> 1 1 1 1 1 1
>
> and x variable is:1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
> 19 20 21 22 23 24
>
> The plot looks like this:http://i.stack.imgur.com/FrIKo.png and I want
> to fit a weibull curve to it. I am using the nls function in R like
> this: nls(y ~ ((a/b) * ((x/b)^(a-1)) * exp(- (x/b)^a)))
>
> This function always throws up an error saying: Error in
> numericDeriv(form[[3L]], names(ind), env) :
> Missing value or an infinity produced when evaluating the model
> In addition: Warning message:
> In nls(y ~ ((a/b) * ((x/b)^(a - 1)) * exp(-(x/b)^a))) :
> No starting values specified for some parameters.
> Initializing ‘a’, ‘b’ to '1.'.
> Consider specifying 'start' or using a selfStart model
>
> So first I tried different starting values without any success. I
> cannot understand how to make a "good" guess at the starting values.
> Then I went with the SSweibull(x, Asym, Drop, lrc, pwr) function which
> is a selfStart function. Now the SSWeibull function expects values for
> Asym,Drop,lrc and pwr and I don't have any clue as to what those
> values might be.
>
> I would appreciate if someone could help me figure out how to proceed.
>
> Background of the data: I have taken some data from bugzilla and my
> "y" variable is number of bugs reported in a particular month and "x"
> variable is the month number after release.
>
> ______________________________________________
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd....@cbs.dk Priv: pda...@gmail.com
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
