Logo ROOT  
Reference Guide
rs302_JeffreysPriorDemo.C File Reference

Detailed Description

tutorial demonstrating and validates the RooJeffreysPrior class

View in nbviewer Open in SWAN Jeffreys's prior is an 'objective prior' based on formal rules. It is calculated from the Fisher information matrix.

Read more: http://en.wikipedia.org/wiki/Jeffreys_prior

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
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_p_FOR_OBS_x with 0 entries
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (u)
**********
** 1 **SET PRINT 1
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 mu 1.00000e+02 1.99000e+01 1.00000e+00 2.00000e+02
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 1
**********
**********
** 5 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
** 6 **MIGRAD 500 1
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
FCN=-360.517 FROM MIGRAD STATUS=INITIATE 4 CALLS 5 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 mu 1.00000e+02 1.99000e+01 2.01361e-01 -5.66619e-05
ERR DEF= 0.5
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=-360.517 FROM MIGRAD STATUS=CONVERGED 24 CALLS 25 TOTAL
EDM=4.72209e-14 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 mu 1.00000e+02 9.98317e+00 1.31866e-03 -2.16215e-06
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 1 ERR DEF=0.5
1.000e+02
**********
** 7 **SET ERR 0.5
**********
**********
** 8 **SET PRINT 1
**********
**********
** 9 **HESSE 500
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=-360.517 FROM HESSE STATUS=OK 7 CALLS 32 TOTAL
EDM=9.50228e-17 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 mu 1.00000e+02 9.98317e+00 2.63731e-04 -5.02515e-03
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 1 ERR DEF=0.5
1.000e+02
[#1] INFO:Fitting -- RooAbsPdf::fitTo(p) Calculating sum-of-weights-squared correction matrix for covariance matrix
**********
** 10 **SET ERR 0.5
**********
**********
** 11 **SET PRINT 1
**********
**********
** 12 **HESSE 500
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=-360.517 FROM HESSE STATUS=OK 5 CALLS 37 TOTAL
EDM=8.98159e-17 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 mu 1.00000e+02 9.98311e+00 1.05492e-05 -5.02515e-03
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 1 ERR DEF=0.5
1.000e+02
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
RooDataHist::genData[x] = 100 bins (100 weights)
RooFitResult: minimized FCN value: -360.517, estimated distance to minimum: 8.98159e-17
covariance matrix quality: Unknown, matrix was externally provided
Status : MINIMIZE=0 HESSE=0 HESSE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu 1.0000e+02 +/- 1.00e+01
variance = 100.001
stdev = 10.0001
jeffreys = 0.0999994
[#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()
{
RooWorkspace w("w");
w.factory("Uniform::u(x[0,1])");
w.factory("mu[100,1,200]");
w.factory("ExtendPdf::p(u,mu)");
RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
asimov->Print();
res->Print();
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()
{
RooWorkspace w("w");
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();
RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
asimov->Print();
res->Print();
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
RooWorkspace w("w");
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);
RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
asimov->Print();
res->Print();
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
RooWorkspace w("w");
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);
RooDataHist *asimov = w.pdf("p")->generateBinned(*w.var("x"), ExpectedData());
RooFitResult *res = w.pdf("p")->fitTo(*asimov, Save(), SumW2Error(kTRUE));
asimov->Print();
res->Print();
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");
}
Author
Kyle Cranmer

Definition in file rs302_JeffreysPriorDemo.C.

RooAbsRealLValue::frame
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Definition: RooAbsRealLValue.cxx:199
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
TH1F.h
RooNumIntConfig.h
RooJeffreysPrior.h
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:207
TLegend.h
TObject::DrawClone
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad for instance with: gROOT->SetSelectedPad(gPad...
Definition: TObject.cxx:221
RooAbsData::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:182
TMatrixTSym
Definition: TMatrixDSymfwd.h:22
TCanvas.h
RooFit::YVar
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
Definition: RooGlobalFunc.cxx:240
RooFit::Binning
RooCmdArg Binning(const RooAbsBinning &binning)
Definition: RooGlobalFunc.cxx:82
RooGenericPdf.h
RooWorkspace::set
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
Definition: RooWorkspace.cxx:977
test
Definition: test.py:1
RooFitResult
Definition: RooFitResult.h:40
RooWorkspace::factory
RooFactoryWSTool & factory()
Return instance to factory tool.
Definition: RooWorkspace.cxx:2166
TMatrixDSym.h
kDashDotted
@ kDashDotted
Definition: TAttLine.h:48
RooDataHist
Definition: RooDataHist.h:39
RooFit
Definition: RooCFunction1Binding.h:29
RooAbsPdf.h
RooDataHist.h
RooPlot.h
RooPlot
Definition: RooPlot.h:44
RooWorkspace::pdf
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1277
sqrt
double sqrt(double)
RooRealVar.h
RooPlot::BuildLegend
std::unique_ptr< TLegend > BuildLegend() const
Build a legend that contains all objects that have been drawn on the plot.
Definition: RooPlot.cxx:1401
kRed
@ kRed
Definition: Rtypes.h:66
RooFitResult.h
RooWorkspace
Definition: RooWorkspace.h:43
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:56
TCanvas
Definition: TCanvas.h:23
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:73
TMatrixTSym::Invert
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
Definition: TMatrixTSym.cxx:961
TH1
Definition: TH1.h:57
RooGenericPdf
Definition: RooGenericPdf.h:25
RooFitResult::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
RooAbsRealLValue::setConstant
void setConstant(Bool_t value=kTRUE)
Definition: RooAbsRealLValue.h:115
RooWorkspace::var
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1295
RooAbsPdf::generateBinned
virtual RooDataHist * generateBinned(const RooArgSet &whatVars, Double_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none()) const
As RooAbsPdf::generateBinned(const RooArgSet&, const RooCmdArg&,const RooCmdArg&, const RooCmdArg&,...
Definition: RooAbsPdf.h:104
TMatrixTSym::Determinant
virtual Double_t Determinant() const
Definition: TMatrixTSym.cxx:935
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:57
RooRealVar::setBins
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:367
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooJeffreysPrior
Definition: RooJeffreysPrior.h:17
RooArgSet
Definition: RooArgSet.h:28
RooFitResult::covarianceMatrix
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
Definition: RooFitResult.cxx:1108
RooFit::ExpectedData
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:232
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2997
c1
return c1
Definition: legend1.C:41
RooAbsPdf::fitTo
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1261
RooWorkspace::defineSet
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
Definition: RooWorkspace.cxx:855