What I don't like about binomial errors and how significant is the performance difference of two BDTs?
A topic I already discussed in my master thesis is the determination of efficiencies, commonly denoted as \(\varepsilon=k/N\). Even assuming absolute certainty on \(k\) and \(N\) (that means: no fit is involved, it’s just counting, the only uncertainty one would assume is the standard \(\sqrt(N)\) when referring to the uncertainty on a Poisson mean) the error propagation seems hard and is often done wrong. The naive propagation (relative error is quadratic sum of relative errors of \(k\) and \(N\)) is clearly wrong because it neglects the correlation of \(k\) and \(N\). Once realising this, one uses the so-called binomial error, \(\sqrt{\varepsilon * (1-\varepsilon)/N}\), which I also dislike.
Where I think, things are done wrong
The formula \(\frac{\varepsilon \cdot(1-\varepsilon)}{N}\) describes (correctly) the variance of the random variable \(\frac{k}{N}\) around its true expectation value \(\varepsilon\):
$$E\left(\left(\frac{k}{N}-E\left(\frac{k}{N}\right)\right)^2\right) = E\left(\left(\frac{k}{N}\right)^2\right)-\left(E\left(\frac{k}{N}\right)\right)^2$$$$=\frac{\varepsilon\cdot(1-\varepsilon)}{N}$$
But, let’s be pedantic, in real life we don’t know \(\varepsilon\), we only know \(k\) and \(N\) and have a good estimator \(\hat{\varepsilon}=\frac{k}{N}\). It’s common practise to (ab)use the estimate for a value in formulas for other estimators (the most famous example of fixing this is the \(\sqrt{N-1}\) factor in the observed standard deviation). But what I’d like to point out is that the variance of \(\frac{k}{N}\) around a true \(\varepsilon\) is the width of the probability \(p\left(\frac{k}{N}|\varepsilon,N\right)\), i.e. for fixed given \(N\) and \(\varepsilon\). But what we actually want to know is the width of the likelihood \(p\left(\varepsilon|\frac{k}{N},N\right)\). I have to admit that sounds ununderstandably pedantic. A more intuitive way of describing it is the following.
The function \(p\left(\frac{k}{N},\varepsilon,N\right)\) has three arguments (and we ignore which are fixed and known and which we want to know). This function has a bump whenever \(\varepsilon=\frac{k}{N}\) and we want to know its width. The standard binomial error measures the width in \(\frac{k}{N}\) direction, but when one wants to report the uncertainty on an efficiency estimated by \(\frac{k}{N}\), then one really wants to know the width in \(\varepsilon\) direction - something completely different.
Friends of type theory or category theory or maybe just programmers…
…will recognize that \(p\) has a real or floating point input \(\varepsilon\) and an integer \(k\). I’m disussing “which real value of \(\varepsilon\) maximises \(p\) and how bad is it if we’re a bit off?” and I’m under the impression the commonly people instead answer “which integer value of \(k\) maximises \(p\) and how bad is it if we’re a bit off?”
Literature
The paper I most often refer to is arXiv:0908.0130, which also told me that the process of interpreting \(p(x,y)\) once as probability \(p(x|y)\) and once as likelihood \(p(y|x)\) is something which literature calls the conjugate prior. For the special case of selection efficiencies the relevant wikipedia pages are those of the binomial and the beta distribution.
Oh no! Bayesian methods!
Yes, I’m not a huge fan of having to decide for priors either. But here things turn out rather nice. The Beta distribution has two parameters, \(\alpha\) and \(\beta\), and they essentially are the two counts for classes of events, selected and rejected: \(\alpha=k+x\) and \(\beta = (N-k)+y\). Where \(x\) and \(y\) are parameters of the prior distribution (here some nice mathematical property of combinations of Beta functions jumps in) and due to symmetry considerations (swap selected and rejected) one wants \(x=y\). Moreover, the maximum likelihood estimate for the efficiency becomes:
$$\hat{\varepsilon}_{ml}=\frac{\alpha-1}{\alpha+\beta-2}$$Which is the standard formula for \(x=y=1\). Now Casadei explains that one can also use Jeffrey’s prior - sounds quite convincing, then \(x=y=1/2\), but after all one only adds a little bit in numerator and denominator, let’s hope half an event more or less doesn’t matter. The mean of the likelihood is not coinciding with the mode (maximum likelihood value):
$$E(\varepsilon)=\frac{\alpha}{\alpha+\beta}$$Result
Luckily Wikipedia knows the variance of the Beta distribution because I wouldn’t want to solve the integral myself:
$$\sigma_\varepsilon = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$$Once you use the priors \(x=y=1\), the denominator \(\alpha+\beta\) becomes \(N\) and one can factorise:
$$\sigma^2_\varepsilon = \frac{(k+1)(N-k+1)}{(N+1)^2(N+2)} = \left(\frac{k+1}{N+1}\right)\left(\frac{N-k+1}{N+1}\right) \frac{1}{N+2}$$Almost the binomial error (up to all the \(+1\) on both sides of the fractions which approximate \(\varepsilon\) and \((1-\varepsilon)\) and the \(+2\) in the last term.
The important difference
In contrast to the binomial error, this one does not have pathologic behaviour in the extreme cases \(\varepsilon=0\) or \(1\)! While a process of exactly \(100\,\%\) efficiency will accept each and every event (and thus \(k=N\)), one must not infer \(100\,\%\) efficiency from \(k=N\). (Assume \(N=1\) and this single event gets accepted. \(100\,\%\) should still be your best guess, but not with vanishing uncertainty.)
More serious stuff: efficiency differences
Now a more serious problem is this: I have two selections and want to know if one is really better than the other. This happens (to me) in machine learning a lot in the sense that I end up with two BDTs, optimise a cut value on their responses, measure the efficiency on some left-over = unbiassed MC and see: “There’s not much MC statistics. Just drawing the ROC curve (with error bars from the above formula) It doesn’t seem a significant difference. But hey, the errors on the two curves are not independent! I measured both on the same sample!”
Reading on on Wikipedia, I realised this can be modelled with a multinomial process with four categories:
- event rejected by both BDTs
- event accepted by BDT1 and rejected by BDT2
- event accepted by BDT2 and rejected by BDT1
- event accepted by both BDTs
The corresponding probabilities \(p_i (1\leq i\leq 4)\) are not the physically interesting one. Instead one would call:
$$\varepsilon_1 = p_2+p_4$$$$\varepsilon_2 = p_3+p_4$$
$$\varepsilon_{1\wedge 2} = p_4$$
$$\overline{\varepsilon} = p_4 + \frac{(p_2+p_3)}{2}$$
$$\Delta(\varepsilon) = p_2-p_3$$
Wikipedia helps again with the conjugate prior distribution and its variance and covariance. For event counts \(k_i\) in each category:
$$\mathrm{Var}(p_i)=\frac{\alpha_i(-\alpha_i+\sum_j\alpha_j)}{(\sum_j\alpha_j)^2(1+\sum_j\alpha_j)}$$$$\mathrm{Cov}(p_i,p_j)=\frac{-\alpha_i\alpha_j}{(\sum_m\alpha_m)^2(1+\sum_m\alpha_m)}$$
Again \(\alpha_i=k_i+x\) (and \(x\) can reasonably be \(1\) or \(\frac{1}{2}\)). Error propagation from here on (considering the covariance term) is straight forward and we can tell if an efficiency difference is significant!
Correlation
I left out the correlation between the two classifiers. To make the formulas
more like proper math, define a random variable \(X_1\) which takes the value 1
(selected) or 0 (rejected). Fun fact: It’s a Bernoulli distributed variable
(corner case of Binomial). Similarly may \(X_2\) be the same for BDT2 and
\(X_{12}\) for the and of both BDTs. (Here nicely integer arithmetic and bool
logic meet because \(X_{12}=X_1\cdot X_2\))
The efficiencies are the same as usual:
$$\varepsilon_1 = E(X_1)$$$$\varepsilon_2 = E(X_2)$$
$$\varepsilon_{12} = E(X_{12})$$
The variances:
$$\mathrm{Var}(X_1) = \varepsilon_1\cdot(1-\varepsilon_1)$$$$\mathrm{Var}(X_2) = \varepsilon_2\cdot(1-\varepsilon_2)$$
The Covariance:
$$\mathrm{Cov}(X_1,X_2) = E(X_1X_2)-E(X_1)E(X_2)$$$$ = E(X_{12})-E(X_1)E(X_2)$$
$$ = \varepsilon_{12}-\varepsilon_1\varepsilon_2$$
As desired, the covariance is zero for uncorrelated classifiers. The Correlation is consequently
$$\mathrm{Cor}(X_1,X_2) = \frac{\varepsilon_{12}-\varepsilon_1\varepsilon_2}{ \sqrt{ \varepsilon_1\cdot(1-\varepsilon_1)\varepsilon_2\cdot(1-\varepsilon_2)}}$$Not utterly helpful though, since it’s really the correlation of the classifiers and not the correlation of their efficiency estimates.
Confidence intervals vs. Variance and Software!
One can not only determine variances but also confidence intervals (central or
shortests; or those corresponding to likelihood values; ;r from Feldman-Cousins
constructions). Papers exist (at least I’ve seen one for the shortest Interval.
Unfortunately not on arXiv, so no link.) Then again you don’t need to derive
all of that, as it exists in TEfficiency (use the static functions!). E.g.:
1auto low = TEfficiency::Bayesian(N,k,alpha,beta,0.683,false,true);
Also, in RooFit RooEfficiency is helpful to model an efficiency.
Coverage?
I didn’t mention coverages, it’s lenghtly discussed in the TEfficiency
documentation. The basic problem is that discrete variables don’t allow perfect
coverage - one will usually overcover.