Logo ROOT   6.14/05
Reference Guide
JeffreysPriorDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook -js
4 /// 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 ///
9 /// Read more:
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 "TCanvas.h"
40 #include "RooPlot.h"
41 #include "RooFitResult.h"
42 #include "TMatrixDSym.h"
43 #include "RooRealVar.h"
44 #include "RooAbsPdf.h"
45 #include "RooNumIntConfig.h"
46 #include "TH1F.h"
47 
48 using namespace RooFit;
49 
50 void JeffreysPriorDemo(){
51  RooWorkspace w("w");
52  w.factory("Uniform::u(x[0,1])");
53  w.factory("mu[100,1,200]");
54  w.factory("ExtendPdf::p(u,mu)");
55 
56  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
57 
58  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));
59 
60  asimov->Print();
61  res->Print();
62  TMatrixDSym cov = res->covarianceMatrix();
63  cout << "variance = " << (cov.Determinant()) << endl;
64  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
65  cov.Invert();
66  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
67 
68  w.defineSet("poi","mu");
69  w.defineSet("obs","x");
70 
71  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
72 
73  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));
74 
75  TCanvas* c1 = new TCanvas;
76  RooPlot* plot = w.var("mu")->frame();
77  pi.plotOn(plot);
78  test->plotOn(plot,LineColor(kRed));
79  plot->Draw();
80 
81 }
82 
83 //_________________________________________________
84 void TestJeffreysGaussMean(){
85  RooWorkspace w("w");
86  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
87  w.factory("n[10,.1,200]");
88  w.factory("ExtendPdf::p(g,n)");
89  w.var("sigma")->setConstant();
90  w.var("n")->setConstant();
91 
92  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
93 
94  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));
95 
96  asimov->Print();
97  res->Print();
98  TMatrixDSym cov = res->covarianceMatrix();
99  cout << "variance = " << (cov.Determinant()) << endl;
100  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
101  cov.Invert();
102  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
103 
104  w.defineSet("poi","mu");
105  w.defineSet("obs","x");
106 
107  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
108 
109  const RooArgSet* temp = w.set("poi");
110  pi.getParameters(*temp)->Print();
111 
112  // return;
113  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));
114 
115  TCanvas* c1 = new TCanvas;
116  RooPlot* plot = w.var("mu")->frame();
117  pi.plotOn(plot);
118  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
119  plot->Draw();
120 }
121 
122 //_________________________________________________
123 void TestJeffreysGaussSigma(){
124  // this one is VERY sensitive
125  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
126  // and you get really bizarre shapes
127  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
128  // and the PDF falls off too fast at high sigma
129  RooWorkspace w("w");
130  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
131  w.factory("n[100,.1,2000]");
132  w.factory("ExtendPdf::p(g,n)");
133  // w.var("sigma")->setConstant();
134  w.var("mu")->setConstant();
135  w.var("n")->setConstant();
136  w.var("x")->setBins(301);
137 
138  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
139 
140  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));
141 
142  asimov->Print();
143  res->Print();
144  TMatrixDSym cov = res->covarianceMatrix();
145  cout << "variance = " << (cov.Determinant()) << endl;
146  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
147  cov.Invert();
148  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
149 
150  w.defineSet("poi","sigma");
151  w.defineSet("obs","x");
152 
153  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
154  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);
155 
156  const RooArgSet* temp = w.set("poi");
157  pi.getParameters(*temp)->Print();
158 
159  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));
160 
161  TCanvas* c1 = new TCanvas;
162  RooPlot* plot = w.var("sigma")->frame();
163  pi.plotOn(plot);
164  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
165  plot->Draw();
166 
167 
168 }
169 
170 //_________________________________________________
171 void TestJeffreysGaussMeanAndSigma(){
172  // this one is VERY sensitive
173  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
174  // and you get really bizarre shapes
175  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
176  // and the PDF falls off too fast at high sigma
177  RooWorkspace w("w");
178  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
179  w.factory("n[100,.1,2000]");
180  w.factory("ExtendPdf::p(g,n)");
181 
182  w.var("n")->setConstant();
183  w.var("x")->setBins(301);
184 
185  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
186 
187  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));
188 
189  asimov->Print();
190  res->Print();
191  TMatrixDSym cov = res->covarianceMatrix();
192  cout << "variance = " << (cov.Determinant()) << endl;
193  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
194  cov.Invert();
195  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
196 
197  w.defineSet("poi","mu,sigma");
198  w.defineSet("obs","x");
199 
200  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
201  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);
202 
203  const RooArgSet* temp = w.set("poi");
204  pi.getParameters(*temp)->Print();
205  // return;
206 
207  TCanvas* c1 = new TCanvas;
208  TH1* Jeff2d = pi.createHistogram("2dJeffreys",*w.var("mu"),Binning(10),YVar(*w.var("sigma"),Binning(10)));
209  Jeff2d->Draw("surf");
210 }
static constexpr double pi
RooCmdArg LineColor(Color_t color)
return c1
Definition: legend1.C:41
Definition: Rtypes.h:59
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
Definition: test.py:1
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Plot (project) PDF on specified frame.
Definition: RooAbsPdf.h:105
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
double sqrt(double)
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
RooCmdArg LineStyle(Style_t style)
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
virtual Double_t Determinant() const
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2974
RooCmdArg SumW2Error(Bool_t flag)
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:41
The Canvas class.
Definition: TCanvas.h:31
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:159
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
The TH1 histogram class.
Definition: TH1.h:56
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
RooCmdArg Save(Bool_t flag=kTRUE)
RooJeffreysPrior.
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Definition: RooGenericPdf.h:25
const Bool_t kTRUE
Definition: RtypesCore.h:87
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:558
RooCmdArg Binning(const RooAbsBinning &binning)