Regression with local confidence estimation

I want to create a regression model that not only predicts the value of interest, but also outputs an estimate on how accurate that value is. With expectation that for some feature values model will be more certain about the prediction, and in another cases less.

A toy example: given that someone's height is $2$ meters, I want to predict that his European shoe size is $45 \pm 3$, where $\pm 3$ is some estimate of accuracy, say std of prediction errors that model made locally in the interval $[1.8, 2.2]$.

It seams that this kind of thing should be doable, for example, by training a regressor, computing smoothed standard deviation of test set errors, and training another regressor on that. However, I can't find any references about anyone doing that, and a question of very similar spirit has no answers.

Am I trying to do something very unusual? Otherwise, could you please give some pointers to existing literature, as I struggle to find anything myself.

• How do you intent to smooth the standard deviation? Smoothing requires some assumptions about how similar inputs should have similar errors, do you have any basis on which to formulate such assumptions? For the moment I can add that when predicting categorical variables, the results of a final soft-max can be interpreted as a confidence in certain results. Nov 9, 2016 at 14:14
• That is indeed a tricky part, and while I can image some ad-hoc strategies, I don't have a satisfactory answer to this now. I was hoping to get some insights from the (non existing?) references. Nov 9, 2016 at 14:20
• A second problem is that as soon as you train a second model on the size of the mistakes of your first model, the next natural question is: how big are the mistakes that your second model makes? Nov 9, 2016 at 14:25
• Indeed, although I was hoping for this to be less of a concern, as I don't expect standard deviation of errors to be very noisy (given I do the smoothing right). Nov 9, 2016 at 14:28
• I don't really see how smoothing is going to help you with noise in the error-rate. If the error-rate of your original regression model is noisy, smoothing those errors isn't going to change that. As I understand, you will get reasoning like: "for people with height between 1.95-2.05 meter, we can get on average the shoe size with very little error, however the shoe size of some with specifically height 1.973 meter may be way off". In the end I would only be interested in the error-rate of my actual predictions, not in the error-rate of similar inputs. Nov 9, 2016 at 14:33

Answering my own question: apparently this is indeed doable and already researched. I kind of stumbled upon things over time, but have not tried any in production as I have moved on to doing other things.

Quantile Random Forests

The simplest answer to my question is probably Quantile Random Forests. Here is a readable blog post with references to papers, as well as R and Python implementations. The rough idea is to modify the random forest a little to remember the target values in the leaves and to derive quantiles from them at prediction time. Very nice!

Uncertainty estimating neural nets

The less simple but lovely (at least for my taste) idea, is to train two neural networks - one for mean and one for variance - that fit the conditional distribution. This one I had to try to believe that it works, so here's how it goes briefly:

Take some synthetic data with $x \sim \text{Uniform}[-1, 1]$, and $y | x \sim \text{Normal}(x, |x|)$:

And try to fit with two neural networks. They need not to be special, and here is a pytorch example with both combined into one for comfort:

class TwoOutputNet(nn.Module):

def __init__(self):
super().__init__()

self.m1 = nn.Linear(1, 10)
self.m2 = nn.Linear(10, 1)
self.s1 = nn.Linear(1, 10)
self.s2 = nn.Linear(10, 1)

def forward(self, x):

m = F.tanh(self.m1(x))
s = F.tanh(self.s1(x))

return self.m2(m), torch.clamp(self.s2(s), min=-5, max=5)


The key is the loss function that takes true value (y), estimated mean (m), estimated standard deviation (s) and returns a loss that penalizes both error weighted by uncertainty and the uncertainty. For example, the negative log likelihood of the normal distribution:

class NegativeLogLikelihood(nn.Module):

def __init__(self):
super().__init__()

def forward(self, y, m, s):

sq_part = ((y - m) / torch.exp(s)) ** 2
lg_part = s

return (1/2 * sq_part + lg_part).mean()


Another point to note, is that fitting for standard deviation directly seems to be hard (convergence is erratic for me), so it is worth experimenting with different parameterizations. The pasted code takes the exponent (that's why the second output of the net is clamped between -5 and 5 -- the standard deviation gets clamped between 0.005 and 150).

In the end we get the both the estimates of the mean and standard deviation at each value of x:

This idea and different parameterizations are explored more in a recent paper. It also notes, that the parametric assumptions, such as used above, are not necessary, and the second net in general can output the expected loss of the first one, which I think is very nice.

Bayesian neural nets

And then there is the whole Bayesian neural networks thing, which I did not investigate further. I quite liked this blog post for inspiration, and I guess edwardlib and this thesis are good starting points.