let's consider some random variables collected in the vector $X$ following the distribution $f_X(X)$. We want to compute the probability that:
$$ p = \textrm{Pr} [G(X) < 0] $$
where $G(X)$ is some arbitrary function. In general, this probability can be computed by pure Monte Carlo simulation. Just draw $N$ samples and count how many lie beyond the failure boundary $G(X) = 0$:
$$ \hat{p} = \frac{N_{G<0}}{N} $$
The variance of this estimate is given by: $$ \textrm{Var}[\hat{p}] = \frac{p(1-p)}{N} \approx \frac{\hat{p}(1-\hat{p})}{N} $$
In our problem, we substitute the expensive-to-evaluate function $G(X)$ with a surrogate model $\hat{G}(X)$ which means that the number of failure samples $N_{G<0}$ might not be correct, i.e. it has its own variance (which we can compute). Using the Delta method, the variance of the variance of the probability estimate can be computed:
$$ \textrm{Var}[\textrm{Var}[\hat{p}]] $$
Now the final question: How can we use this to determine some confidence interval of $\hat{p}$?
Thanks so much!