Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs302_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 "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
50using namespace RooFit;
51
52void 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//_________________________________________________
90void 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//_________________________________________________
133void 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//_________________________________________________
182void 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}
const Bool_t kTRUE
Definition RtypesCore.h:91
@ kRed
Definition Rtypes.h:66
@ kDashDotted
Definition TAttLine.h:48
double sqrt(double)
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition RooAbsData.h:193
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
The RooDataHist is a container class to hold N-dimensional binned data.
Definition RooDataHist.h:37
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Implementation of Jeffrey's prior.
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:44
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
std::unique_ptr< TLegend > BuildLegend() const
Build a legend that contains all objects that have been drawn on the plot.
Definition RooPlot.cxx:1406
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:691
The RooWorkspace is a persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
virtual Double_t Determinant() const
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 YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
RooCmdArg SumW2Error(Bool_t flag)
RooCmdArg Save(Bool_t flag=kTRUE)
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
RooCmdArg Binning(const RooAbsBinning &binning)
RooCmdArg LineColor(Color_t color)
RooCmdArg LineStyle(Style_t style)
return c1
Definition legend1.C:41
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition test.py:1