Uncertainty propagation for the solution of an integral equation I have a dataset and I use Maximum Likelihood Estimation to estimate the values of the parameters of a weibull distribution. The MLE theory provides with theoretical Confidence Intervals (asymptotical, or for $n$ samples).
Then, I use the fitted Weibull distribution in an expression which is currently optimised numerically :
$Y(t_0) = h(t_0) . \int_{0}^{t_0} S(t) dt + S(t_0)$
Where $t_0$ is unknown and $h$ and $S$ are the hazard function and the survival function of the distribution, and therefore are functions of the parameters.
I would like to propagate uncertainty on the fitted weibull parameters to estimate confidence intervales or quantiles for Y(t_0), how could I do that (numerically or analytically) ?
Thanks !
 A: Let's break this down into easier problems.  To keep the post reasonably short, I will only sketch a good confidence interval procedure without going into all the details.
What is interesting about this situation is that because $Y$ varies in such a complex, nonlinear fashion with the distribution parameters, a careful analysis and special solution are needed to obtain valid, unbiased confidence intervals.

The Weibull model and parameterization
To begin, we have to get into the details of the model because we need to know how $Y$ depends on the parameter estimates.  The basic Weibull distribution of shape $k\gt 0$ is determined by the survival function
$$S(x;k) = \exp(-(x^k)),\quad x \ge 0.$$
It extends to a family of distributions by introducing a rate parameter $\theta\gt 0$ to multiply $x:$
$$S(x;k,\theta) = S(\theta x;k).$$
Its hazard function is defined as the negative logarithmic derivative of $S,$
$$h(x;k,\theta) = -\frac{\mathrm{d}}{\mathrm{d}x}\log S(x;k,\theta) = \frac{\mathrm{d}}{\mathrm{d}x} (\theta x)^k = k\, \theta^k x^{k-1},$$
a particularly simple form.
The integral
Thus, the integral in the question is
$$s(x,k,\theta)=\int_0^x S(t;k,\theta)\,\mathrm{d}t = \int_0^x \exp(-(\theta t)^k)\,\mathrm{d}t$$
which we may integrate via the (strictly increasing) substitution $t=(u/\theta)^{1/k},$ $\mathrm{d}t=\theta^{-1/k}u^{1/k-1}\mathrm{d}u/k:$
$$s(x,k,\theta) = \int_0^u \exp(-u)\,\theta^{-1/k}u^{1/k-1}\mathrm{d}u/k=\frac{1}{k\,\theta^{1/k}}\Gamma\left(\frac{1}{k}, (\theta x)^k\right).$$
$\Gamma$ is the incomplete Gamma function, widely available in statistical software as a multiple of the Gamma CDF of shape $1/k.$
An explicit representation of $Y$
The foregoing results yield
$$\begin{aligned}
Y(x;k,\theta) &= h(x;k,\theta) s(x;k,\theta) + S(x;k,\theta)\\
&= \theta^{k-1/k}x^{k-1} \Gamma\left(\frac{1}{k}, (\theta x)^k\right) + \exp(-(\theta x)^k).
\end{aligned}$$

This example for $x=2$ shows $Y$ may have a saddle point.  Here, that point is near $(k,\theta)=(1.7, 0.6).$  For this reason I chose to study these particular parameter values in detail below.
A confidence interval for $Y$
At this point the situation gets complicated because

*

*$Y$ is a function of two parameters, not just a transformation of one.


*Even when you fix one of the parameters, $Y$ is not necessary a one-to-one transformation of the other.
What we can do is explore the values of $Y$ that are consistent with the data.   What that means is variations in the parameters $(k,\theta)$ can only decrease the likelihood of the data.  When they decrease it too much (more about that in an instant), their combined values have to be considered implausible.
Theory (based on the asymptotic distribution of the log likelihood) says that when you allow $p$ parameters to vary, you should allow the log likelihood to decrease by up to one-half a percentile of a $\chi^2(p)$ distribution: anything smaller is implausible.  Doing this determines a region (in the parameter space, a subset of $\mathbb{R}^p$) called a confidence set.  The confidence level of this confidence set is the chosen percentile.  For instance, for 95% confidence with $p=1$ parameter you would let the log likelihood to fall by up to $1.92$ because there is a 95% chance that a $\chi^2(1)$ variable will be $2\times 1.92 = 3.84$ or less.  When varying $p=2$ parameters simultaneously, you would let the log likelihood fall by up to $3.0.$
Because $Y$ cannot necessarily be used as a parameter, we must vary the two parameters $k$ and $\theta$ to explore how the log likelihood depends on them, while examining the range of values of $Y=Y(x,k,\theta)$ that arise within the confidence region.  But what value should we use for $p:$ $1$ to reflect our focus on a single value $Y$ or $2$ to reflect the need to vary two parameters?
Simulations indicate the right value may be neither.  I studied the case $k=1.7,$ $\theta=0.6,$ $x=2$ intensively.  For sample sizes of $51$ and $300$ I found that assuming $p=1$ produces an interval for $Y$ having around $92\%$ confidence.  Here is a plot of the intervals for 500 datasets of $51$ observations each:

The true value of $Y$ is marked with a horizontal axis at $1.456.$  The datasets sorted by the lengths of the confidence intervals they produced.  Estimated values of $Y$ are shown with dots (which tend to be near the upper ends of the confidence intervals).  Intervals that do not cover $Y$ are shown in red.  There are too many of them and they tend to be biased low.  (This bias persists with sample sizes of $300.$)
Assuming $p=2$ produces an interval having around 98% confidence (based on the same simulated datasets):

(Notice the change of scale on the vertical axis.)
Now there aren't enough red intervals: if you set $p=2,$ your procedure will have higher confidence than you want.  (That's not a good thing, because it implies you spent too much to obtain your data.  Roughly, the sample size is $40\%$ greater than needed to achieve a decision procedure that meets your requirements.)
A solution: bootstrapping
These potential problems with bias (in the estimates of $Y$ and in the confidence interval coverage) suggest bootstrapping the confidence interval.  Two forms of bootstrap are attractive: the usual nonparametric method in which the data are resampled from the raw dataset and a parametric method in which the data are sampled from the distribution defined by the Maximum Likelihood parameter estimates.
I experimented with both methods, but recommend the parametric method because it is likelier to compensate well for the bias in using MLEs in the first place.

This is the default output of the boot::boot function in R after $50000$ parametric iterations.  The original dataset consisted of $300$ observations this time.  "$t$" is the bootstrap value of $Y.$  The skewed bootstrap distribution shown here indicates the desirability of the bias correction.
To summarize,

*

*The Maximum Likelihood nominal $95\%$ confidence intervals are $[1.431, 1.459]$ ($p=1$) and $[1.423, 1.462]$ ($p=2$).  Remember, though, that the former is likely too short and the latter too long.


*The 95% BCa (bias corrected and accelerated) confidence interval was estimated from these results as $[1.453, 1.497].$  This interval is shifted noticeably higher than the MLE intervals.  This is (mild) confirmation of the expectation that bootstrapping will remove at least some of the bias in the MLE estimator.


*Unfortunately, BCa intervals tend to be "unstable" in the sense that they often use extreme quantiles of the bootstrap distribution.  Three other bootstrap intervals ("Normal", "Basic", and "Percentile") run from $1.446$ to $1.449$ on the lower end to $1.469$ on the upper end.  These, too, are shifted but not by as much.  They are also narrower than the MLE intervals.  If this pattern persists, narrower intervals are good: they provide more precision.
One could identify which interval is best to use via simulation, as in the first two figures above, but since this would require days of computation, I haven't bothered.

Bootstrapping code
#
# The log likelihood for data array `x`, as a function of the shape parameter `k` 
# and the log of the rate parameter.  (Log rates or scales are better estimation targets
#  than the rates or scales themselves.)
#
Lambda <- function(beta, x) sum(dweibull(x, beta[1], exp(-beta[2]), log=TRUE))
#
# `Y` as a function of the shape parameter `k`, rate parameter `theta`, and
# data vector `x`.
#
Y <- function(k, theta, x) {
  z <- (k - 1/k) * log(theta) +
    (k-1) * log(x) + 
    pgamma((x*theta)^k, 1/k, log.p=TRUE) + lgamma(1/k) - log(k) - log(theta)
  exp(z) + exp(-(theta * x)^k)
}
#
# A synthetic dataset.
#
k <- 1.7
theta <- 0.6
t0 <- 2                # Endpoint of integral defining `Y`
print(Y(k, theta, t0)) # True value of `Y`

n <- 300
set.seed(17)
x <- rweibull(n, k, 1/theta)
fit <- maxLik(Lambda, start=c(1, 0), x=x)
#
# The maximum likelihood estimates.
#
k.hat <- coefficients(fit)[1]
theta.hat <- exp(-coefficients(fit)[2])
print(Y(k.hat, theta.hat, t0)) # MLE of `Y`
#
# The function to bootstrap.
#
f <- function(ds, i, method="Parametric") {
  if (method=="Parametric") {
    x <- rweibull(length(i), k.hat, 1/theta.hat) # Parametric
  } else {
    x <- ds[i] # Nonparametric
  }
  fit <- maxLik(Lambda, start=c(1, 0), x=x)
  Y(coefficients(fit)[1], exp(-coefficients(fit)[2]), t0)
}
#
# The bootstrap.
# (Requires perhaps 0.005 sec per iteration.)
#
library(boot)
B <- boot(x, f, 5e4)
plot(B)
boot.ci(B)           # Prints four CIs for comparison

