Measuring uncertainty of a fitted calibration curve

Question

My question is: how can you compute the prediction interval of a calibration curve, just like you might for any other regression model.

A calibration curve maps estimated probabilities to empirical probabilities. A fairly simple but common method seems to be fitting a piecewise constant "histogram" curve, as shown here. However, I imagine this may not work very well when the amount of data is small. On the other hand, R seems to fit some sort of spline instead for the curve.

I don't know exactly how R fits the spline curve, but fitting estimated probabilities to ground truth labels $\{0, 1\}$ corresponding to "incorrect prediction" and "correct prediction" seems to be the obvious thing to do. The problem comes when I want to measure the uncertainty in the calibration curve: if the calibration curve says that a predicted confidence of 70% really maps to 73% empirically with an 1-std prediction interval of [0.23, 1.23], the interval is NOT an interval on the empirical probability, but rather an interval on the binary labels.

In the case of a histogram calibration curve, I can get the uncertainty on the estimated empirical probability by treating it as a bernoulli random variable and computing the variance. This solution does not carry over to the case where I use another arbitrary regression model.

Answer 1

In the typical bayesian setting, when we have some predictors $X = \{x_i\}$ (in this case, uncalibrated predictions), labels $R = \{r_i\}$ (in this case, the binary result) and model parameters $\theta$, we can define a posterior over the parameters as:

$$P(\theta | R, X) \propto P(R|X,\theta)P(\theta)$$

Where $P(\theta)$ is some simple prior over the parameters. Usually, we model $R|X,\theta \sim D(X,\theta)$ for some distribution $D(X,\theta)$. A popular choice is $D = \mathcal{N}(\mu = f(X;\theta), \sigma^2)$, or for binary tasks, $D = \text{Bernoulli}(p = f(X;\theta))$.

However, neither of these choices would be valid in our case: we want our model to define a probability distribution over calibrated predictions $Y$, not the results $R$. Suppose we define a distribution $Y \sim \text{Beta}(\alpha, \beta = f(X;\theta))$ -- I chose the beta distribution as its support is the unit interval, although other choices can be made. Then we can write

$$D = \text{Bernoulli}(p \sim Y) = \text{Bernoulli}(p \sim \text{Beta}(\alpha, \beta = f(X;\theta)))$$

Computing $P(R|X,\theta)$ now requires integrating out $Y$ in the general case, but it just so happens that in this case we can show $D = \text{Bernoulli}(p = \frac{\alpha}{\alpha + \beta})$.

Note that the MLE version of this model -- $\theta^* = \text{argmax}_\theta P(R|X,\theta)$ -- is equivalent to training the model to regress $R$ from $X$, as mentioned. Yet by specifying the beta-bernoulli distribution, we can perform full bayesian inference. Now that we have defined both the prior and likelihood terms, all which remains is to sample from the posterior (which can be done via MCMC, variational inference, or other bayesian methods depending on the exact form of $f$), and then once we have samples of $\theta$, to sample from $Y$ in order to obtain the predictive distribution.