From ce8b9c3d3c86d4c2f476b2621b2827ded684f7b2 Mon Sep 17 00:00:00 2001
From: Stephan Hageboeck <stephan.hageboeck@cern.ch>
Date: Mon, 16 Dec 2019 15:48:20 +0100
Subject: [PATCH] [RF] Stabilise RooJeffreysPrior.
Since it is running an internal fit, RooJeffreysPrior is unstable when a
parameter is close to its boundary. To fix this, the prior PDF is now
cloned, and the ranges of its parameters are extended to make fits converge
even at the boundary.
The JeffreysPrior tutorial has been renamed and cleaned up.
---
roofit/roofit/inc/RooJeffreysPrior.h | 33 ++--
roofit/roofit/src/RooJeffreysPrior.cxx | 162 ++++++------------
...sPriorDemo.C => rs302_JeffreysPriorDemo.C} | 38 ++--
3 files changed, 102 insertions(+), 131 deletions(-)
rename tutorials/roostats/{JeffreysPriorDemo.C => rs302_JeffreysPriorDemo.C} (86%)
diff --git a/roofit/roofit/inc/RooJeffreysPrior.h b/roofit/roofit/inc/RooJeffreysPrior.h
index 8326079a1b1..7c5177db01d 100644
--- a/roofit/roofit/inc/RooJeffreysPrior.h
+++ b/roofit/roofit/inc/RooJeffreysPrior.h
@@ -17,7 +17,7 @@ class RooArgList ;
class RooJeffreysPrior : public RooAbsPdf {
public:
- RooJeffreysPrior() ;
+ RooJeffreysPrior() { };
RooJeffreysPrior(const char *name, const char *title, RooAbsPdf& nominal, const RooArgList& paramSet, const RooArgList& obsSet) ;
virtual ~RooJeffreysPrior() ;
@@ -27,22 +27,31 @@ public:
const RooArgList& lowList() const { return _obsSet ; }
const RooArgList& paramList() const { return _paramSet ; }
- Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
- Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
-
protected:
- RooRealProxy _nominal; // The nominal value
- //RooAbsPdf* _nominal; // The nominal value
- RooArgList _ownedList ; // List of owned components
- RooListProxy _obsSet ; // Low-side variation
- RooListProxy _paramSet ; // interpolation parameters
- mutable TIterator* _paramIter ; //! Iterator over paramSet
- mutable TIterator* _obsIter ; //! Iterator over lowSet
+ RooPdfProxy _nominal; // Proxy to the PDF for this prior.
+ RooListProxy _obsSet ; // Observables of the PDF.
+ RooListProxy _paramSet ; // Parameters of the PDF.
Double_t evaluate() const;
- ClassDef(RooJeffreysPrior,1) // Sum of RooAbsReal objects
+private:
+ struct CacheElem : public RooAbsCacheElement {
+ public:
+ virtual ~CacheElem() = default;
+ // Payload
+ std::unique_ptr<RooAbsPdf> _pdf;
+ std::unique_ptr<RooArgSet> _pdfVariables;
+
+ virtual RooArgList containedArgs(Action) override {
+ RooArgList list(*_pdf);
+ list.add(*_pdfVariables, true);
+ return list;
+ }
+ };
+ mutable RooObjCacheManager _cacheMgr; //!
+
+ ClassDef(RooJeffreysPrior,2) // Sum of RooAbsReal objects
};
#endif
diff --git a/roofit/roofit/src/RooJeffreysPrior.cxx b/roofit/roofit/src/RooJeffreysPrior.cxx
index 7095346ec3d..b604ac11a8f 100644
--- a/roofit/roofit/src/RooJeffreysPrior.cxx
+++ b/roofit/roofit/src/RooJeffreysPrior.cxx
@@ -1,17 +1,17 @@
/** \class RooJeffreysPrior
\ingroup Roofit
-RooJeffreysPrior
-**/
-
+Implementation of Jeffrey's prior. This class estimates the fisher information matrix by generating
+a binned Asimov dataset from the supplied PDFs, fitting it, retrieving the covariance matrix and inverting
+it. It returns the square root of the determinant of this matrix.
+Numerical integration is used to normalise. Since each integration step requires fits to be run,
+evaluating complicated PDFs may take long.
-#include "RooFit.h"
-
-#include "Riostream.h"
-#include "Riostream.h"
-#include <math.h>
+Check the tutorial rs302_JeffreysPriorDemo.C for a demonstration with a simple PDF.
+**/
#include "RooJeffreysPrior.h"
+
#include "RooAbsReal.h"
#include "RooAbsPdf.h"
#include "RooErrorHandler.h"
@@ -23,8 +23,7 @@ RooJeffreysPrior
#include "RooFitResult.h"
#include "RooNumIntConfig.h"
#include "RooRealVar.h"
-
-#include "TError.h"
+#include "RooHelpers.h"
using namespace std;
@@ -33,54 +32,44 @@ ClassImp(RooJeffreysPrior);
using namespace RooFit;
////////////////////////////////////////////////////////////////////////////////
-///_obsSet("!obsSet","obs-side variation",this),
+/// Construct a new JeffreysPrior.
+/// \param[in] name Name of this object.
+/// \param[in] title Title (for plotting)
+/// \param[in] nominal The PDF to base Jeffrey's prior on.
+/// \param[in] paramSet Parameters of the PDF.
+/// \param[in] obsSet Observables of the PDF.
RooJeffreysPrior::RooJeffreysPrior(const char* name, const char* title,
RooAbsPdf& nominal,
const RooArgList& paramSet,
const RooArgList& obsSet) :
RooAbsPdf(name, title),
- _nominal("nominal","nominal",this,nominal,kFALSE,kFALSE),
- _obsSet("!obsSet","obs-side variation",this,kFALSE,kFALSE),
- _paramSet("!paramSet","high-side variation",this)
+ _nominal("nominal","nominal",this, nominal, false, false),
+ _obsSet("!obsSet","Observables",this, false, false),
+ _paramSet("!paramSet","Parameters",this),
+ _cacheMgr(this, 1, true, false)
{
- _obsIter = _obsSet.createIterator() ;
- _paramIter = _paramSet.createIterator() ;
-
-
- TIterator* inputIter1 = obsSet.createIterator() ;
- RooAbsArg* comp ;
- while((comp = (RooAbsArg*)inputIter1->Next())) {
+ for (const auto comp : obsSet) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
- << " in first list is not of type RooAbsReal" << endl ;
+ << " in observable list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_obsSet.add(*comp) ;
- // if (takeOwnership) {
- // _ownedList.addOwned(*comp) ;
- // }
}
- delete inputIter1 ;
- TIterator* inputIter3 = paramSet.createIterator() ;
- while((comp = (RooAbsArg*)inputIter3->Next())) {
+ for (const auto comp : paramSet) {
if (!dynamic_cast<RooAbsReal*>(comp)) {
coutE(InputArguments) << "RooJeffreysPrior::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
- << " in first list is not of type RooAbsReal" << endl ;
+ << " in parameter list is not of type RooAbsReal" << endl ;
RooErrorHandler::softAbort() ;
}
_paramSet.add(*comp) ;
- // if (takeOwnership) {
- // _ownedList.addOwned(*comp) ;
- // }
}
- delete inputIter3 ;
// use a different integrator by default.
if(paramSet.getSize()==1)
this->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D") ;
-
}
////////////////////////////////////////////////////////////////////////////////
@@ -90,21 +79,9 @@ RooJeffreysPrior::RooJeffreysPrior(const RooJeffreysPrior& other, const char* na
RooAbsPdf(other, name),
_nominal("!nominal",this,other._nominal),
_obsSet("!obsSet",this,other._obsSet),
- _paramSet("!paramSet",this,other._paramSet)
-{
- _obsIter = _obsSet.createIterator() ;
- _paramIter = _paramSet.createIterator() ;
-
- // Member _ownedList is intentionally not copy-constructed -- ownership is not transferred
-}
-
-////////////////////////////////////////////////////////////////////////////////
-/// Default constructor
-
-RooJeffreysPrior::RooJeffreysPrior()
+ _paramSet("!paramSet",this,other._paramSet),
+ _cacheMgr(this, 1, true, false)
{
- _obsIter = NULL;
- _paramIter = NULL;
}
@@ -113,8 +90,7 @@ RooJeffreysPrior::RooJeffreysPrior()
RooJeffreysPrior::~RooJeffreysPrior()
{
- if (_obsIter) delete _obsIter ;
- if (_paramIter) delete _paramIter ;
+
}
////////////////////////////////////////////////////////////////////////////////
@@ -122,65 +98,41 @@ RooJeffreysPrior::~RooJeffreysPrior()
Double_t RooJeffreysPrior::evaluate() const
{
- RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
- RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
- // create Asimov dataset
- // _paramSet.Print("v");
- // return sqrt(_paramSet.getRealValue("mu"));
- *(_nominal.arg().getVariables()) = _paramSet;
- /*
- cout << "_______________"<<endl;
- _paramSet.Print("v");
- _nominal->getVariables()->Print("v");
- cout << "_______________"<<endl;
- */
- RooDataHist* data = ((RooAbsPdf&)(_nominal.arg())).generateBinned(_obsSet,ExpectedData());
- // data->Print("v");
- //RooFitResult* res = _nominal->fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kTRUE));
- RooFitResult* res = ((RooAbsPdf&)(_nominal.arg())).fitTo(*data, Save(),PrintLevel(-1),Minos(kFALSE),SumW2Error(kFALSE));
- TMatrixDSym cov = res->covarianceMatrix();
- cov.Invert();
- double ret = sqrt(cov.Determinant());
-
- /*
- // for 1 parameter can avoid making TMatrix etc.
- // but number of params may be > 1 with others held constant
- if(_paramSet.getSize()==1){
- RooRealVar* var = (RooRealVar*) _paramSet.first();
- // also, the _paramSet proxy one does not pick up a different value
- cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
- // need to get the actual variable instance out of the pdf like below
- var = (RooRealVar*) _nominal->getVariables()->find(var->GetName());
- cout << "eval at "<< ret << " " << 1/(var->getError()) << endl;
- }
- */
+ RooHelpers::LocalChangeMsgLevel msgLvlRAII(RooFit::WARNING);
+
+
+ CacheElem* cacheElm = (CacheElem*) _cacheMgr.getObj(nullptr);
+ if (!cacheElm) {
+ //Internally, we have to enlarge the range of fit parameters to make
+ //fits converge even if we are close to the limit of a parameter. Therefore, we clone the pdf and its
+ //observables here. If something happens to the external PDF, the cache is wiped,
+ //and we start to clone again.
+ auto& pdf = _nominal.arg();
+ RooAbsPdf* clonePdf = static_cast<RooAbsPdf*>(pdf.cloneTree());
+ auto vars = clonePdf->getParameters(_obsSet);
+ for (auto varTmp : *vars) {
+ auto& var = static_cast<RooRealVar&>(*varTmp);
+ auto range = var.getRange();
+ double span = range.second - range.first;
+ var.setRange(range.first - 0.1*span, range.second + 0.1 * span);
+ }
- // res->Print();
- delete data;
- delete res;
- RooMsgService::instance().setGlobalKillBelow(msglevel);
+ cacheElm = new CacheElem;
+ cacheElm->_pdf.reset(clonePdf);
+ cacheElm->_pdfVariables.reset(vars);
- // cout << "eval at "<< ret << endl;
- // _paramSet.Print("v");
- return ret;
+ _cacheMgr.setObj(nullptr, cacheElm);
+ }
-}
+ auto& cachedPdf = *cacheElm->_pdf;
+ auto& pdfVars = *cacheElm->_pdfVariables;
+ pdfVars = _paramSet;
-////////////////////////////////////////////////////////////////////////////////
-/// if (matchArgs(allVars,analVars,x)) return 1 ;
-/// if (matchArgs(allVars,analVars,mean)) return 2 ;
-/// return 1;
+ std::unique_ptr<RooDataHist> data( cachedPdf.generateBinned(_obsSet,ExpectedData()) );
+ std::unique_ptr<RooFitResult> res( cachedPdf.fitTo(*data, Save(),PrintLevel(-1),Minos(false),SumW2Error(false)) );
+ TMatrixDSym cov = res->covarianceMatrix();
+ cov.Invert();
-Int_t RooJeffreysPrior::getAnalyticalIntegral(RooArgSet& /*allVars*/, RooArgSet& /*analVars*/, const char* /*rangeName*/) const
-{
- return 0 ;
+ return sqrt(cov.Determinant());
}
-////////////////////////////////////////////////////////////////////////////////
-
-Double_t RooJeffreysPrior::analyticalIntegral(Int_t code, const char* /*rangeName*/) const
-{
- R__ASSERT(code==1 );
- //cout << "evaluating analytic integral" << endl;
- return 1.;
-}
diff --git a/tutorials/roostats/JeffreysPriorDemo.C b/tutorials/roostats/rs302_JeffreysPriorDemo.C
similarity index 86%
rename from tutorials/roostats/JeffreysPriorDemo.C
rename to tutorials/roostats/rs302_JeffreysPriorDemo.C
index 0904ad8d5a1..c2f720a5f92 100644
--- a/tutorials/roostats/JeffreysPriorDemo.C
+++ b/tutorials/roostats/rs302_JeffreysPriorDemo.C
@@ -36,14 +36,16 @@
#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"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMatrixDSym.h"
using namespace RooFit;
@@ -71,22 +73,24 @@ void JeffreysPriorDemo()
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"));
+ RooGenericPdf *test = new RooGenericPdf("Expected", "Expected = 1/#sqrt{#mu}", "1./sqrt(mu)",
+ *w.set("poi"));
TCanvas *c1 = new TCanvas;
- // The method to compute the Jeffreys prior becomes unstable at the boundaries.
- // Therefore, we don't plot it all the way.
- RooPlot *plot = w.var("mu")->frame(Range(2, 199));
+ RooPlot *plot = w.var("mu")->frame();
pi.plotOn(plot);
- test->plotOn(plot, LineColor(kRed));
+ 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("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();
@@ -113,13 +117,16 @@ void TestJeffreysGaussMean()
pi.getParameters(*temp)->Print();
// return;
- RooGenericPdf *test = new RooGenericPdf("test", "test", "1", *w.set("poi"));
+ 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(kDotted));
+ test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
plot->Draw();
+
+ auto legend = plot->BuildLegend();
+ legend->DrawClone();
}
//_________________________________________________
@@ -159,13 +166,16 @@ void TestJeffreysGaussSigma()
const RooArgSet *temp = w.set("poi");
pi.getParameters(*temp)->Print();
- RooGenericPdf *test = new RooGenericPdf("test", "test", "sqrt(2.)/sigma", *w.set("poi"));
+ 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(kDotted));
+ test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));
plot->Draw();
+
+ auto legend = plot->BuildLegend();
+ legend->DrawClone();
}
//_________________________________________________
@@ -177,7 +187,7 @@ void TestJeffreysGaussMeanAndSigma()
// 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("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)");
@@ -206,6 +216,6 @@ void TestJeffreysGaussMeanAndSigma()
// return;
TCanvas *c1 = new TCanvas;
- TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10), YVar(*w.var("sigma"), Binning(10)));
+ TH1 *Jeff2d = pi.createHistogram("2dJeffreys", *w.var("mu"), Binning(10, -5., 5), YVar(*w.var("sigma"), Binning(10, 1., 5.)));
Jeff2d->Draw("surf");
}
--
GitLab