ROOT   Reference Guide
rs302_JeffreysPriorDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook -js
4 /// \brief tutorial demonstrating and validates the RooJeffreysPrior class
5 ///
6 /// Jeffreys's prior is an 'objective prior' based on formal rules.
7 /// It is calculated from the Fisher information matrix.
8 ///
10 /// http://en.wikipedia.org/wiki/Jeffreys_prior
11 ///
12 /// The analytic form is not known for most PDFs, but it is for
13 /// simple cases like the Poisson mean, Gaussian mean, Gaussian sigma.
14 ///
15 /// This class uses numerical tricks to calculate the Fisher Information Matrix
16 /// efficiently. In particular, it takes advantage of a property of the
17 /// 'Asimov data' as described in
18 /// Asymptotic formulae for likelihood-based tests of new physics
19 /// Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells
20 /// http://arxiv.org/abs/arXiv:1007.1727
21 ///
22 /// This Demo has four parts:
23 /// 1. TestJeffreysPriorDemo -- validates Poisson mean case 1/sqrt(mu)
24 /// 2. TestJeffreysGaussMean -- validates Gaussian mean case
25 /// 3. TestJeffreysGaussSigma -- validates Gaussian sigma case 1/sigma
26 /// 4. TestJeffreysGaussMeanAndSigma -- demonstrates 2-d example
27 ///
28 /// \macro_image
29 /// \macro_output
30 /// \macro_code
31 ///
32 /// \author Kyle Cranmer
33
34 #include "RooJeffreysPrior.h"
35
36 #include "RooWorkspace.h"
37 #include "RooDataHist.h"
38 #include "RooGenericPdf.h"
39 #include "RooPlot.h"
40 #include "RooFitResult.h"
41 #include "RooRealVar.h"
42 #include "RooAbsPdf.h"
43 #include "RooNumIntConfig.h"
44
45 #include "TH1F.h"
46 #include "TCanvas.h"
47 #include "TLegend.h"
48 #include "TMatrixDSym.h"
49
50 using namespace RooFit;
51
52 void rs302_JeffreysPriorDemo()
53 {
54  RooWorkspace w("w");
55  w.factory("Uniform::u(x[0,1])");
56  w.factory("mu[100,1,200]");
57  w.factory("ExtendPdf::p(u,mu)");
58
59  RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
60
61  RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
62
63  asimov->Print();
64  res->Print();
65  TMatrixDSym cov = res->covarianceMatrix();
66  cout << "variance = " << (cov.Determinant()) << endl;
67  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
68  cov.Invert();
69  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
70
71  w.defineSet("poi", "mu");
72  w.defineSet("obs", "x");
73
74  RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
75
76  RooGenericPdf *test = new RooGenericPdf("Expected", "Expected = 1/#sqrt{#mu}", "1./sqrt(mu)",
77  *w.set("poi"));
78
79  TCanvas *c1 = new TCanvas;
80  RooPlot *plot = w.var("mu")->frame();
81  pi.plotOn(plot);
82  test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
83  plot->Draw();
84
85  auto legend = plot->BuildLegend();
86  legend->DrawClone();
87 }
88
89 //_________________________________________________
90 void TestJeffreysGaussMean()
91 {
92  RooWorkspace w("w");
93  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5.,5],sigma[1,0,10])");
94  w.factory("n[10,.1,200]");
95  w.factory("ExtendPdf::p(g,n)");
96  w.var("sigma")->setConstant();
97  w.var("n")->setConstant();
98
99  RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
100
101  RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
102
103  asimov->Print();
104  res->Print();
105  TMatrixDSym cov = res->covarianceMatrix();
106  cout << "variance = " << (cov.Determinant()) << endl;
107  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
108  cov.Invert();
109  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
110
111  w.defineSet("poi", "mu");
112  w.defineSet("obs", "x");
113
114  RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
115
116  const RooArgSet *temp = w.set("poi");
117  pi.getParameters(*temp)->Print();
118
119  // return;
120  RooGenericPdf *test = new RooGenericPdf("constant", "Expected = constant", "1", *w.set("poi"));
121
122  TCanvas *c1 = new TCanvas;
123  RooPlot *plot = w.var("mu")->frame();
124  pi.plotOn(plot);
125  test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
126  plot->Draw();
127
128  auto legend = plot->BuildLegend();
129  legend->DrawClone();
130 }
131
132 //_________________________________________________
133 void TestJeffreysGaussSigma()
134 {
135  // this one is VERY sensitive
136  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
137  // and you get really bizarre shapes
138  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
139  // and the PDF falls off too fast at high sigma
140  RooWorkspace w("w");
141  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
142  w.factory("n[100,.1,2000]");
143  w.factory("ExtendPdf::p(g,n)");
144  // w.var("sigma")->setConstant();
145  w.var("mu")->setConstant();
146  w.var("n")->setConstant();
147  w.var("x")->setBins(301);
148
149  RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
150
151  RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
152
153  asimov->Print();
154  res->Print();
155  TMatrixDSym cov = res->covarianceMatrix();
156  cout << "variance = " << (cov.Determinant()) << endl;
157  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
158  cov.Invert();
159  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
160
161  w.defineSet("poi", "sigma");
162  w.defineSet("obs", "x");
163
164  RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
165
166  const RooArgSet *temp = w.set("poi");
167  pi.getParameters(*temp)->Print();
168
169  RooGenericPdf *test = new RooGenericPdf("test", "Expected = #sqrt{2}/#sigma", "sqrt(2.)/sigma", *w.set("poi"));
170
171  TCanvas *c1 = new TCanvas;
172  RooPlot *plot = w.var("sigma")->frame();
173  pi.plotOn(plot);
174  test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
175  plot->Draw();
176
177  auto legend = plot->BuildLegend();
178  legend->DrawClone();
179 }
180
181 //_________________________________________________
182 void TestJeffreysGaussMeanAndSigma()
183 {
184  // this one is VERY sensitive
185  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
186  // and you get really bizarre shapes
187  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
188  // and the PDF falls off too fast at high sigma
189  RooWorkspace w("w");
190  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1.,5.])");
191  w.factory("n[100,.1,2000]");
192  w.factory("ExtendPdf::p(g,n)");
193
194  w.var("n")->setConstant();
195  w.var("x")->setBins(301);
196
197  RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
198
199  RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
200
201  asimov->Print();
202  res->Print();
203  TMatrixDSym cov = res->covarianceMatrix();
204  cout << "variance = " << (cov.Determinant()) << endl;
205  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
206  cov.Invert();
207  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
208
209  w.defineSet("poi", "mu,sigma");
210  w.defineSet("obs", "x");
211
212  RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
213
214  const RooArgSet *temp = w.set("poi");
215  pi.getParameters(*temp)->Print();
216  // return;
217
218  TCanvas *c1 = new TCanvas;
219  TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10, -5., 5), YVar(*w.var("sigma"), Binning(10, 1., 5.)));
220  Jeff2d->Draw("surf");
221 }
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TH1F.h
RooNumIntConfig.h
RooJeffreysPrior.h
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:207
TLegend.h
RooAbsData::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:182
TMatrixTSym< Double_t >
TCanvas.h
RooFit::YVar
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
Definition: RooGlobalFunc.cxx:240
RooFit::Binning
RooCmdArg Binning(const RooAbsBinning &binning)
Definition: RooGlobalFunc.cxx:82
RooGenericPdf.h
test
Definition: test.py:1
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
TMatrixDSym.h
RooPlot::frame
static RooPlot * frame(const RooAbsRealLValue &var, Double_t xmin, Double_t xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition: RooPlot.cxx:249
kDashDotted
@ kDashDotted
Definition: TAttLine.h:48
RooDataHist
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
TMatrixTSym::Determinant
virtual Double_t Determinant() const
Definition: TMatrixTSym.cxx:935
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
RooAbsPdf.h
RooDataHist.h
RooPlot.h
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
sqrt
double sqrt(double)
RooRealVar.h
RooPlot::BuildLegend
std::unique_ptr< TLegend > BuildLegend() const
Build a legend that contains all objects that have been drawn on the plot.
Definition: RooPlot.cxx:1401
kRed
@ kRed
Definition: Rtypes.h:66
RooFitResult.h
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:56
TCanvas
The Canvas class.
Definition: TCanvas.h:23
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:67
TH1
TH1 is the base class of all histogramm classes in ROOT.
Definition: TH1.h:58
RooGenericPdf
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Definition: RooGenericPdf.h:25
RooFitResult::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:57
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
TMatrixTSym::Invert
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
Definition: TMatrixTSym.cxx:961
RooJeffreysPrior
Implementation of Jeffrey's prior.
Definition: RooJeffreysPrior.h:17
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooFitResult::covarianceMatrix
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
Definition: RooFitResult.cxx:1108
RooFit::ExpectedData
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:232
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3050
c1
return c1
Definition: legend1.C:41