ROOT   Reference Guide
Searching...
No Matches
rs302_JeffreysPriorDemo.C File Reference

## Detailed Description

tutorial demonstrating and validates the RooJeffreysPrior class

Jeffreys's prior is an 'objective prior' based on formal rules. It is calculated from the Fisher information matrix.

The analytic form is not known for most PDFs, but it is for simple cases like the Poisson mean, Gaussian mean, Gaussian sigma.

This class uses numerical tricks to calculate the Fisher Information Matrix efficiently. In particular, it takes advantage of a property of the 'Asimov data' as described in Asymptotic formulae for likelihood-based tests of new physics Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells http://arxiv.org/abs/arXiv:1007.1727

This Demo has four parts:

1. TestJeffreysPriorDemo – validates Poisson mean case 1/sqrt(mu)
2. TestJeffreysGaussMean – validates Gaussian mean case
3. TestJeffreysGaussSigma – validates Gaussian sigma case 1/sigma
4. TestJeffreysGaussMeanAndSigma – demonstrates 2-d example
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Fitting -- RooAbsPdf::fitTo(p) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- using CPU computation library compiled with -mavx2
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_p_genData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
Minuit2Minimizer: Minimize with max-calls 500 convergence for edm < 1 strategy 1
Minuit2Minimizer : Valid minimum - status = 0
FVAL = -360.517018598809159
Edm = 4.13877261906962124e-17
Nfcn = 22
mu = 100 +/- 9.98317 (limited)
[#1] INFO:Fitting -- RooAbsPdf::fitTo(p) Calculating sum-of-weights-squared correction matrix for covariance matrix
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
RooDataHist::genData[x] = 100 bins (100 weights)
RooFitResult: minimized FCN value: -360.517, estimated distance to minimum: 3.66543e-18
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0 HESSE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu 1.0000e+02 +/- 1.00e+01
variance = 100.016
stdev = 10.0008
jeffreys = 0.0999921
[#1] INFO:NumericIntegration -- RooRealIntegral::init(jeffreys_Int[mu]) using numeric integrator RooAdaptiveGaussKronrodIntegrator1D to calculate Int(mu)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(Expected_Int[mu]) using numeric integrator RooIntegrator1D to calculate Int(mu)
#include "RooWorkspace.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMatrixDSym.h"
using namespace RooFit;
void rs302_JeffreysPriorDemo()
{
w.factory("Uniform::u(x[0,1])");
w.factory("mu[100,1,200]");
w.factory("ExtendPdf::p(u,mu)");
std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
std::unique_ptr<RooFitResult> res{w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};
asimov->Print();
res->Print();
TMatrixDSym cov = res->covarianceMatrix();
cout << "variance = " << (cov.Determinant()) << endl;
cout << "stdev = " << sqrt(cov.Determinant()) << endl;
cov.Invert();
cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
w.defineSet("poi", "mu");
w.defineSet("obs", "x");
RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
RooGenericPdf *test = new RooGenericPdf("Expected", "Expected = 1/#sqrt{#mu}", "1./sqrt(mu)",
*w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
plot->Draw();
auto legend = plot->BuildLegend();
legend->DrawClone();
}
//_________________________________________________
void TestJeffreysGaussMean()
{
w.factory("Gaussian::g(x[0,-20,20],mu[0,-5.,5],sigma[1,0,10])");
w.factory("n[10,.1,200]");
w.factory("ExtendPdf::p(g,n)");
w.var("sigma")->setConstant();
w.var("n")->setConstant();
std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
std::unique_ptr<RooFitResult> res{w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};
asimov->Print();
res->Print();
TMatrixDSym cov = res->covarianceMatrix();
cout << "variance = " << (cov.Determinant()) << endl;
cout << "stdev = " << sqrt(cov.Determinant()) << endl;
cov.Invert();
cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
w.defineSet("poi", "mu");
w.defineSet("obs", "x");
RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
// return;
RooGenericPdf *test = new RooGenericPdf("constant", "Expected = constant", "1", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
plot->Draw();
auto legend = plot->BuildLegend();
legend->DrawClone();
}
//_________________________________________________
void TestJeffreysGaussSigma()
{
// this one is VERY sensitive
// if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
// and you get really bizarre shapes
// if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
// and the PDF falls off too fast at high sigma
w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
w.factory("n[100,.1,2000]");
w.factory("ExtendPdf::p(g,n)");
// w.var("sigma")->setConstant();
w.var("mu")->setConstant();
w.var("n")->setConstant();
w.var("x")->setBins(301);
std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
std::unique_ptr<RooFitResult> res{w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};
asimov->Print();
res->Print();
TMatrixDSym cov = res->covarianceMatrix();
cout << "variance = " << (cov.Determinant()) << endl;
cout << "stdev = " << sqrt(cov.Determinant()) << endl;
cov.Invert();
cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
w.defineSet("poi", "sigma");
w.defineSet("obs", "x");
RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
RooGenericPdf *test = new RooGenericPdf("test", "Expected = #sqrt{2}/#sigma", "sqrt(2.)/sigma", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("sigma")->frame();
pi.plotOn(plot);
plot->Draw();
auto legend = plot->BuildLegend();
legend->DrawClone();
}
//_________________________________________________
void TestJeffreysGaussMeanAndSigma()
{
// this one is VERY sensitive
// if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
// and you get really bizarre shapes
// if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
// and the PDF falls off too fast at high sigma
w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1.,5.])");
w.factory("n[100,.1,2000]");
w.factory("ExtendPdf::p(g,n)");
w.var("n")->setConstant();
w.var("x")->setBins(301);
std::unique_ptr<RooDataHist> asimov{w.pdf("p")->generateBinned(*w.var("x"), ExpectedData())};
std::unique_ptr<RooFitResult> res{w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};
asimov->Print();
res->Print();
TMatrixDSym cov = res->covarianceMatrix();
cout << "variance = " << (cov.Determinant()) << endl;
cout << "stdev = " << sqrt(cov.Determinant()) << endl;
cov.Invert();
cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;
w.defineSet("poi", "mu,sigma");
w.defineSet("obs", "x");
RooJeffreysPrior pi("jeffreys", "jeffreys", *w.pdf("p"), *w.set("poi"), *w.set("obs"));
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
// return;
TCanvas *c1 = new TCanvas;
TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10, -5., 5), YVar(*w.var("sigma"), Binning(10, 1., 5.)));
Jeff2d->Draw("surf");
}
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
@ 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
Definition RooAbsPdf.h:124
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
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:45
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:3066
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

Definition in file rs302_JeffreysPriorDemo.C.