What is the proper way to report confidence intervals (or std, variance) for machine learning - especially wrt the loss and the accuracy?

Question

I think this should be a basic question, if not at least a fundamental one. I think it might be a bit more nuanced than it might look at first sight - or at least it was to me.

Say we have a data set $D = \{ (x_i, y_i \}^N_{i=1}$ (say the validation set) and we want to intuitively compute the loss of our model $\hat f$ (e.g. a neural network). In classification, want to have a good accuracy on average on any point:

$$ \hat acc(D, \hat f) = \hat{ E}_{x_i, y_i \in D}[ 1\{ f(x_i) = y_i \} ] $$

so the random variable R.V. is $1\{ f(x_i) = y_i \}$. Simirly we'd compute the variance, std, 95%-confidence interval of r.v. $1\{ f(x_i) = y_i \}$ so we'd do:

$$ \hat std(D, \hat f) = \hat{Std}_{x_i, y_i \in D}[ 1\{ f(x_i) = y_i \} ] $$

Those are the two numbers that make sense to report to me. In general for any loss function we'd do the same but instead the r.v. would be $Loss(f(x_i), y_i)$ e.g. square loss for regression, cross-entropy loss for classification, etc.

My guess is that it is not correct to report the confidence interval, std of $\hat acc(D, \hat f)$ i.e. my guess is that it is wrong to compute the confidence interval of the r.v. $\hat acc(D, \hat f)$. For example, since this is the std of a mean, "the real" std can be shrunk arbitrarily by making $N$ large since it's std is $\approx \frac{ \sigma }{N}$. So this obviously seems wrong to me. It also makes sense intuitively from the throwing darts analogy - if we throw darts and get 5 means, they will on average be very similar since we are averging out the variation by considering the darts independently. This is relevant because one could try to compute the confidence interval on the validation set by getting say $K$ different batches of size $N$ and then compute the confidence interval, std, error bars with the $K$ batch errors instead of the $K*N$ individual losses (or e.g. 0-1 loss for the accuracy).

With this motivation in mind I decided to get function that give me tensors of size $[N]$ (with the loss for each example, including 0 1 values when computing the accuracy) to compute the error bars (std and confidence intervals). But when I actually used them I got very unexpected values for my std and 95% confidence intervals, opposite to what I'd expect.

Thus I am concretely wondering these things:

1. am I computing the expected loss/accuracy and the error bars (confidence interval, std) of the right random variables. I claim $L(f(x), y)$ - even for the 0-1 loss $1\{ f(x_i) = y_i \}$ - is correct, and computing the mean, ci, std of $\hat acc(D, \hat f)$ is wrong.
2. I am also wondering why my stds are larger than my 95% confidence intervals for the 0-1 (assuming 1 is indeed correct). Using the intuition of approximating normals with enough data, I'd expect that that 95% requires more confidence and should be like 1.96 standard deviations. So I expect $\sigma_{95\%} > \sigma_{std}$ but my simulations/unit tests reveal that is not true, which is confusing me.

To make the code concrete and self contained here is my code and it's output:

import torch
from torch import FloatTensor, nn, Tensor

P_CI = {0.90: 1.64,
        0.95: 1.96,
        0.98: 2.33,
        0.99: 2.58,
        }

# from uutils.torch_uu.metrics.confidence_intervals import torch_compute_confidence_interval_classification_torch

def torch_compute_confidence_interval_classification_torch(data: Tensor,
                                                   by_pass_30_data_points: bool = False,
                                                   p_confidence: float = 0.95
                                                   ) -> Tensor:
    """
    Computes CI interval
        [B] -> [1]
    According to [1] CI the confidence interval for classification error can be calculated as follows:
        error +/- const * sqrt( (error * (1 - error)) / n)

    The values for const are provided from statistics, and common values used are:
        1.64 (90%)
        1.96 (95%)
        2.33 (98%)
        2.58 (99%)
    Assumptions:
    Use of these confidence intervals makes some assumptions that you need to ensure you can meet. They are:

    Observations in the validation data set were drawn from the domain independently (e.g. they are independent and
    identically distributed).
    At least 30 observations were used to evaluate the model.
    This is based on some statistics of sampling theory that takes calculating the error of a classifier as a binomial
    distribution, that we have sufficient observations to approximate a normal distribution for the binomial
    distribution, and that via the central limit theorem that the more observations we classify, the closer we will get
    to the true, but unknown, model skill.

    Ref:
        - computed according to: https://machinelearningmastery.com/report-classifier-performance-confidence-intervals/

    todo:
        - how does it change for other types of losses
    """
    B: int = data.size(0)
    assert (data >= 0.0).all(), f'Data has to be positive for this CI to work but you have some negative value.'
    assert B >= 30 or by_pass_30_data_points, f' Not enough data for CI calc to be valid and approximate a' \
                                              f'normal, you have: {B=} but needed 30.'
    const: float = P_CI[p_confidence]
    error: Tensor = data.mean()
    val = torch.sqrt((error * (1 - error)) / B)
    ci_interval: float = const * val
    return ci_interval

def accuracy(output: Tensor,
             target: Tensor,
             topk: tuple[int] = (1,),
             reduction: str = 'mean'
             ) -> tuple[FloatTensor]:
    """
    Computes the accuracy over the k top predictions for the specified values of k
    In top-5 accuracy you give yourself credit for having the right answer
    if the right answer appears in your top five guesses.

    if reduction is
        - "none" computes the 1 or 0 topk acc for each example so each entry is a tensor of size [B]
        - "mean" (default) compute the usual topk acc where each entry is acck of size [], single value tensor

    ref:
    - https://stackoverflow.com/questions/51503851/calculate-the-accuracy-every-epoch-in-pytorch/63271002#63271002
    - https://pytorch.org/docs/stable/generated/torch.topk.html
    - https://discuss.pytorch.org/t/imagenet-example-accuracy-calculation/7840
    - https://gist.github.com/weiaicunzai/2a5ae6eac6712c70bde0630f3e76b77b
    - https://discuss.pytorch.org/t/top-k-error-calculation/48815/2
    - https://stackoverflow.com/questions/59474987/how-to-get-top-k-accuracy-in-semantic-segmentation-using-pytorch

    :param output: output is the prediction of the model e.g. scores, logits, raw y_pred before normalization or getting classes
    :param target: target is the truth
    :param topk: tuple of topk's to compute e.g. (1, 2, 5) computes top 1, top 2 and top 5.
    e.g. in top 2 it means you get a +1 if your models's top 2 predictions are in the right label.
    So if your model predicts cat, dog (0, 1) and the true label was bird (3) you get zero
    but if it were either cat or dog you'd accumulate +1 for that example.
    :return: list of topk accuracies [top1st, top2nd, ...] depending on your topk input. Size [] or [B] depending on
        reduction type.
    """
    with torch.no_grad():
        # ---- get the topk most likely labels according to your model
        # get the largest k \in [n_classes] (i.e. the number of most likely probabilities we will use)
        maxk = max(topk)  # max number labels we will consider in the right choices for out model
        batch_size = target.size(0)

        # get top maxk indicies that correspond to the most likely probability scores
        # (note _ means we don't care about the actual top maxk scores just their corresponding indicies/labels)
        _, y_pred = output.topk(k=maxk, dim=1)  # _, [B, n_classes] -> [B, maxk]
        y_pred = y_pred.t()  # [B, maxk] -> [maxk, B] Expects input to be <= 2-D tensor and transposes dimensions 0 and 1.

        # - get the credit for each example if the models predictions is in maxk values (main crux of code)
        # for any example, the model will get credit if it's prediction matches the ground truth
        # for each example we compare if the model's best prediction matches the truth. If yes we get an entry of 1.
        # if the k'th top answer of the model matches the truth we get 1.
        # Note: this for any example in batch we can only ever get 1 match (so we never overestimate accuracy <1)
        target_reshaped = target.view(1, -1).expand_as(y_pred)  # [B] -> [B, 1] -> [maxk, B]
        # compare every topk's model prediction with the ground truth & give credit if any matches the ground truth
        correct = (
                y_pred == target_reshaped)  # [maxk, B] were for each example we know which topk prediction matched truth
        # original: correct = pred.eq(target.view(1, -1).expand_as(pred))

        # -- get topk accuracy
        list_topk_accs = []  # idx is topk1, topk2, ... etc
        for k in topk:
            # get tensor of which topk answer was right
            ind_which_topk_matched_truth = correct[:k]  # [maxk, B] -> [k, B]
            # accuracy for the current topk for the whole batch,  [k, B] -> [B]
            indicator_which_topk_matched_truth = ind_which_topk_matched_truth.float().sum(dim=0)
            assert indicator_which_topk_matched_truth.size() == torch.Size([batch_size])
            # put a 1 in the location of the topk we allow if we got it right, only 1 of the k for each B can be 1.
            # Important: you can only have 1 right in the k dimension since the label will only have 1 label and our
            if reduction == 'none':
                topk_acc = indicator_which_topk_matched_truth
                assert topk_acc.size() == torch.Size([batch_size])
            elif reduction == 'mean':
                # compute topk accuracies - the model's ability to get it right within it's top k guesses/preds
                topk_acc = indicator_which_topk_matched_truth.mean()  # topk accuracy for entire batch
                assert topk_acc.size() == torch.Size([])
            else:
                raise ValueError(f'Invalid reduction type, got: {reduction=}')
            list_topk_accs.append(topk_acc)
        return tuple(list_topk_accs)  # list of topk accuracies for entire batch [topk1, topk2, ... etc]


def acc_test():
    B = 30
    Dx, Dy = 2, 10
    mdl = nn.Linear(Dx, Dy)
    x = torch.randn(B, Dx)
    y_logits = mdl(x)
    y = torch.randint(high=Dy, size=(B,))
    print(y.size())
    acc1, acc5 = accuracy(output=y_logits, target=y, topk=(1, 5))
    print(f'{acc1=}')
    print(f'{acc5=}')

    accs1, accs5 = accuracy(output=y_logits, target=y, topk=(1, 5), reduction='none')
    print(f'{accs1=}')
    print(f'{accs5=}')
    print(f'{accs1.mean()=}')
    print(f'{accs5.mean()=}')
    print(f'{accs1.std()=}')
    print(f'{accs5.std()=}')
    print(f'{torch_compute_confidence_interval_classification_torch(accs1)=}')
    print(f'{torch_compute_confidence_interval_classification_torch(accs5)=}')


if __name__ == '__main__':
    acc_test()
    print('Done, success! \a')

output

Connected to pydev debugger (build 213.5744.248)
torch.Size([30])
acc1=tensor(0.2000)
acc5=tensor(0.5667)
accs1=tensor([0., 0., 1., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 1., 0., 0., 0., 1., 1., 1., 0., 0., 0., 0.])
accs5=tensor([0., 1., 1., 0., 1., 0., 0., 0., 1., 0., 1., 0., 1., 1., 0., 0., 1., 1.,
        1., 1., 0., 0., 1., 1., 1., 1., 1., 0., 0., 1.])
accs1.mean()=tensor(0.2000)
accs5.mean()=tensor(0.5667)
accs1.std()=tensor(0.4068)
accs5.std()=tensor(0.5040)
torch_compute_confidence_interval_classification_torch(accs1)=tensor(0.1431)
torch_compute_confidence_interval_classification_torch(accs5)=tensor(0.1773)
Done, success! 

how should I be computing confidence intervals (ci)? is the 0-1 loss a special problem with 95 confidence invervals? Since I assume ci tries to approximate a normal but the vector is a bunch of 0's and 1's which seems really ackward to me to approximate as a gaussian...what is going on and what is the right way to compute confidence intervals on the validation set for any loss function?

---

related:

• https://stackoverflow.com/questions/70356922/what-is-the-proper-way-to-compute-95-confidence-intervals-with-pytorch-for-clas
• What is the proper way to report confidence intervals (or std, variance) for machine learning - especially wrt the loss and the accuracy?
• https://www.reddit.com/r/learnmachinelearning/comments/rh5h2f/what_is_the_proper_way_to_report_confidence/
• https://www.quora.com/unanswered/What-is-the-proper-way-to-report-confidence-intervals-or-std-variance-for-machine-learning-especially-with-regards-to-the-loss-and-the-accuracy
• MIT course on CI: https://www.youtube.com/watch?v=MzvRQFYUEFU&list=PLUl4u3cNGP60hI9ATjSFgLZpbNJ7myAg6&index=205

Answer 1

tldr;

Confidence intervals (ci) compute:

• the probability that the true mean is in the given interval (usually written mu_n +- ci

Assumptions:

• traditional confidence intervals statements only hold for statements about the value (parameter, random quantitiy, etc) we want to estimate being the mean
• you have enough samples so that the analysis holds (e.g. the mean $mu_n = 1/n sum_i x_i$, where n>=30 is recommended)

If those assumptions hold (**i.e. your esitmating the true mean via the sample mean with a +- value **) then use the code bellow that I provided called torch_compute_confidence_interval for regression, classification, anything you want.

---

First, asfaik confidence intervals (ci) is an open research problem in deep learning (DL) - so more sophisticated answers probably exist. But I will provide a practical answer that I plan to use (and see others using when reporting results in DL).

To compute confidence intervals we have to understand a little bit of ci first. What they are is a probabilistic statement over the random surveys/samples of data sets that the mean you are trying to report is withing the reported interval. So when people say:

mean_error +- CI for p=95%

it means if you sampled 95 data sets you'd expect the true mean to lie in that interval 95 of the time (but you wouldn't know which ones, so you can't say for any specific interval you compute that the mean will be there).

This means you can only use it for reporting means. This is because the maths that goes behind it (which isn't very hard) approximates the computation of the probability that the bound holds (or the confidence interval holds) by taking advantage that we can compute probabilities analytically for sample means because the approximate a normal according to the central limit theorem CLT. So the specific CI that is computed assumes the quanity you want to compute is a sample mean and computes your +- numbers using this normal approximation. Thus, usually it's recomended to have n>=30 data points for the specific data set you are using but things can still work out nicely since ci can be computed with a t distribution instead of a normal (denoted z in stats software).

Given those assumptions you can simply do the following:

def torch_compute_confidence_interval(data: Tensor,
                                           confidence: float = 0.95
                                           ) -> Tensor:
    """
    Computes the confidence interval for a given survey of a data set.
    """
    n = len(data)
    mean: Tensor = data.mean()
    # se: Tensor = scipy.stats.sem(data)  # compute standard error
    # se, mean: Tensor = torch.std_mean(data, unbiased=True)  # compute standard error
    se: Tensor = data.std(unbiased=True) / (n**0.5)
    t_p: float = float(scipy.stats.t.ppf((1 + confidence) / 2., n - 1))
    ci = t_p * se
    return mean, ci

I've tested it and compared it to things specialized for classification and they agree in values up to 1e-2 so the code works. Output:

Connected to pydev debugger (build 213.5744.248)
x_bernoulli.std()=tensor(0.5040)
ci_95=0.1881992999915952
ci_95_cls=tensor(0.1850)
ci_95_anything=tensor(0.1882)
x_bernoulli.std()=tensor(0.5085, grad_fn=<StdBackward>)
ci_95_torch=tensor(0.1867, grad_fn=<MulBackward0>)
x.std()=tensor(0.9263)
ci_95=0.3458867459004733
ci_95_torch=tensor(0.3459)
x.std()=tensor(1.0181, grad_fn=<StdBackward>)
ci_95_torch=tensor(0.3802, grad_fn=<MulBackward0>)

For more details see my ultimate-utils library where I comment on the maths in the docs: https://github.com/brando90/ultimate-utils/blob/e81a8c3c4425b33e00b3ade172705f20b626b2b1/ultimate-utils-proj-src/uutils/torch_uu/metrics/confidence_intervals.py#L1

---

Comments on DL

If you are reporting the error of a specific model e.g. neural net, like this you are more or less reporting that the true mean error for that very specific neural net and weights lies in those bounds. But as I said this is an open research area so fancier things must be available e.g. consider some layers are actually random, etc.

Answer 2

Well, as I see you have something like 2 populations, where [ones/whole_population] 6/30 & 17/30 are computed as weights, being weighted average... but you didn't prove the Normality of your distribution - with any test - { moreover, your experimental data more suits binomial distribution - having either 0 or 1 output, not Normal!!}... - & as so as C.I.=mean+/-3*sqrt(variance) IS VALID only for Normal distributed variables -- you're getting strange results