10

I have been doing some data analysis in R and I am trying to figure out how to fit my data to a 3 parameter Weibull distribution. I found how to do it with a 2 parameter Weibull but have come up short in finding how to do it with a 3 parameter.

Here is how I fit the data using the fitdistr function from the MASS package:

y <- fitdistr(x[[6]], 'weibull')

x[[6]] is a subset of my data and y is where I am storing the result of the fitting.

Julius Vainora
  • 47,421
  • 9
  • 90
  • 102
Matthew Crews
  • 4,105
  • 7
  • 33
  • 57
  • Perhaps if you made a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) that demonstrates your question / problem, people would find it easier to answer. Specifically, what does `x[[6]]` look like. At a minimum, post `str(x[[6]]` or preferably the results of `dput(x[[6]])`. – Andrie Aug 05 '12 at 16:07
  • 2
    You can't use the builtin `weibull` distribution available in R, because it's a two parameters weibull distribution. You have to compute custom probability density function (3 parameters) and use it instead. – dickoa Aug 05 '12 at 16:17

2 Answers2

8

First, you might want to look at FAdist package. However, that is not so hard to go from rweibull3 to rweibull:

> rweibull3
function (n, shape, scale = 1, thres = 0) 
thres + rweibull(n, shape, scale)
<environment: namespace:FAdist>

and similarly from dweibull3 to dweibull

> dweibull3
function (x, shape, scale = 1, thres = 0, log = FALSE) 
dweibull(x - thres, shape, scale, log)
<environment: namespace:FAdist>

so we have this

> x <- rweibull3(200, shape = 3, scale = 1, thres = 100)
> fitdistr(x, function(x, shape, scale, thres) 
       dweibull(x-thres, shape, scale), list(shape = 0.1, scale = 1, thres = 0))
      shape          scale          thres    
    2.42498383     0.85074556   100.12372297 
 (  0.26380861) (  0.07235804) (  0.06020083)

Edit: As mentioned in the comment, there appears various warnings when trying to fit the distribution in this way

Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  non-finite finite-difference value [3]
There were 20 warnings (use warnings() to see them)
Error in optim(x = c(60.7075705026659, 60.6300379017397, 60.7669410153573,  : 
  L-BFGS-B needs finite values of 'fn'
In dweibull(x, shape, scale, log) : NaNs produced

For me at first it was only NaNs produced, and that is not the first time when I see it so I thought that it isn't so meaningful since estimates were good. After some searching it seemed to be quite popular problem and I couldn't find neither cause nor solution. One alternative could be using stats4 package and mle() function, but it seemed to have some problems too. But I can offer you to use a modified version of code by danielmedic which I have checked a few times:

thres <- 60
x <- rweibull(200, 3, 1) + thres

EPS = sqrt(.Machine$double.eps) # "epsilon" for very small numbers

llik.weibull <- function(shape, scale, thres, x)
{ 
  sum(dweibull(x - thres, shape, scale, log=T))
}

thetahat.weibull <- function(x)
{ 
  if(any(x <= 0)) stop("x values must be positive")

  toptim <- function(theta) -llik.weibull(theta[1], theta[2], theta[3], x)

  mu = mean(log(x))
  sigma2 = var(log(x))
  shape.guess = 1.2 / sqrt(sigma2)
  scale.guess = exp(mu + (0.572 / shape.guess))
  thres.guess = 1

  res = nlminb(c(shape.guess, scale.guess, thres.guess), toptim, lower=EPS)

  c(shape=res$par[1], scale=res$par[2], thres=res$par[3])
}

thetahat.weibull(x)
    shape     scale     thres 
 3.325556  1.021171 59.975470
Julius Vainora
  • 47,421
  • 9
  • 90
  • 102
  • I do that and I get an error with the following message: Error in fitdistr(x, function(x, shape, scale, thres) dweibull(x - thres, : optimization failed In addition: Warning messages: 1: In dweibull(x - thres, shape, scale) : NaNs produced 2: In dweibull(x - thres, shape, scale) : NaNs produced 3: In dweibull(x - thres, shape, scale) : NaNs produced 4: In dweibull(x - thres, shape, scale) : NaNs produced – Matthew Crews Aug 07 '12 at 02:24
  • @Wallhood, I edited the answer, now it seems to work perfectly, but unfortunately it doesn't provide information about variance. – Julius Vainora Aug 07 '12 at 09:57
  • 1
    Wow, I can not tell you how awesome that is and how grateful I am. If you are ever in Portland, Oregon I would happily buy you a beer. – Matthew Crews Aug 09 '12 at 02:24
1

An alternative: package "lmom". The estimative by L-moments technique

library(lmom)
thres <- 60
x <- rweibull(200, 3, 1) + thres
moments = samlmu(x, sort.data = TRUE)
log.moments <- samlmu( log(x), sort.data = TRUE )
weibull_3parml <- pelwei(moments)
weibull_3parml
zeta      beta     delta 
59.993075  1.015128  3.246453

But I don´t know how to do some Goodness-of-fit statistics in this package or in the solution above. Others packages you can do Goodness-of-fit statistics easily. Anyway, you can use alternatives like: ks.test or chisq.test

Marcel
  • 147
  • 4