Prediction interval for a Weibull distribution 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.)
 A: 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.
