Derivation of confidence and prediction intervals of predictions for probit and logit (and GLMs in general) The derivation of the prediction interval for the linear model is quite simple: Obtaining a formula for prediction limits in a linear model .
How to derive the confidence and prediction intervals for the fitted values  of the logit and probit regressions (and GLMs in general)?
 A: In GLM, prediction is a non-linear function $f$ of the product of covariates $X$ with estimated coefficient vector $\hat{\beta}$:
$$\hat{y} = f(X\hat{\beta})$$
Finite-sample distribution of $\hat{\beta}$ is generally unknown, but as long as  $\hat{\beta}$ is a maximum likelihood estimate, it has asymptotic normal distribution $\mathcal{N}(\beta, -H^{-1})$, where $H$ is the Hessian matrix of the likelihood function in its maximum. The p-values of $\beta$ that are shown as an output of a regression are nearly always based on this asymptotics. But if you feel your sample is too small for asymptotics, use numerical distribution (e.g. bootstrapping). 
When you use asymptotic normal distribution of $\hat{\beta}$ (and therefore $X\hat{\beta}$), distribution of $\hat{y}$ is still non-normal due to non-linear $f$. You can ignore it - get normal confidence bounds $(z_{lower}, z_{upper})$ for  $X\beta$, and plug them into $f$, getting bounds for $y$ as $(y_{lower}, y_{upper}) = (f(z_{lower}), f(z_{upper}))$.
Another strategy (called delta method) is to take a Taylor expansion of $f$ around $X\hat{\beta}$ - it will be linear in $\hat{\beta}$. Therefore, you can approximate distribution of $f(X\hat{\beta})$ as 
$$f(X\hat{\beta}) \sim \mathcal{N}\left(f(X\beta), -(f^{'}(X\beta))^2 X H^{-1} X^T \right)$$
Then the asymptotic 95% confidence interval for $f(X\beta)$ would look like
$$ f(X\hat{\beta}) \pm 1.96 \sqrt{(f^{'}(X\hat{\beta}))^2 X H(\hat{\beta})^{-1} X^T}$$
Now you need only to find expression for Hessian matrices for particular models, like logistic regression in this question. And this question presents practical comparison of bootstrap, transformed normal bounds, and delta method for logistic regression.
A: When all else fails, you can always construct bootstrapped CIs for any statistic. Here's a simple algorithm:


*

*Draw $N$ samples with replacement from $X$ (where $N$ is the number of rows in $X$). You'll find that about 2/3rds of your observations will appear in such a sample.

*Use these samples to fit a model

*Use this model to generate predictions for the observations in $X$ that weren't used in training.

*Repeat this process 100 or so times (the more the merrier) to accumulate a collection of predictions for each observation. This collection is an approximation to the distribution of your predictions. Call these your "bootstrapped predictions".

*Construct confidence intervals by taking quantiles on the predictions. E.g. for a particular observation, calculate the .025 and .975 quantiles for a 95% confidence interval.

