TH2 TEfficiency example

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?

https://root.cern.ch/doc/v614/classTH2.html → 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)

1 Like

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:

  1. unrelated to the actual issue, and not really wrong here but I noticed
  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)

  1. That TH1::Divide( … “B”) and TH2::Divide( … “B”) compute different
    errors screams bug to me - even more, the values you report seem plain
    wrong. So if the policy is that “B” shouldn’t be used, then it shouldn’t
    exist (or report an error or use TEfficiency internally). For all I know “B”
    is very popular and often recommended.

Cheers,
Paul

Hi Paul, thank you.

  1. I have fixed the typo in the loop over bins. For completeness, you can find the updated macro below.

  2. 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.