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 std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
60
61 std::unique_ptr<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);
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 std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
100
101 std::unique_ptr<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);
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 std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
150
151 std::unique_ptr<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);
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 std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
198
199 std::unique_ptr<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}
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
@ kRed
Definition Rtypes.h:66
@ kDashDotted
Definition TAttLine.h:48
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1={}, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}, const RooCmdArg &arg9={}, const RooCmdArg &arg10={}) const override
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition RooAbsPdf.h:124
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Implementation of a probability density function that takes a RooArgList of servers and a C++ express...
Implementation of Jeffrey's prior.
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
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:59
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
Double_t Determinant() const override
TMatrixTSym< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:299
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg={})
RooCmdArg Save(bool flag=true)
RooCmdArg SumW2Error(bool flag)
RooCmdArg ExpectedData(bool flag=true)
RooCmdArg Binning(const RooAbsBinning &binning)
RooCmdArg LineColor(Color_t color)
RooCmdArg LineStyle(Style_t style)
return c1
Definition legend1.C:41
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26