Hi,
I am trying to obtain 2D efficiencies based on two TH2 histograms (h_pass, h_total) with non-equal bins. It would be great if you can provide a working example, in which the efficiencies are drawn using TEfficiency (or another tool). So far I have used TH2::Divide, but the errors are wrong.
Thanks and best regards,
Jonatan
ROOT Version: 6.08/07
Platform: linuxx8664gcc
Did you try the binomial option to get the errors right when using Divide?
ROOT: TH2 Class Reference → Public Member Functions inherited from TH1 → Divide (const TH1 *h1, const TH1 *h2, Double_t c1=1, Double_t c2=1, Option_t *option=“”)
If errors are defined (see TH1::Sumw2), errors are also recalculated Note that if h1 or h2 have Sumw2 set, Sumw2 is automatically called for this if not already set. The resulting errors are calculated assuming uncorrelated histograms. However, if option =“B” is specified, Binomial errors are computed. In this case c1 and c2 do not make real sense and they are ignored.
Hi,
In the example below (th2EfficiencyTest.C) you can see a comparison between the TH2 and TH1 divisions when binomial errors are computed. This simple macro prints the pass, total and efficiency values for both TH1 and TH2, together with their errors in brackets. Clearly the TH2 binomial errors are not properly computed.
eff th2 0.91(0.96) th1 0.91(0.02)
eff th2 0.82(0.90) th1 0.82(0.05)
eff th2 1.00(1.00) th1 1.00(0.00)
eff th2 0.73(0.86) th1 0.73(0.02)
eff th2 0.71(0.84) th1 0.71(0.03)
eff th2 0.76(0.87) th1 0.76(0.02)
Thanks,
Jonatan
https://piedra.web.cern.ch/piedra/root/th2_efficiency.root
https://piedra.web.cern.ch/piedra/root/th2EfficiencyTest.C
Hi,
You should use TEfficiency if you want correct error (i.e. not using the normal approximation).
Attached here is an example of using TEfficiency with 2D histogram
Lorenzo
example_TEfficiency2D.C (1.0 KB)
Thank you Lorenzo.
I have run your example and expanded it. The version below compares the TH1 results with the TEfficiency output. I find more reasonable the TH1 errors. Besides, the efficiency values are not exactly equal.
Best,
Jonatan
http://piedra.web.cern.ch/piedra/root/example_TEfficiency2D.C
Hi,
two things:
for (Int_t i=1; i<=h_eff->GetNbinsX(); i++) {
for (Int_t j=1; j<=h_eff->GetNbinsX(); j++) {
the second here should probably be Y instead of X.
(In your example GetNbins* returns 14 in any case, so all fine for the sake
of what we’re discussing)
Cheers,
Paul
Hi Paul, thank you.
I have fixed the typo in the loop over bins. For completeness, you can find the updated macro below.
With these tests that I have performed I not only get different errors between TH1::Divide and TH2::Divide when “B” is used; I also get different errors (and even efficiencies) between TH1::Divide and TEfficiency.
I can understand that “B” is popular for TH1. But from these findings, either I’m making some mistake, or the TH2 (efficiency) errors have a bug somewhere.
All the best,
Jonatan
https://piedra.web.cern.ch/piedra/root/example_TEfficiency2D.C
Hi Jonatan,
Since this puzzled me (TH2::Divide isn’t implemented, it just calls exactly the
same code as TH1::Divide) I tried around and noticed that calling Sumw2
manually on h_pass is enough to get consistent errors:
TH2F* h_pass = (TH2F*)file->Get("h_pass");
TH2F* h_total = (TH2F*)file->Get("h_total");
h_pass ->Sumw2();
I’m not 100% sure but I suspect your calls to
TH1::SetDefaultSumw2();
TH2::SetDefaultSumw2();
are semi-futile: (I think) they don’t affect already existing histograms (such
as those coming from your .root file) but affect only newly created histograms.
Therefore the TH1 versions work (you create these from scratch and they get
created with sumw2) but the TH2 versions don’t
((*h_pass->GetSumw2()).GetSize() returns 0).
Cheers,
Paul
PS: Let me be clear. This fix shows that technically Divide(…“b”) does the
right thing for 2d histograms (or at least the same as TH1 [discussion of the
binomial divide vs TEfficiency goes here]). But this was too easy to run into
and too hard to spot. Looking at the Divide code, the binomial-or-not logic
isn’t reached when Sumw2 isn’t set.
if (fSumw2.fN) {
if (b2 == 0) { fSumw2.fArray[i] = 0; continue; }
Double_t b1sq = b1 * b1; Double_t b2sq = b2 * b2;
Double_t c1sq = c1 * c1; Double_t c2sq = c2 * c2;
Double_t e1sq = h1->GetBinErrorSqUnchecked(i);
Double_t e2sq = h2->GetBinErrorSqUnchecked(i);
if (binomial) {
if (b1 != b2) {
// in the case of binomial statistics c1 and c2 must be 1 otherwise it does not make sense
// c1 and c2 are ignored
//fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));//this is the formula in Hbook/Hoper1
//fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/b2); // old formula from G. Flucke
// formula which works also for weighted histogram (see http://root-forum.cern.ch/viewtopic.php?t=3753 )
fSumw2.fArray[i] = TMath::Abs( ( (1. - 2.* b1 / b2) * e1sq + b1sq * e2sq / b2sq ) / b2sq );
} else {
//in case b1=b2 error is zero
//use TGraphAsymmErrors::BayesDivide for getting the asymmetric error not equal to zero
fSumw2.fArray[i] = 0;
}
} else {
fSumw2.fArray[i] = c1sq * c2sq * (e1sq * b2sq + e2sq * b1sq) / (c2sq * c2sq * b2sq * b2sq);
}
}
And here one could issue a warning or error because if you provide the “b”
option, then clearly you are interested in the binomial error computation
and it should not be skipped.
Thanks a lot Paul. You are right, my call to
TH2::SetDefaultSumw2();
was useless. I have modified the macro [1] and I get exactly the same, as it should be, for both TH1 and TH2. The puzzle remains for the second example [2]. When I run it I get
bin 1,1 | pass th2 121(11) th1 121(11) | total th2 131(11) th1 131(11) | eff th2 0.917 +- 0.958 th1 0.924 +- 0.023
bin 1,2 | pass th2 28( 5) th1 28( 5) | total th2 29( 5) th1 29( 5) | eff th2 0.935 +- 0.967 th1 0.966 +- 0.034
bin 2,1 | pass th2 194(14) th1 194(14) | total th2 221(15) th1 221(15) | eff th2 0.874 +- 0.935 th1 0.878 +- 0.022
bin 2,2 | pass th2 72( 8) th1 72( 8) | total th2 77( 9) th1 77( 9) | eff th2 0.924 +- 0.961 th1 0.935 +- 0.028
bin 3,1 | pass th2 254(16) th1 254(16) | total th2 290(17) th1 290(17) | eff th2 0.873 +- 0.934 th1 0.876 +- 0.019
bin 3,2 | pass th2 78( 9) th1 78( 9) | total th2 89( 9) th1 89( 9) | eff th2 0.868 +- 0.932 th1 0.876 +- 0.035
bin 4,1 | pass th2 105(10) th1 105(10) | total th2 119(11) th1 119(11) | eff th2 0.876 +- 0.936 th1 0.882 +- 0.030
bin 4,2 | pass th2 34( 6) th1 34( 6) | total th2 36( 6) th1 36( 6) | eff th2 0.921 +- 0.960 th1 0.944 +- 0.038
This means that both the efficiency and the error differ between TEfficiency and TH1. And clearly, the right efficiency value is given by TH1. In addition, when I do
TH2::SetDefaultSumw2();
the TEfficiency errors are all of them zero.
Cheers,
Jonatan
[1] https://piedra.web.cern.ch/piedra/root/th2EfficiencyTest.C
[2] https://piedra.web.cern.ch/piedra/root/example_TEfficiency2D.C
[3] https://piedra.web.cern.ch/piedra/root/th2_efficiency.root
Hi Jonatan,
thanks for the feedback.
I changed the printout a bit:
printf(" bin %d,%d | teff %.3f - %.3f + %.3f th from teff %.3f +- %.3f (sanity check sqrt(N) is: %.3f)\n",
i, j,
pEff ->GetEfficiency(heff->GetBin(i,j)),
pEff ->GetEfficiencyErrorLow(heff->GetBin(i,j)),
pEff ->GetEfficiencyErrorUp(heff->GetBin(i,j)),
heff ->GetBinContent(i,j), heff ->GetBinError(i,j),
TMath::Sqrt(heff ->GetBinContent(i,j))
);
To compare the errors computed by TEfficiency to the errors in the TH2
extracted from TEfficiency. I think the values speak for themselves.
bin 1,1 | teff 0.917 - 0.070 + 0.068 th from teff 0.917 +- 0.957 (sanity check sqrt(N) is: 0.957)
bin 1,2 | teff 0.833 - 0.140 + 0.132 th from teff 0.833 +- 0.913 (sanity check sqrt(N) is: 0.913)
bin 2,1 | teff 0.848 - 0.061 + 0.061 th from teff 0.848 +- 0.921 (sanity check sqrt(N) is: 0.921)
bin 2,2 | teff 0.818 - 0.111 + 0.109 th from teff 0.818 +- 0.905 (sanity check sqrt(N) is: 0.905)
bin 3,1 | teff 0.846 - 0.069 + 0.069 th from teff 0.846 +- 0.920 (sanity check sqrt(N) is: 0.920)
bin 3,2 | teff 0.889 - 0.094 + 0.090 th from teff 0.889 +- 0.943 (sanity check sqrt(N) is: 0.943)
bin 4,1 | teff 0.846 - 0.095 + 0.094 th from teff 0.846 +- 0.920 (sanity check sqrt(N) is: 0.920)
bin 4,2 | teff 0.800 - 0.168 + 0.157 th from teff 0.800 +- 0.894 (sanity check sqrt(N) is: 0.894)
The exported histogram only contains central values but no errors.
(to be fair, with the bayesian statistics option, the TEfficiency errors are
asymmetric and TH only have symmetric errors, so the conversion is not straight
forward).
For the other things you wrote:
The efficiencies reported by TEfficiency are not passed/total. This is desired and correct.
It’s unfortunately a statistics lecture on its own and often overlooked, but see https://root.cern.ch/doc/master/classTEfficiency.html IV statistics options → Bayesian methods.
I myself liked the read of https://arxiv.org/abs/0908.0130. (I think on wikipedia there’s also a dive into the topic once you go to the conjugate prior of the binomial distribution).
It basically boils down to the fact that the likelihood for a true ε given an
observed k/N has a different functional form than the probabilities for k/N
given a true ε.
From the mathematical side you’ll notice that k/N can only take rational
values, while ε may well be an irrational number.
From the practical side, if ε is perfectly 100% there’s no chance k will ever
be smaller than N. BUT a true Bayesian (with a reasonable prior) would never
assume a coin will always show tails, just because out of ten tosses it always
showed tails. They will expect some value, probably somewhere between 90 and
100 percent.
Furthermore you said “In addition, when I do TH2::SetDefaultSumw2(); the TEfficiency errors are all of them zero.”.
puh, haven’t looked into that yet. Sounds like a bug to me.
Thank you Paul.
I have taken your printout, and surprisingly I get different (and not asymmetric) errors. The macro has been updated in [1]. I’m using ROOT 6.08/07 built for linuxx8664gcc.
bin 1,1 | teff 0.917 - 0.023 + 0.023 th from teff 0.917 +- 0.958 (sanity check sqrt(N) is: 0.958)
bin 1,2 | teff 0.935 - 0.041 + 0.041 th from teff 0.935 +- 0.967 (sanity check sqrt(N) is: 0.967)
bin 2,1 | teff 0.874 - 0.022 + 0.022 th from teff 0.874 +- 0.935 (sanity check sqrt(N) is: 0.935)
bin 2,2 | teff 0.924 - 0.029 + 0.029 th from teff 0.924 +- 0.961 (sanity check sqrt(N) is: 0.961)
bin 3,1 | teff 0.873 - 0.019 + 0.019 th from teff 0.873 +- 0.934 (sanity check sqrt(N) is: 0.934)
bin 3,2 | teff 0.868 - 0.035 + 0.035 th from teff 0.868 +- 0.932 (sanity check sqrt(N) is: 0.932)
bin 4,1 | teff 0.876 - 0.030 + 0.030 th from teff 0.876 +- 0.936 (sanity check sqrt(N) is: 0.936)
bin 4,2 | teff 0.921 - 0.042 + 0.042 th from teff 0.921 +- 0.960 (sanity check sqrt(N) is: 0.960)
[1] https://piedra.web.cern.ch/piedra/root/example_TEfficiency2D.C
Thanks a lot for the statistics explanation. I will read the paper that you pointed me to.
Hi,
ah, sorry I missed that I made more changes which I forgot when writing the previous posts. With the script you just posted I get asymmetric errors when I reduce the statistics:
for(int i=0; i<100; ++i)
{
f2.GetRandom2(x,y);
Maybe a word of caution. The paper and TEfficiency are a bit on the rigorous side.
When reading the paper that the typical binomial efficiency estimator and error formula are k/N and sqrt(ε(1-ε)/N)=sqrt(k(N-k)/N^3). And the Bayesian estimators only increase k, N, and (N-k) by numbers at the order of 1.
So as long as N=N+2, k=k+1, (N-k)=(N-k)+1 are good approximations, the typical binomial errors are okayish.
And if you’re in the regime where these don’t hold … you can argue that your real problem is a small data sample and not the formulae you’re using to analyse it.
Cheers,
Paul
Thank you! Numbers agree perfectly now.
Thanks also for the wise word of caution.
Cheers,
Jonatan
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.