Logo ROOT  

Reference Guide
rs302_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
#include "RooJeffreysPrior.h"
#include "RooWorkspace.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooNumIntConfig.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TMatrixDSym.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("Expected", "Expected = 1/#sqrt{#mu}", "1./sqrt(mu)",
*w.set("poi"));
TCanvas *c1 = new TCanvas;
RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
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();
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);
test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
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();
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);
test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
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();
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");
}
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
kRed
@ kRed
Definition: Rtypes.h:64
kDashDotted
@ kDashDotted
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:166
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:1108
RooGenericPdf
RooGenericPdf is a concrete implementation of a probability density function, which takes a RooArgLis...
Definition: RooGenericPdf.h:25
RooJeffreysPrior
Implementation of Jeffrey's prior.
Definition: RooJeffreysPrior.h:17
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
RooPlot::frame
static RooPlot * frame(const RooAbsRealLValue &var, Double_t xmin, Double_t xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition: RooPlot.cxx:277
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:1384
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
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:2998
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
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
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:239
RooFit::SumW2Error
RooCmdArg SumW2Error(Bool_t flag)
Definition: RooGlobalFunc.cxx:206
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:186
RooFit::ExpectedData
RooCmdArg ExpectedData(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:231
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 rs302_JeffreysPriorDemo.C.