Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
50 using namespace RooFit;
51 
52 void 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 }