14
$\begingroup$

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.

$\endgroup$
  • 5
    $\begingroup$ The Weibull is a continuous distribution so that the probability of getting exactly zero has probability zero. If you're getting many zeros in your data, that's an immediate clue that the Weibull is inappropriate. At any rate, your data look like count data (or at least, are discrete) and so a Weibull is probably not the best choice. $\endgroup$ – cardinal Dec 15 '11 at 14:10
  • $\begingroup$ Adding some context as to where the data came from will help anyone trying to answer tremendously. $\endgroup$ – cardinal Dec 15 '11 at 14:11
9
$\begingroup$

(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
$\endgroup$
4
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ The OP explained that a previous researcher chose to fit a Weibull distribution, even though the data are real-world counts -- non-negative integer counts of the number of occurrences of something. It's unclear how your three cases relate to such a situation. $\endgroup$ – r.e.s. Dec 15 '11 at 18:46
  • $\begingroup$ Oh, good note! Fitting to the Weibull distribution is egregiously wrong. It has continuous support and is never used to model counts but survival times. Negative binomial distributions would be a sort of equivalent two parameter distribution for modeling counts, which of course depends on the nature of the data-generating process (of which we have 0 information, as the problem is stated). Thanks for pointing that out to me. $\endgroup$ – AdamO Dec 15 '11 at 21:24
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

[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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.