Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rs302_JeffreysPriorDemo.py
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 (C++ version), and P. P. (Python translation)
33
34import ROOT
35
36
38
39 w = ROOT.RooWorkspace("w")
40 w.factory("Uniform::u(x[0,1])")
41 w.factory("mu[100,1,200]")
42 w.factory("ExtendPdf::p(u,mu)")
43
44 asimov = w.pdf("p").generateBinned(w["x"], ExpectedData=True)
45
46 res = w.pdf("p").fitTo(asimov, Save=True, SumW2Error=True)
47
48 asimov.Print()
49 res.Print()
50 cov = res.covarianceMatrix()
51 print("variance = ", (cov.Determinant()))
52 print("stdev = ", cov.Determinant() ** 0.5)
53 cov.Invert()
54 print("jeffreys = ", cov.Determinant() ** 0.5)
55
56 w.defineSet("poi", "mu")
57 w.defineSet("obs", "x")
58
59 pi = ROOT.RooJeffreysPrior("jeffreys", "jeffreys", w.pdf("p"), w.set("poi"), w.set("obs"))
60
61 test = ROOT.RooGenericPdf("Expected", "Expected = 1/#sqrt#mu", "1./sqrt(mu)", w.set("poi"))
62
63 c1 = ROOT.TCanvas()
64 plot = w["mu"].frame()
65 pi.plotOn(plot)
66 test.plotOn(plot, LineColor="r", LineStyle=ROOT.kDashDotted)
67 plot.Draw()
68
69 legend = plot.BuildLegend()
70 legend.DrawClone()
71
72 c1.Update()
73 c1.Draw()
74 c1.SaveAs("rs302_JeffreysPriorDemo.1.png")
75
76
77# _________________________________________________
78def TestJeffreysGaussMean():
79
80 w = ROOT.RooWorkspace("w")
81 w.factory("Gaussian.g(x[0,-20,20],mu[0,-5.,5],sigma[1,0,10])")
82 w.factory("n[10,.1,200]")
83 w.factory("ExtendPdf.p(g,n)")
84 w.var("sigma").setConstant()
85 w.var("n").setConstant()
86
87 asimov = w.pdf("p").generateBinned(w.var("x"), ExpectedData())
88
89 resw.pdf("p").fitTo(asimov, Save(), SumW2Error(True))
90
91 asimov.Print()
92 res.Print()
93 cov = res.covarianceMatrix()
94 print("variance = ", (cov.Determinant()))
95 print("stdev = ", sqrt(cov.Determinant()))
96 cov.Invert()
97 print("jeffreys = ", sqrt(cov.Determinant()))
98
99 w.defineSet("poi", "mu")
100 w.defineSet("obs", "x")
101
102 pi = RooJeffreysPrior("jeffreys", "jeffreys", w.pdf("p"), w.set("poi"), w.set("obs"))
103
104 temp = w.set("poi")
105 pi.getParameters(temp).Print()
106
107 # return;
108 test = RooGenericPdf("constant", "Expected = constant", "1", w.set("poi"))
109
110 c2 = TCanvas()
111 plot = w.var("mu").frame()
112 pi.plotOn(plot)
113 test.plotOn(plot, LineColor(kRed), LineStyle(kDashDotted))
114 plot.Draw()
115
116 legend = plot.BuildLegend()
117 legend.DrawClone()
118 c2.Update()
119 c2.Draw()
120 c2.SaveAs("rs302_JeffreysPriorDemo.2.png")
121
122
123# _________________________________________________
124def TestJeffreysGaussSigma():
125
126 # this one is VERY sensitive
127 # if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
128 # and you get really bizarre shapes
129 # if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
130 # and the PDF falls off too fast at high sigma
131 w = RooWorkspace("w")
132 w.factory("Gaussian.g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])")
133 w.factory("n[100,.1,2000]")
134 w.factory("ExtendPdf.p(g,n)")
135 # w.var("sigma").setConstant()
136 w.var("mu").setConstant()
137 w.var("n").setConstant()
138 w.var("x").setBins(301)
139
140 asimov = w.pdf("p").generateBinned(w.var("x"), ExpectedData())
141
142 resw.pdf("p").fitTo(asimov, Save(), SumW2Error(True))
143
144 asimov.Print()
145 res.Print()
146 cov = res.covarianceMatrix()
147 print("variance = ", (cov.Determinant()))
148 print("stdev = ", sqrt(cov.Determinant()))
149 cov.Invert()
150 print("jeffreys = ", sqrt(cov.Determinant()))
151
152 w.defineSet("poi", "sigma")
153 w.defineSet("obs", "x")
154
155 pi = RooJeffreysPrior("jeffreys", "jeffreys", w.pdf("p"), w.set("poi"), w.set("obs"))
156
157 temp = w.set("poi")
158 pi.getParameters(temp).Print()
159
160 test = RooGenericPdf("test", "Expected = #sqrt2/#sigma", "sqrt(2.)/sigma", w.set("poi"))
161
162 c3 = TCanvas()
163 plot = w.var("sigma").frame()
164 pi.plotOn(plot)
165 test.plotOn(plot, LineColor(kRed), LineStyle(kDashDotted))
166 plot.Draw()
167
168 legend = plot.BuildLegend()
169 legend.DrawClone()
170 c3.Update()
171 c3.Draw()
172 c3.SaveAs("rs302_JeffreysPriorDemo.3.png")
173
174
175# _________________________________________________
176def TestJeffreysGaussMeanAndSigma():
177
178 # this one is VERY sensitive
179 # if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
180 # and you get really bizarre shapes
181 # if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
182 # and the PDF falls off too fast at high sigma
183 w = RooWorkspace("w")
184 w.factory("Gaussian.g(x[0,-20,20],mu[0,-5,5],sigma[1,1.,5.])")
185 w.factory("n[100,.1,2000]")
186 w.factory("ExtendPdf.p(g,n)")
187
188 w.var("n").setConstant()
189 w.var("x").setBins(301)
190
191 asimov = w.pdf("p").generateBinned(w.var("x"), ExpectedData())
192
193 resw.pdf("p").fitTo(asimov, Save(), SumW2Error(True))
194
195 asimov.Print()
196 res.Print()
197 cov = res.covarianceMatrix()
198 print("variance = ", (cov.Determinant()))
199 print("stdev = ", sqrt(cov.Determinant()))
200 cov.Invert()
201 print("jeffreys = ", sqrt(cov.Determinant()))
202
203 w.defineSet("poi", "mu,sigma")
204 w.defineSet("obs", "x")
205
206 pi = RooJeffreysPrior("jeffreys", "jeffreys", w.pdf("p"), w.set("poi"), w.set("obs"))
207
208 temp = w.set("poi")
209 pi.getParameters(temp).Print()
210 # return;
211
212 c4 = TCanvas()
213 Jeff2d = pi.createHistogram(
214 "2dJeffreys", w.var("mu"), Binning(10, -5.0, 5), YVar(w.var("sigma"), Binning(10, 1.0, 5.0))
215 )
216 Jeff2d.Draw("surf")
217 c4.Update()
218 c4.Draw()
219 c4.SaveAs("rs302_JeffreysPriorDemo.4.png")
220
221
void Print(GNN_Data &d, std::string txt="")
Implementation of a probability density function that takes a RooArgList of servers and a C++ express...
Implementation of Jeffrey's prior.
Persistable container for RooFit projects.
The Canvas class.
Definition TCanvas.h:23