Logo ROOT   6.18/05
Reference Guide
JeffreysPriorDemo.C File Reference
Tutorials » RooStats Tutorials

Detailed Description

View in nbviewer Open in SWAN 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.

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
␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m
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)
[#0] ERROR:Integration -- RooAdaptiveGaussKronrodIntegrator1D::integral() ERROR: maximum number of subdivisions reached
[#1] INFO:NumericIntegration -- RooRealIntegral::init(test_Int[mu]) using numeric integrator RooIntegrator1D to calculate Int(mu)
#include "RooJeffreysPrior.h"
#include "RooWorkspace.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "TMatrixDSym.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooNumIntConfig.h"
#include "TH1F.h"
using namespace RooFit;
void 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();
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("test", "test", "1./sqrt(mu)", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
test->plotOn(plot, LineColor(kRed));
plot->Draw();
}
//_________________________________________________
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();
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("test", "test", "1", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
test->plotOn(plot, LineColor(kRed), LineStyle(kDotted));
plot->Draw();
}
//_________________________________________________
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();
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"));
pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps", 3);
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
RooGenericPdf *test = new RooGenericPdf("test", "test", "sqrt(2.)/sigma", *w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("sigma")->frame();
pi.plotOn(plot);
test->plotOn(plot, LineColor(kRed), LineStyle(kDotted));
plot->Draw();
}
//_________________________________________________
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();
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"));
pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps", 3);
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), YVar(*w.var("sigma"), Binning(10)));
Jeff2d->Draw("surf");
}
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
kRed
@ kRed
Definition: Rtypes.h:64
kDotted
@ kDotted
Definition: TAttLine.h:48
sqrt
double sqrt(double)
RooAbsData::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:162
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooDataHist
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooFitResult::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
RooFitResult::covarianceMatrix
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
Definition: RooFitResult.cxx:1089
RooGenericPdf
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Definition: RooGenericPdf.h:25
RooJeffreysPrior
RooJeffreysPrior.
Definition: RooJeffreysPrior.h:17
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:41
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:558
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
TCanvas
The Canvas class.
Definition: TCanvas.h:31
TH1
The TH1 histogram class.
Definition: TH1.h:56
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
TMatrixTSym::Determinant
virtual Double_t Determinant() const
Definition: TMatrixTSym.cxx:934
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:960
c1
return c1
Definition: legend1.C:41
RooFit
Template specialisation used in RooAbsArg:
Definition: RooCFunction1Binding.h:29
RooFit::Binning
RooCmdArg Binning(const RooAbsBinning &binning)
Definition: RooGlobalFunc.cxx:81
RooFit::YVar
RooCmdArg YVar(const RooAbsRealLValue &var, const RooCmdArg &arg=RooCmdArg::none())
Definition: RooGlobalFunc.cxx:236
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:204
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:184
RooFit::ExpectedData
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:228
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:55
RooFit::LineStyle
RooCmdArg LineStyle(Style_t style)
Definition: RooGlobalFunc.cxx:56
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:67
test
Definition: test.py:1
Author
Kyle Cranmer

Definition in file JeffreysPriorDemo.C.