How to fit a Weibull distribution to input data containing zeroes?

Question

I'm trying to reproduce an existing prediction algorithm, handed down by a retired researcher. The first step is to fit some observed data to a Weibull distribution, to obtain a shape and scale which will be used for predicting future values. I'm using R to do this. Here's an example of my code:

x<-c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,51,77,78,144,34,29,45,16,15,37,218,170,44,121)
f<-fitdistr(x, 'weibull')

This works fine unless there are any zeroes in the input array, which causes it to fail completely. The same thing happens in SAS. As I understand it, this is because one of the steps in calculating the Weibull distribution is taking the natural log, which is undefined for 0. Is there a reasonable way to work around this?

The best I've found so far is to add 1 to all of my input values, fit the curve, and then subtract one from my predicted values ("shift" the curve up and then back down by 1). This fits the previously predicted data fairly well, but seems like it must be a wrong way of doing so.

edit: The values in the input array are observed, real-world data (the number of occurrences of something) for a range of years. So in some years the number of occurrences was zero. Whether it's the best way or not (I agree that it may not be), the original algorithm author claims to have used the Weibull distribution, and I have to try to replicate their process.

Answer 1

(As others have pointed out, a Weibull distribution is not likely to be an appropriate approximation when the data are integers only. The following is intended just to help you determine what that previous researcher did, rightly or wrongly.)

There are several alternative methods that are not affected by zeros in the data, such as using various method-of moments estimators. These typically require numerical solution of equations involving the gamma function, because the moments of the Weibull distribution are given in terms of this function. I'm not familiar with R, but here's a Sage program that illustrates one of the simpler methods -- maybe it can be adapted to R? (You can read about this and other such methods in, e.g., "The Weibull distribution: a handbook" by Horst Rinne, p. 455ff -- however, there is a typo in his eq.12.4b, as the '-1' is redundant).

"""
Blischke-Scheuer method-of-moments estimation of (a,b)
for the Weibull distribution F(t) = 1 - exp(-(t/a)^b)
""" 

x = [23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,16,15,37,218,170,44,121]
xbar = mean(x)
varx = variance(x)
var("b"); f(b) = gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2
bhat = find_root(f, 0.01, 100)
ahat = xbar/gamma(1+1/bhat)
print "Estimates: (ahat, bhat) = ", (ahat, bhat)

This produced the output

Estimates: (ahat, bhat) =  (81.316784310814455, 1.3811394719075942)

---

If the above data are modified (just for illustration) by replacing the three smallest values by $0$, i.e.

x = [23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121]

then the same procedure produces the output

Estimates: (ahat, bhat) =  (78.479354097488923, 1.2938352346035282)

---

EDIT: I just installed R to give it a try. At the risk of making this answer over-long, for anyone interested here's my R-code for the Blischke-Scheuer method:

fit_weibull <- function(x)
{
    xbar <- mean(x)
    varx <- var(x)
    f <- function(b){return(gamma(1+2/b)/gamma(1+1/b)^2 - 1 - varx/xbar^2)}
    bhat <- uniroot(f,c(0.02,50))$root
    ahat <- xbar/gamma(1+1/bhat)
    return(c(ahat,bhat))
}

This reproduces (to five significant digits) the two Sage examples above:

x <- c(23,19,37,38,40,36,172,48,113,90,54,104,90,54,157,
     51,77,78,144,34,29,45,16,15,37,218,170,44,121)
fit_weibull(x)
[1] 81.316840  1.381145

x <- c(23,0,37,38,40,36,172,48,113,90,54,104,90,54,157,
      51,77,78,144,34,29,45,0,0,37,218,170,44,121)
fit_weibull(x)
[1] 78.479180  1.293821

Answer 2

It should fail, you should be grateful that it failed.

Your observations showed that failures occurred at the very moment you started observing them. If this is a real process, coming from real (and not simulated data), you need to somehow account for the reason why you're getting zeros. I've seen survival studies where 0 times show up as a consequence of one of several things:

1. The data are actually truncated: objects were at risk and failed before the study started and you want to pretend you had observed them all along.
2. The instruments are poorly calibrated: you don't have enough measurement precision for the study and so failures occurring near the start time were coded as exactly zero.
3. The thing coded as a zero is not a zero. They're people or objects that were excluded from the analysis one way or another. The zero just shows up in the data as a consequence of merging, sorting, or otherwise recoding missing values.

So for case 1: you need to use proper censoring methods, even if that means retrospectively pulling records. Case 2 means you can use the EM algorithm because you have a precision issue. Bayesian methods work similarly here as well. Case 3 means you just need to exclude the values that were supposed to be missing.

Answer 3

You could also try fitting a three-parameter Weibull, where the third parameter is a location parameter, let us say $\theta$. This amounts to estimating the constant that you ought to add to the data to get you the best fit to the Weibull. You might do this using a profile likelihood approach by putting a "wrapper" around fitdistr, where the wrapper takes a value of $\theta$ and the data, adds $\theta$ to the data, calls the fitdistr function, and returns the associated logliklihood:

foo <- function(theta, x)
{
  if (theta <= -min(x)) return(Inf);
  f <- fitdistr(x+theta, 'weibull')
  -2*f$loglik
}

Then minimize this function using one-dimensional optimization:

bar <- optimize(foo, lower=-min(x)+0.001, upper=-min(x)+10, x=x)

where I have just made up the "+10" based on nothing at all.

For the data with the three smallest values replaced with zeroes, we get:

> bar
$minimum
[1] 2.878442

$objective
[1] 306.2792

> fitdistr(x+bar$minimum, 'weibull')
     shape        scale   
   1.2836432   81.1678283 
 ( 0.1918654) (12.3101211)
> 

bar$minimum is the MLE of $\theta$, and the fitdistr outputs are the MLEs of the Weibull parameters, jointly with $\theta$ that is. As you can see, they are pretty close to the method-of-moments estimators @r.e.s. demonstrated above.

Answer 4

I agree with cardinal's answer above. However, it is also quite common to add a constant to avoid zeros. Another value commonly used is 0.5, but any positive constant might have been used. You might try a range of values to see if you can identify the exact value used by the previous researcher. Then you could be confident that you are able to reproduce his results, before going on a search for a better distribution.

Answer 5

[Assuming Weibull is appropriate] Johnson Kotz and Balakrishnan's book has a lot of ways to estimate Weibull parameters. Some of these do not depend on the data not including zeroes (e.g. using the mean and standard deviation, or using certain percentiles).

Johnson, N. L., Kotz, S., and Balakrishnan, N. (1994). Continuous Univariate Distributions. New York: Wiley, roughly on page 632.