Prediction interval for a Weibull distribution

Question

Suppose I am producing units and have tested the failure time of some proportion of them (say 10%), possibly with some right censoring if they didn't fail within the testing period.

I'm able to fit a Weibull distribution that models the failure time of my units fairly accurately, using MLE to get estimates for the shape and scale parameters. This lets me say, for example, that 8% of my units will fail during the second month of operation.

But I'd like to put a prediction interval around that prediction, to be able to say that 90% of the time, between 3% and 12% of the units will fail during that time. Is there a formula or method for computing these?

(What I've found before asking: I've found formulas for computing confidence intervals. These are much too tight. Far more than 10% of the data falls outside the confidence interval in practice. I've also found papers that discuss prediction intervals around time-to-failure, but what I'm looking for is a pointwise confidence interval for reliability at each possible time interval.)

Answer 1

If I understand your question correctly, what you want is an interval containing $$ \theta=P(t_1 < T \le t_2)=e^{-(t_1/a)^b}-e^{-(t_2/a)^b}, $$ with a nominal probability $1-\alpha$ for given values of $t_1$ and $t_2$. Since $\theta$ is a parameter and not a random variable, this interval will be a confidence and not a prediction interval.

One approach would be to fit the model using e.g. MASS::fitdistr, compute the approximate asymptotic standard error of $\theta$ using the delta method, and then construct an approximate confidence interval for $\theta$ based on asymptotic normality of $\hat\theta$.

An approach that usually works better is be to construct the confidence interval based on the profile log likelihood. After reparameterising the model in terms of $\theta$ and perhaps $a$, the profile log likelihood for $\theta$ would be given by $$ l_p(\theta) = l(\theta,\hat a(\theta)), $$ where $l$ is the log likelihood and $\hat a(\theta)$ is the MLE of $a$ for a given value of $\theta$. I think you'll need to compute $\hat a(\theta)$ numerically. The approximate $1-\alpha$ confidence limits are then equal to the values of $\theta$ at which $l_p(\theta)$ falls $\frac{\chi_{1,\alpha}^2}2$ below the maximum log likelihood.