Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TestBinomial.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_fit
3 /// \notebook -js
4 /// Perform a fit to a set of data with binomial errors
5 /// like those derived from the division of two histograms.
6 /// Three different fits are performed and compared:
7 ///
8 /// - simple least square fit to the divided histogram obtained
9 /// from TH1::Divide with option b
10 /// - least square fit to the TGraphAsymmErrors obtained from
11 /// TGraphAsymmErrors::BayesDivide
12 /// - likelihood fit performed on the dividing histograms using
13 /// binomial statistics with the TBinomialEfficiency class
14 ///
15 /// The first two methods are biased while the last one is statistical correct.
16 /// Running the script passing an integer value n larger than 1, n fits are
17 /// performed and the bias are also shown.
18 /// To run the script :
19 ///
20 /// to show the bias performing 100 fits for 1000 events per "experiment"
21 ///
22 /// ~~~{.cpp}
23 /// root[0]: .x TestBinomial.C+
24 /// ~~~
25 ///
26 /// to show the bias performing 100 fits for 1000 events per "experiment"
27 ///
28 /// ~~~{.cpp}
29 /// .x TestBinomial.C+(100, 1000)
30 /// ~~~
31 ///
32 /// \macro_image
33 /// \macro_output
34 /// \macro_code
35 ///
36 /// \author Rene Brun
37 
38 #include "TBinomialEfficiencyFitter.h"
39 #include "TVirtualFitter.h"
40 #include "TH1.h"
41 #include "TRandom3.h"
42 #include "TF1.h"
43 #include "TFitResult.h"
44 #include "TStyle.h"
45 #include "TCanvas.h"
46 #include "TLegend.h"
47 #include "TPaveStats.h"
48 #include "Math/IntegratorOptions.h"
49 #include <cassert>
50 #include <iostream>
51 
52 void TestBinomial(int nloop = 100, int nevts = 100, bool plot = false, bool debug = false, int seed = 111)
53 {
54  gStyle->SetMarkerStyle(20);
55  gStyle->SetLineWidth(2.0);
56  gStyle->SetOptStat(11);
57 
58  TObjArray hbiasNorm;
59  hbiasNorm.Add(new TH1D("h0Norm", "Bias Histogram fit",100,-5,5));
60  hbiasNorm.Add(new TH1D("h1Norm","Bias Binomial fit",100,-5,5));
61  TObjArray hbiasThreshold;
62  hbiasThreshold.Add(new TH1D("h0Threshold", "Bias Histogram fit",100,-5,5));
63  hbiasThreshold.Add(new TH1D("h1Threshold","Bias Binomial fit",100,-5,5));
64  TObjArray hbiasWidth;
65  hbiasWidth.Add(new TH1D("h0Width", "Bias Histogram fit",100,-5,5));
66  hbiasWidth.Add(new TH1D("h1Width","Bias Binomial fit",100,-5,5));
67  TH1D* hChisquared = new TH1D("hChisquared",
68  "#chi^{2} probability (Baker-Cousins)", 200, 0.0, 1.0);
69 
70  TVirtualFitter::SetDefaultFitter("Minuit2");
71  ROOT::Math::IntegratorOneDimOptions::SetDefaultIntegrator("Gauss");
72 
73  // Note: in order to be able to use TH1::FillRandom() to generate
74  // pseudo-experiments, we use a trick: generate "selected"
75  // and "non-selected" samples independently. These are
76  // statistically independent and therefore can be safely
77  // added to yield the "before selection" sample.
78 
79 
80  // Define (arbitrarily?) a distribution of input events.
81  // Here: assume a x^(-2) distribution. Boundaries: [10, 100].
82 
83  Double_t xmin =10, xmax = 100;
84  TH1D* hM2D = new TH1D("hM2D", "x^(-2) denominator distribution",
85  45, xmin, xmax);
86  TH1D* hM2N = new TH1D("hM2N", "x^(-2) numerator distribution",
87  45, xmin, xmax);
88  TH1D* hM2E = new TH1D("hM2E", "x^(-2) efficiency",
89  45, xmin, xmax);
90 
91  TF1* fM2D = new TF1("fM2D", "(1-[0]/(1+exp(([1]-x)/[2])))/(x*x)",
92  xmin, xmax);
93  TF1* fM2N = new TF1("fM2N", "[0]/(1+exp(([1]-x)/[2]))/(x*x)",
94  xmin, xmax);
95  TF1* fM2Fit = new TF1("fM2Fit", "[0]/(1+exp(([1]-x)/[2]))",
96  xmin, xmax);
97  TF1* fM2Fit2 = 0;
98 
99  TRandom3 rb(seed);
100 
101  // First try: use a single set of parameters.
102  // For each try, we need to find the overall normalization
103 
104  Double_t normalization = 0.80;
105  Double_t threshold = 25.0;
106  Double_t width = 5.0;
107 
108  fM2D->SetParameter(0, normalization);
109  fM2D->SetParameter(1, threshold);
110  fM2D->SetParameter(2, width);
111  fM2N->SetParameter(0, normalization);
112  fM2N->SetParameter(1, threshold);
113  fM2N->SetParameter(2, width);
114  Double_t integralN = fM2N->Integral(xmin, xmax);
115  Double_t integralD = fM2D->Integral(xmin, xmax);
116  Double_t fracN = integralN/(integralN+integralD);
117  Int_t nevtsN = rb.Binomial(nevts, fracN);
118  Int_t nevtsD = nevts - nevtsN;
119 
120  std::cout << nevtsN << " " << nevtsD << std::endl;
121 
122  gStyle->SetOptFit(1111);
123 
124  // generate many times to see the bias
125  for (int iloop = 0; iloop < nloop; ++iloop) {
126 
127  // generate pseudo-experiments
128  hM2D->Reset();
129  hM2N->Reset();
130  hM2D->FillRandom(fM2D->GetName(), nevtsD);
131  hM2N->FillRandom(fM2N->GetName(), nevtsN);
132  hM2D->Add(hM2N);
133 
134  // construct the "efficiency" histogram
135  hM2N->Sumw2();
136  hM2E->Divide(hM2N, hM2D, 1, 1, "b");
137 
138  // Fit twice, using the same fit function.
139  // In the first (standard) fit, initialize to (arbitrary) values.
140  // In the second fit, use the results from the first fit (this
141  // makes it easier for the fit -- but the purpose here is not to
142  // show how easy or difficult it is to obtain results, but whether
143  // the CORRECT results are obtained or not!).
144 
145  fM2Fit->SetParameter(0, 0.5);
146  fM2Fit->SetParameter(1, 15.0);
147  fM2Fit->SetParameter(2, 2.0);
148  fM2Fit->SetParError(0, 0.1);
149  fM2Fit->SetParError(1, 1.0);
150  fM2Fit->SetParError(2, 0.2);
151  TH1 * hf = fM2Fit->GetHistogram();
152  // std::cout << "Function values " << std::endl;
153  // for (int i = 1; i <= hf->GetNbinsX(); ++i)
154  // std::cout << hf->GetBinContent(i) << " ";
155  // std::cout << std::endl;
156 
157  TCanvas* cEvt;
158  if (plot) {
159  cEvt = new TCanvas(Form("cEnv%d",iloop),
160  Form("plots for experiment %d", iloop),
161  1000, 600);
162  cEvt->Divide(1,2);
163  cEvt->cd(1);
164  hM2D->DrawCopy("HIST");
165  hM2N->SetLineColor(kRed);
166  hM2N->DrawCopy("HIST SAME");
167  cEvt->cd(2);
168  }
169  for (int fit = 0; fit < 2; ++fit) {
170  Int_t status = 0;
171  switch (fit) {
172  case 0:
173  {
174  // TVirtualPad * pad = gPad;
175  // new TCanvas();
176  // fM2Fit->Draw();
177  // gPad = pad;
178  TString optFit = "RN";
179  if (debug) optFit += TString("SV");
180  TFitResultPtr res = hM2E->Fit(fM2Fit, optFit);
181  if (plot) {
182  hM2E->DrawCopy("E");
183  fM2Fit->SetLineColor(kBlue);
184  fM2Fit->DrawCopy("SAME");
185  }
186  if (debug) res->Print();
187  status = res;
188  break;
189  }
190  case 1:
191  {
192  // if (fM2Fit2) delete fM2Fit2;
193  // fM2Fit2 = dynamic_cast<TF1*>(fM2Fit->Clone("fM2Fit2"));
194  fM2Fit2 = fM2Fit; // do not clone/copy the function
195  if (fM2Fit2->GetParameter(0) >= 1.0)
196  fM2Fit2->SetParameter(0, 0.95);
197  fM2Fit2->SetParLimits(0, 0.0, 1.0);
198 
199  // TVirtualPad * pad = gPad;
200  // new TCanvas();
201  // fM2Fit2->Draw();
202  // gPad = pad;
203 
204  TBinomialEfficiencyFitter bef(hM2N, hM2D);
205  TString optFit = "RI S";
206  if (debug) optFit += TString("V");
207  TFitResultPtr res = bef.Fit(fM2Fit2,optFit);
208  status = res;
209  if (status !=0) {
210  std::cerr << "Error performing binomial efficiency fit, result = "
211  << status << std::endl;
212  res->Print();
213  continue;
214  }
215  if (plot) {
216  fM2Fit2->SetLineColor(kRed);
217  fM2Fit2->DrawCopy("SAME");
218 
219  bool confint = (status == 0);
220  if (confint) {
221  // compute confidence interval on fitted function
222  auto htemp = fM2Fit2->GetHistogram();
223  ROOT::Fit::BinData xdata;
224  ROOT::Fit::FillData(xdata, fM2Fit2->GetHistogram() );
225  TGraphErrors gr(fM2Fit2->GetHistogram() );
226  res->GetConfidenceIntervals(xdata, gr.GetEY(), 0.68, false);
227  gr.SetFillColor(6);
228  gr.SetFillStyle(3005);
229  gr.DrawClone("4 same");
230  }
231  }
232  if (debug) {
233  res->Print();
234  }
235  }
236  }
237 
238  if (status != 0) break;
239 
240  Double_t fnorm = fM2Fit->GetParameter(0);
241  Double_t enorm = fM2Fit->GetParError(0);
242  Double_t fthreshold = fM2Fit->GetParameter(1);
243  Double_t ethreshold = fM2Fit->GetParError(1);
244  Double_t fwidth = fM2Fit->GetParameter(2);
245  Double_t ewidth = fM2Fit->GetParError(2);
246  if (fit == 1) {
247  fnorm = fM2Fit2->GetParameter(0);
248  enorm = fM2Fit2->GetParError(0);
249  fthreshold = fM2Fit2->GetParameter(1);
250  ethreshold = fM2Fit2->GetParError(1);
251  fwidth = fM2Fit2->GetParameter(2);
252  ewidth = fM2Fit2->GetParError(2);
253  hChisquared->Fill(fM2Fit2->GetProb());
254  }
255 
256  TH1D* h = dynamic_cast<TH1D*>(hbiasNorm[fit]);
257  h->Fill((fnorm-normalization)/enorm);
258  h = dynamic_cast<TH1D*>(hbiasThreshold[fit]);
259  h->Fill((fthreshold-threshold)/ethreshold);
260  h = dynamic_cast<TH1D*>(hbiasWidth[fit]);
261  h->Fill((fwidth-width)/ewidth);
262  }
263  }
264 
265 
266  TCanvas* c1 = new TCanvas("c1",
267  "Efficiency fit biases",10,10,1000,800);
268  c1->Divide(2,2);
269 
270  TH1D *h0, *h1;
271  c1->cd(1);
272  h0 = dynamic_cast<TH1D*>(hbiasNorm[0]);
273  h0->Draw("HIST");
274  h1 = dynamic_cast<TH1D*>(hbiasNorm[1]);
275  h1->SetLineColor(kRed);
276  h1->Draw("HIST SAMES");
277  TLegend* l1 = new TLegend(0.1, 0.75, 0.5, 0.9,
278  "plateau parameter", "ndc");
279  l1->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
280  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
281  l1->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
282  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
283  l1->Draw();
284 
285  c1->cd(2);
286  h0 = dynamic_cast<TH1D*>(hbiasThreshold[0]);
287  h0->Draw("HIST");
288  h1 = dynamic_cast<TH1D*>(hbiasThreshold[1]);
289  h1->SetLineColor(kRed);
290  h1->Draw("HIST SAMES");
291  TLegend* l2 = new TLegend(0.1, 0.75, 0.5, 0.9,
292  "threshold parameter", "ndc");
293  l2->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
294  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
295  l2->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
296  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
297  l2->Draw();
298 
299  c1->cd(3);
300  h0 = dynamic_cast<TH1D*>(hbiasWidth[0]);
301  h0->Draw("HIST");
302  h1 = dynamic_cast<TH1D*>(hbiasWidth[1]);
303  h1->SetLineColor(kRed);
304  h1->Draw("HIST SAMES");
305  TLegend* l3 = new TLegend(0.1, 0.75, 0.5, 0.9, "width parameter", "ndc");
306  l3->AddEntry(h0, Form("histogram: mean = %4.2f RMS = \
307  %4.2f", h0->GetMean(), h0->GetRMS()), "l");
308  l3->AddEntry(h1, Form("binomial : mean = %4.2f RMS = \
309  %4.2f", h1->GetMean(), h1->GetRMS()), "l");
310  l3->Draw();
311 
312  c1->cd(4);
313  hChisquared->Draw("HIST");
314 }
315 
316 int main() {
317  TestBinomial();
318 }