Logo ROOT  
Reference Guide
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
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"
49#include <cassert>
50#include <iostream>
51
52void TestBinomial(int nloop = 100, int nevts = 100, bool plot = false, bool debug = false, int seed = 111)
53{
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
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]);
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]);
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]);
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
316int main() {
317 TestBinomial();
318}
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
char * Form(const char *fmt,...)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
void GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double *x, double *ci, double cl=0.95, bool norm=false) const
get confidence intervals for an array of n points x.
Definition: FitResult.cxx:545
static void SetDefaultIntegrator(const char *name)
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
Binomial fitter for the division of two histograms.
The Canvas class.
Definition: TCanvas.h:27
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:701
1-Dim function class
Definition: TF1.h:210
virtual void SetParError(Int_t ipar, Double_t error)
Set error for parameter number ipar.
Definition: TF1.cxx:3475
virtual TH1 * GetHistogram() const
Return a pointer to the histogram used to visualise the function.
Definition: TF1.cxx:1567
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1913
virtual Double_t GetProb() const
Return the fit probability.
Definition: TF1.cxx:1938
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2505
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3500
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes.
Definition: TF1.cxx:1350
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:628
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:506
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:31
virtual void Print(Option_t *option="") const
Print result of the fit, by default chi2, parameter values and errors.
Definition: TFitResult.cxx:44
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
Double_t * GetEY() const
Definition: TGraphErrors.h:68
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9741
The TH1 histogram class.
Definition: TH1.h:56
virtual void FillRandom(const char *fname, Int_t ntimes=5000)
Fill histogram following distribution in function fname.
Definition: TH1.cxx:3445
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7086
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3808
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition: TH1.cxx:778
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
Double_t GetRMS(Int_t axis=1) const
Definition: TH1.h:311
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:3045
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Bool_t Divide(TF1 *f1, Double_t c1=1)
Performs the operation: this = this/(c1*f1) if errors are defined (see TH1::Sumw2),...
Definition: TH1.cxx:2753
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8476
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
An array of TObjects.
Definition: TObjArray.h:37
void Add(TObject *obj)
Definition: TObjArray.h:74
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad for instance with: gROOT->SetSelectedPad(gPad...
Definition: TObject.cxx:219
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1165
Random number generator class based on M.
Definition: TRandom3.h:27
Basic string class.
Definition: TString.h:131
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1590
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1542
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
int main(int argc, char **argv)
return c1
Definition: legend1.C:41
TGraphErrors * gr
Definition: legend1.C:25
TH1F * h1
Definition: legend1.C:5
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.