ROOT logo
// Perform a fit to a set of data with binomial errors 
// like those derived from the division of two histograms. 
// Three different fits are performed and compared:  
//
//   -  simple least square fit to the divided histogram obtained 
//      from TH1::Divide with option b
//   -  least square fit to the TGraphAsymmErrors obtained from 
//      TGraphAsymmErrors::BayesDivide
//   -  likelihood fit performed on the dividing histograms using 
//      binomial statistics with the TBinomialEfficiency class
// 
// The first two methods are biased while the last one  is statistical correct. 
// Running the script passing an integer value n larger than 1, n fits are 
// performed and the bias are also shown.   
// To run the script : 
// 
//  to show the bias performing 100 fits for 1000 events per "experiment"
//  root[0]: .x TestBinomial.C+
// 
//  to show the bias performing 100 fits for 1000 events per "experiment"
//           .x TestBinomial.C+(100, 1000)
// 
//
#include "TBinomialEfficiencyFitter.h"
#include "TVirtualFitter.h"
#include "TH1.h"
#include "TRandom3.h"
#include "TF1.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include <cassert>
#include <iostream>

void TestBinomial(int nloop = 100, int nevts = 100, bool plot=false)
{
   gStyle->SetMarkerStyle(20);
   gStyle->SetLineWidth(2.0);
   gStyle->SetOptStat(11);

   TObjArray hbiasNorm;
   hbiasNorm.Add(new TH1D("h0Norm", "Bias Histogram fit",100,-5,5));
   hbiasNorm.Add(new TH1D("h1Norm","Bias Binomial fit",100,-5,5));
   TObjArray hbiasThreshold;
   hbiasThreshold.Add(new TH1D("h0Threshold", "Bias Histogram fit",100,-5,5));
   hbiasThreshold.Add(new TH1D("h1Threshold","Bias Binomial fit",100,-5,5));
   TObjArray hbiasWidth;
   hbiasWidth.Add(new TH1D("h0Width", "Bias Histogram fit",100,-5,5));
   hbiasWidth.Add(new TH1D("h1Width","Bias Binomial fit",100,-5,5));
   TH1D* hChisquared = new TH1D("hChisquared", 
      "#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);

   //TVirtualFitter::SetDefaultFitter("Minuit2");

   // Note: in order to be able to use TH1::FillRandom() to generate
   //       pseudo-experiments, we use a trick: generate "selected"
   //       and "non-selected" samples independently. These are
   //       statistically independent and therefore can be safely
   //       added to yield the "before selection" sample.


   // Define (arbitrarily?) a distribution of input events.
   // Here: assume a x^(-2) distribution. Boundaries: [10, 100].

   Double_t xmin =10, xmax = 100;
   TH1D* hM2D = new TH1D("hM2D", "x^(-2) denominator distribution", 
      45, xmin, xmax);
   TH1D* hM2N = new TH1D("hM2N", "x^(-2) numerator distribution",   
      45, xmin, xmax);
   TH1D* hM2E = new TH1D("hM2E", "x^(-2) efficiency",               
      45, xmin, xmax);

   TF1*  fM2D = new TF1("fM2D", "(1-[0]/(1+exp(([1]-x)/[2])))/(x*x)", 
      xmin, xmax);
   TF1*  fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",     
      xmin, xmax);
   TF1*  fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",       
      xmin, xmax);
   TF1*  fM2Fit2 = 0;

   TRandom3 rb(0);

   // First try: use a single set of parameters.
   // For each try, we need to find the overall normalization

   Double_t norm = 0.80;
   Double_t threshold = 25.0;
   Double_t width = 5.0;

   fM2D->SetParameter(0, norm);
   fM2D->SetParameter(1, threshold);
   fM2D->SetParameter(2, width);
   fM2N->SetParameter(0, norm);
   fM2N->SetParameter(1, threshold);
   fM2N->SetParameter(2, width);
   Double_t integralN = fM2N->Integral(xmin, xmax);
   Double_t integralD = fM2D->Integral(xmin, xmax);
   Double_t fracN = integralN/(integralN+integralD);
   Int_t nevtsN = rb.Binomial(nevts, fracN);
   Int_t nevtsD = nevts - nevtsN;

   gStyle->SetOptFit(1111);

   // generate many times to see the bias
   for (int iloop = 0; iloop < nloop; ++iloop) {

     // generate pseudo-experiments
     hM2D->Reset();
     hM2N->Reset();
     hM2D->FillRandom(fM2D->GetName(), nevtsD);
     hM2N->FillRandom(fM2N->GetName(), nevtsN);
     hM2D->Add(hM2N);

     // construct the "efficiency" histogram
     hM2N->Sumw2();
     hM2E->Divide(hM2N, hM2D, 1, 1, "b");

     // Fit twice, using the same fit function.
     // In the first (standard) fit, initialize to (arbitrary) values.
     // In the second fit, use the results from the first fit (this
     // makes it easier for the fit -- but the purpose here is not to
     // show how easy or difficult it is to obtain results, but whether
     // the CORRECT results are obtained or not!).

     fM2Fit->SetParameter(0, 0.5);
     fM2Fit->SetParameter(1, 15.0);
     fM2Fit->SetParameter(2, 2.0);
     fM2Fit->SetParError(0, 0.1);
     fM2Fit->SetParError(1, 1.0);
     fM2Fit->SetParError(2, 0.2);
     TCanvas* cEvt;
     if (plot) {
       cEvt = new TCanvas(Form("cEnv%d",iloop),
			  Form("plots for experiment %d", iloop),
			  1000, 600);
       cEvt->Divide(1,2);
       cEvt->cd(1);
       hM2D->DrawCopy("HIST");
       hM2N->SetLineColor(kRed);
       hM2N->DrawCopy("HIST SAME");
       cEvt->cd(2);
     }
     for (int fit = 0; fit < 2; ++fit) {
       Int_t status = 0;
       switch (fit) {
       case 0:
	 hM2E->Fit(fM2Fit, "RN");
	 if (plot) {
	   hM2E->DrawCopy("E");
	   fM2Fit->DrawCopy("SAME");
	 }
	 break;
       case 1:
	 if (fM2Fit2) delete fM2Fit2;
	 fM2Fit2 = dynamic_cast<TF1*>(fM2Fit->Clone("fM2Fit2"));
	 if (fM2Fit2->GetParameter(0) >= 1.0)
	   fM2Fit2->SetParameter(0, 0.95);
	 fM2Fit2->SetParLimits(0, 0.0, 1.0);
	 TBinomialEfficiencyFitter bef(hM2N, hM2D);
	 status = bef.Fit(fM2Fit2,"RI");
	 if (status!=0) {
	   std::cerr << "Error performing binomial efficiency fit, result = "
		     << status << std::endl;
	   continue;
	 }
	 if (plot) {
	   fM2Fit2->SetLineColor(kRed);
	   fM2Fit2->DrawCopy("SAME");
	 }
       }

       if (status != 0) break;

       Double_t fnorm = fM2Fit->GetParameter(0);
       Double_t enorm = fM2Fit->GetParError(0);
       Double_t fthreshold = fM2Fit->GetParameter(1);
       Double_t ethreshold = fM2Fit->GetParError(1);
       Double_t fwidth = fM2Fit->GetParameter(2);
       Double_t ewidth = fM2Fit->GetParError(2);
       if (fit == 1) {
	 fnorm = fM2Fit2->GetParameter(0);
	 enorm = fM2Fit2->GetParError(0);
	 fthreshold = fM2Fit2->GetParameter(1);
	 ethreshold = fM2Fit2->GetParError(1);
	 fwidth = fM2Fit2->GetParameter(2);
	 ewidth = fM2Fit2->GetParError(2);
	 hChisquared->Fill(fM2Fit2->GetProb());
       }

       TH1D* h = dynamic_cast<TH1D*>(hbiasNorm[fit]);
       h->Fill((fnorm-norm)/enorm);
       h = dynamic_cast<TH1D*>(hbiasThreshold[fit]);
       h->Fill((fthreshold-threshold)/ethreshold);
       h = dynamic_cast<TH1D*>(hbiasWidth[fit]);
       h->Fill((fwidth-width)/ewidth);
     }
   }
   

   TCanvas* c1 = new TCanvas("c1",
      "Efficiency fit biases",10,10,1000,800); 
   c1->Divide(2,2);

   TH1D *h0, *h1;
   c1->cd(1);
   h0 = dynamic_cast<TH1D*>(hbiasNorm[0]);
   h0->Draw("HIST");
   h1 = dynamic_cast<TH1D*>(hbiasNorm[1]);
   h1->SetLineColor(kRed);
   h1->Draw("HIST SAMES");
   TLegend* l1 = new TLegend(0.1, 0.75, 0.5, 0.9, 
      "plateau parameter", "ndc");
   l1->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
      %4.2f", h0->GetMean(), h0->GetRMS()), "l");
   l1->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
      %4.2f", h1->GetMean(), h1->GetRMS()), "l");
   l1->Draw();

   c1->cd(2);
   h0 = dynamic_cast<TH1D*>(hbiasThreshold[0]);
   h0->Draw("HIST");
   h1 = dynamic_cast<TH1D*>(hbiasThreshold[1]);
   h1->SetLineColor(kRed);
   h1->Draw("HIST SAMES");
   TLegend* l2 = new TLegend(0.1, 0.75, 0.5, 0.9, 
      "threshold parameter", "ndc");
   l2->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
      %4.2f", h0->GetMean(), h0->GetRMS()), "l");
   l2->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
      %4.2f", h1->GetMean(), h1->GetRMS()), "l");
   l2->Draw();

   c1->cd(3);
   h0 = dynamic_cast<TH1D*>(hbiasWidth[0]);
   h0->Draw("HIST");
   h1 = dynamic_cast<TH1D*>(hbiasWidth[1]);
   h1->SetLineColor(kRed);
   h1->Draw("HIST SAMES");
   TLegend* l3 = new TLegend(0.1, 0.75, 0.5, 0.9, "width parameter", "ndc");
   l3->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
      %4.2f", h0->GetMean(), h0->GetRMS()), "l");
   l3->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
      %4.2f", h1->GetMean(), h1->GetRMS()), "l");
   l3->Draw();

   c1->cd(4);
   hChisquared->Draw("HIST");
}
 TestBinomial.C:1
 TestBinomial.C:2
 TestBinomial.C:3
 TestBinomial.C:4
 TestBinomial.C:5
 TestBinomial.C:6
 TestBinomial.C:7
 TestBinomial.C:8
 TestBinomial.C:9
 TestBinomial.C:10
 TestBinomial.C:11
 TestBinomial.C:12
 TestBinomial.C:13
 TestBinomial.C:14
 TestBinomial.C:15
 TestBinomial.C:16
 TestBinomial.C:17
 TestBinomial.C:18
 TestBinomial.C:19
 TestBinomial.C:20
 TestBinomial.C:21
 TestBinomial.C:22
 TestBinomial.C:23
 TestBinomial.C:24
 TestBinomial.C:25
 TestBinomial.C:26
 TestBinomial.C:27
 TestBinomial.C:28
 TestBinomial.C:29
 TestBinomial.C:30
 TestBinomial.C:31
 TestBinomial.C:32
 TestBinomial.C:33
 TestBinomial.C:34
 TestBinomial.C:35
 TestBinomial.C:36
 TestBinomial.C:37
 TestBinomial.C:38
 TestBinomial.C:39
 TestBinomial.C:40
 TestBinomial.C:41
 TestBinomial.C:42
 TestBinomial.C:43
 TestBinomial.C:44
 TestBinomial.C:45
 TestBinomial.C:46
 TestBinomial.C:47
 TestBinomial.C:48
 TestBinomial.C:49
 TestBinomial.C:50
 TestBinomial.C:51
 TestBinomial.C:52
 TestBinomial.C:53
 TestBinomial.C:54
 TestBinomial.C:55
 TestBinomial.C:56
 TestBinomial.C:57
 TestBinomial.C:58
 TestBinomial.C:59
 TestBinomial.C:60
 TestBinomial.C:61
 TestBinomial.C:62
 TestBinomial.C:63
 TestBinomial.C:64
 TestBinomial.C:65
 TestBinomial.C:66
 TestBinomial.C:67
 TestBinomial.C:68
 TestBinomial.C:69
 TestBinomial.C:70
 TestBinomial.C:71
 TestBinomial.C:72
 TestBinomial.C:73
 TestBinomial.C:74
 TestBinomial.C:75
 TestBinomial.C:76
 TestBinomial.C:77
 TestBinomial.C:78
 TestBinomial.C:79
 TestBinomial.C:80
 TestBinomial.C:81
 TestBinomial.C:82
 TestBinomial.C:83
 TestBinomial.C:84
 TestBinomial.C:85
 TestBinomial.C:86
 TestBinomial.C:87
 TestBinomial.C:88
 TestBinomial.C:89
 TestBinomial.C:90
 TestBinomial.C:91
 TestBinomial.C:92
 TestBinomial.C:93
 TestBinomial.C:94
 TestBinomial.C:95
 TestBinomial.C:96
 TestBinomial.C:97
 TestBinomial.C:98
 TestBinomial.C:99
 TestBinomial.C:100
 TestBinomial.C:101
 TestBinomial.C:102
 TestBinomial.C:103
 TestBinomial.C:104
 TestBinomial.C:105
 TestBinomial.C:106
 TestBinomial.C:107
 TestBinomial.C:108
 TestBinomial.C:109
 TestBinomial.C:110
 TestBinomial.C:111
 TestBinomial.C:112
 TestBinomial.C:113
 TestBinomial.C:114
 TestBinomial.C:115
 TestBinomial.C:116
 TestBinomial.C:117
 TestBinomial.C:118
 TestBinomial.C:119
 TestBinomial.C:120
 TestBinomial.C:121
 TestBinomial.C:122
 TestBinomial.C:123
 TestBinomial.C:124
 TestBinomial.C:125
 TestBinomial.C:126
 TestBinomial.C:127
 TestBinomial.C:128
 TestBinomial.C:129
 TestBinomial.C:130
 TestBinomial.C:131
 TestBinomial.C:132
 TestBinomial.C:133
 TestBinomial.C:134
 TestBinomial.C:135
 TestBinomial.C:136
 TestBinomial.C:137
 TestBinomial.C:138
 TestBinomial.C:139
 TestBinomial.C:140
 TestBinomial.C:141
 TestBinomial.C:142
 TestBinomial.C:143
 TestBinomial.C:144
 TestBinomial.C:145
 TestBinomial.C:146
 TestBinomial.C:147
 TestBinomial.C:148
 TestBinomial.C:149
 TestBinomial.C:150
 TestBinomial.C:151
 TestBinomial.C:152
 TestBinomial.C:153
 TestBinomial.C:154
 TestBinomial.C:155
 TestBinomial.C:156
 TestBinomial.C:157
 TestBinomial.C:158
 TestBinomial.C:159
 TestBinomial.C:160
 TestBinomial.C:161
 TestBinomial.C:162
 TestBinomial.C:163
 TestBinomial.C:164
 TestBinomial.C:165
 TestBinomial.C:166
 TestBinomial.C:167
 TestBinomial.C:168
 TestBinomial.C:169
 TestBinomial.C:170
 TestBinomial.C:171
 TestBinomial.C:172
 TestBinomial.C:173
 TestBinomial.C:174
 TestBinomial.C:175
 TestBinomial.C:176
 TestBinomial.C:177
 TestBinomial.C:178
 TestBinomial.C:179
 TestBinomial.C:180
 TestBinomial.C:181
 TestBinomial.C:182
 TestBinomial.C:183
 TestBinomial.C:184
 TestBinomial.C:185
 TestBinomial.C:186
 TestBinomial.C:187
 TestBinomial.C:188
 TestBinomial.C:189
 TestBinomial.C:190
 TestBinomial.C:191
 TestBinomial.C:192
 TestBinomial.C:193
 TestBinomial.C:194
 TestBinomial.C:195
 TestBinomial.C:196
 TestBinomial.C:197
 TestBinomial.C:198
 TestBinomial.C:199
 TestBinomial.C:200
 TestBinomial.C:201
 TestBinomial.C:202
 TestBinomial.C:203
 TestBinomial.C:204
 TestBinomial.C:205
 TestBinomial.C:206
 TestBinomial.C:207
 TestBinomial.C:208
 TestBinomial.C:209
 TestBinomial.C:210
 TestBinomial.C:211
 TestBinomial.C:212
 TestBinomial.C:213
 TestBinomial.C:214
 TestBinomial.C:215
 TestBinomial.C:216
 TestBinomial.C:217
 TestBinomial.C:218
 TestBinomial.C:219
 TestBinomial.C:220
 TestBinomial.C:221
 TestBinomial.C:222
 TestBinomial.C:223
 TestBinomial.C:224
 TestBinomial.C:225
 TestBinomial.C:226
 TestBinomial.C:227
 TestBinomial.C:228
 TestBinomial.C:229
 TestBinomial.C:230
 TestBinomial.C:231
 TestBinomial.C:232
 TestBinomial.C:233
 TestBinomial.C:234
 TestBinomial.C:235
 TestBinomial.C:236
 TestBinomial.C:237
 TestBinomial.C:238
 TestBinomial.C:239
 TestBinomial.C:240
 TestBinomial.C:241
 TestBinomial.C:242
 TestBinomial.C:243
 TestBinomial.C:244
 TestBinomial.C:245
 TestBinomial.C:246
 TestBinomial.C:247
 TestBinomial.C:248
 TestBinomial.C:249
 TestBinomial.C:250
thumb