ROOT   Reference Guide
ProfileLikelihoodCalculator.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer 28/07/2008
3
4/*************************************************************************
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class RooStats::ProfileLikelihoodCalculator
13 \ingroup Roostats
14
15The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator
16(the interface class for tools which can produce both a RooStats HypoTestResult
17and ConfInterval). The tool uses the profile likelihood ratio as a test statistic,
18and assumes that Wilks' theorem is valid. Wilks' theorem states that \f$-2 \cdot \ln(\lambda) \f$
19(profile likelihood ratio) is asymptotically distributed as a \f$\chi^2 \f$ distribution
20with \f$N \f$ degrees of freedom. Thus, \f$p\f$-values can be
21constructed, and the profile likelihood ratio can be used to construct a
22LikelihoodInterval. (In the future, this class could be extended to use toy
23Monte Carlo to calibrate the distribution of the test statistic).
24
25Usage: It uses the interface of the CombinedCalculator, so that it can be
26configured by specifying:
27
28 - A model common model (*e.g.* a family of specific models, which includes both
29 the null and alternate)
30 - A data set
31 - A set of parameters of interest. The nuisance parameters will be all other
32 parameters of the model.
33 - A set of parameters which specify the null hypothesis (including values
34 and const/non-const status).
35
36The interface allows one to pass the model, data, and parameters either directly
37or via a ModelConfig class. The alternate hypothesis leaves the parameter free
38to take any value other than those specified by the null hypothesis. There is
39therefore no need to specify the alternate parameters.
40
41After configuring the calculator, one only needs to call GetHypoTest() (which
42will return a HypoTestResult pointer) or GetInterval() (which will return a
43ConfInterval pointer).
44
45This calculator can work with both one-dimensional intervals or multi-
46dimensional ones (contours).
47
48Note that for hypothesis tests, it is often better to use the
49AsymptoticCalculator, which can compute in addition the expected
50\f$p\f$-value using an Asimov data set.
51
52*/
53
55
57
60
61#include "RooFitResult.h"
62#include "RooRealVar.h"
63#include "RooProfileLL.h"
64#include "RooGlobalFunc.h"
65#include "RooMsgService.h"
66
68#include "RooMinimizer.h"
69//#include "RooProdPdf.h"
70
71using namespace std;
72
74
75using namespace RooFit;
76using namespace RooStats;
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// default constructor
81
82ProfileLikelihoodCalculator::ProfileLikelihoodCalculator() :
83 CombinedCalculator(), fFitResult(0), fGlobalFitDone(false)
84{
85}
86
88 double size, const RooArgSet* nullParams ) :
89 CombinedCalculator(data,pdf, paramsOfInterest, size, nullParams ),
90 fFitResult(0), fGlobalFitDone(false)
91{
92 // constructor from pdf and parameters
93 // the pdf must contain eventually the nuisance parameters
94}
95
98 fFitResult(0), fGlobalFitDone(false)
99{
100 // construct from a ModelConfig. Assume data model.GetPdf() will provide full description of model including
101 // constraint term on the nuisances parameters
102 assert(model.GetPdf() );
103}
104
105
106////////////////////////////////////////////////////////////////////////////////
107/// destructor
108/// cannot delete prod pdf because it will delete all the composing pdf's
109/// if (fOwnPdf) delete fPdf;
110/// fPdf = 0;
111
113 if (fFitResult) delete fFitResult;
114}
115
117 // reset and clear fit result
118 // to be called when a new model or data are set in the calculator
119 if (fFitResult) delete fFitResult;
120 fFitResult = 0;
121}
122
124 // perform a global fit of the likelihood letting with all parameter of interest and
125 // nuisance parameters
126 // keep the list of fitted parameters
127
128 DoReset();
129 RooAbsPdf * pdf = GetPdf();
131 if (!data || !pdf ) return 0;
132
133 // get all non-const parameters
134 RooArgSet* constrainedParams = pdf->getParameters(*data);
135 if (!constrainedParams) return 0;
136 RemoveConstantParameters(constrainedParams);
137
138 const auto& config = GetGlobalRooStatsConfig();
140 RooFit::Offset(config.useLikelihoodOffset) );
141
142 // check if global fit has been already done
143 if (fFitResult && fGlobalFitDone) {
144 delete constrainedParams;
145 return nll;
146 }
147
148 // calculate MLE
149 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGLobalFit - find MLE " << std::endl;
150
151 if (fFitResult) delete fFitResult;
153
154 // print fit result
155 if (fFitResult) {
157
158 if (fFitResult->status() != 0)
159 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGlobalFit - Global fit failed - status = " << fFitResult->status() << std::endl;
160 else
161 fGlobalFitDone = true;
162 }
163
164 delete constrainedParams;
165 return nll;
166}
167
169 // Minimizer the given NLL using the default options
170
171 const char * minimType = ""; // empty string to select RooMinimizer default
172 const char * minimAlgo = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str();
174 int level = ROOT::Math::MinimizerOptions::DefaultPrintLevel() -1;// RooFit level starts from -1
176 // do global fit and store fit result for further use
177
178 const auto& config = GetGlobalRooStatsConfig();
179
180 RooMinimizer minim(*nll);
181 minim.setStrategy(strategy);
182 minim.setEps(tolerance);
183 minim.setPrintLevel(level);
184 minim.optimizeConst(2); // to optimize likelihood calculations
185 minim.setEvalErrorWall(config.useEvalErrorWall);
186
187 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoMinimizeNLL - using " << minim.minimizerType()
188 << " / " << minimAlgo << " with strategy " << strategy << std::endl;
189
190 int status = -1;
191 for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
192 status = minim.minimize(minimType,minimAlgo);
193 if (status%1000 == 0) { // ignore erros from Improve
194 break;
195 } else if (tries < maxtries) {
196 cout << " ----> Doing a re-scan first" << endl;
197 minim.minimize(minimType,"Scan");
198 if (tries == 2) {
199 if (strategy == 0 ) {
200 cout << " ----> trying with strategy = 1" << endl;;
201 minim.setStrategy(1);
202 }
203 else
204 tries++; // skip this trial if strategy is already 1
205 }
206 if (tries == 3) {
207 cout << " ----> trying with improve" << endl;;
208 minimType = "Minuit";
210 }
211 }
212 }
213
214 RooFitResult * result = minim.save();
215
216
217 return result;
218}
219
220////////////////////////////////////////////////////////////////////////////////
221/// Main interface to get a RooStats::ConfInterval.
222/// It constructs a profile likelihood ratio, and uses that to construct a RooStats::LikelihoodInterval.
223
225// RooAbsPdf* pdf = fWS->pdf(fPdfName);
226// RooAbsData* data = fWS->data(fDataName);
227 RooAbsPdf * pdf = GetPdf();
229 if (!data || !pdf || fPOI.empty()) return 0;
230
231 RooArgSet* constrainedParams = pdf->getParameters(*data);
232 RemoveConstantParameters(constrainedParams);
233
234
235 /*
236 RooNLLVar* nll = new RooNLLVar("nll","",*pdf,*data, Extended(),Constrain(*constrainedParams));
237 RooProfileLL* profile = new RooProfileLL("pll","",*nll, *fPOI);
238 profile->addOwnedComponents(*nll) ; // to avoid memory leak
239 */
240
241 // do a global fit cloning the data
242 RooAbsReal * nll = DoGlobalFit();
243 if (!nll) return 0;
244
245 if (!fFitResult) {
246 delete nll;
247 return 0;
248 }
249
250 RooAbsReal* profile = nll->createProfile(fPOI);
251 profile->addOwnedComponents(*nll) ; // to avoid memory leak
252
253 // t.b.f. " RooProfileLL should keep and provide possibility to query on global minimum
254 // set POI to fit value (this will speed up profileLL calculation of global minimum)
255 const RooArgList & fitParams = fFitResult->floatParsFinal();
256 for (int i = 0; i < fitParams.getSize(); ++i) {
257 RooRealVar & fitPar = (RooRealVar &) fitParams[i];
258 RooRealVar * par = (RooRealVar*) fPOI.find( fitPar.GetName() );
259 if (par) {
260 par->setVal( fitPar.getVal() );
261 par->setError( fitPar.getError() );
262 }
263 }
264
265 // do this so profile will cache inside the absolute minimum and
266 // minimum values of nuisance parameters
267 // (no need to this here)
268 // profile->getVal();
269 //RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
270 // profile->Print();
271
272 TString name = TString("LikelihoodInterval_");// + TString(GetName() );
273
274 // make a list of fPOI with fit result values and pass to LikelihoodInterval class
275 // bestPOI is a cloned list of POI only with their best fit values
276 RooArgSet fitParSet(fitParams);
277 RooArgSet * bestPOI = new RooArgSet();
278 for (auto const *arg : fPOI){
279 RooAbsArg * p = fitParSet.find( arg->GetName() );
282 }
283 // fPOI contains the parameter of interest of the PL object
284 // and bestPOI contains a snapshot with the best fit values
285 LikelihoodInterval* interval = new LikelihoodInterval(name, profile, &fPOI, bestPOI);
286 interval->SetConfidenceLevel(1.-fSize);
287 delete constrainedParams;
288 return interval;
289}
290
291////////////////////////////////////////////////////////////////////////////////
292/// Main interface to get a HypoTestResult.
293/// It does two fits:
294/// 1. The first lets the null parameters float, so it's a maximum likelihood estimate.
295/// 2. The second is to the null model (fixing null parameters to their specified values): *e.g.* a conditional maximum likelihood.
296/// Since not all parameters are floating, this likelihood will be lower than the unconditional model.
297///
298/// The ratio of the likelihood obtained from the conditional MLE to the MLE is the profile likelihood ratio.
299/// Wilks' theorem is used to get \f$p\f$-values.
300
302// RooAbsPdf* pdf = fWS->pdf(fPdfName);
303// RooAbsData* data = fWS->data(fDataName);
304 RooAbsPdf * pdf = GetPdf();
306
307
308 if (!data || !pdf) return 0;
309
310 if (fNullParams.empty()) return 0;
311
312 // make a clone and ordered list since a vector will be associated to keep parameter values
313 // clone the list since first fit will changes the fNullParams values
314 RooArgList poiList;
315 poiList.addClone(fNullParams); // make a clone list
316
317
318 // do a global fit
319 RooAbsReal * nll = DoGlobalFit();
320 if (!nll) return 0;
321
322 if (!fFitResult) {
323 delete nll;
324 return 0;
325 }
326
327 RooArgSet* constrainedParams = pdf->getParameters(*data);
328 RemoveConstantParameters(constrainedParams);
329
330 double nLLatMLE = fFitResult->minNll();
331 // in case of using offset need to save offset value
332 double nlloffset = (RooStats::IsNLLOffset() ) ? nll->getVal() - nLLatMLE : 0;
333
334 // set POI to given values, set constant, calculate conditional MLE
335 std::vector<double> oldValues(poiList.getSize() );
336 for (unsigned int i = 0; i < oldValues.size(); ++i) {
337 RooRealVar * mytarget = (RooRealVar*) constrainedParams->find(poiList[i].GetName());
338 if (mytarget) {
339 oldValues[i] = mytarget->getVal();
340 mytarget->setVal( ( (RooRealVar&) poiList[i] ).getVal() );
341 mytarget->setConstant(true);
342 }
343 }
344
345
346
347 // perform the fit only if nuisance parameters are available
348 // get nuisance parameters
349 // nuisance parameters are the non const parameters from the likelihood parameters
350 RooArgSet nuisParams(*constrainedParams);
351
352 // need to remove the parameter of interest
353 RemoveConstantParameters(&nuisParams);
354
355 // check there are variable parameter in order to do a fit
356 bool existVarParams = false;
357 for (auto const *myarg : static_range_cast<RooRealVar *> (nuisParams)) {
358 if ( !myarg->isConstant() ) {
359 existVarParams = true;
360 break;
361 }
362 }
363
364 double nLLatCondMLE = nLLatMLE;
365 if (existVarParams) {
366 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit " << std::endl;
367
368 RooFitResult * fit2 = DoMinimizeNLL(nll);
369
370 // print fit result
371 if (fit2) {
372 nLLatCondMLE = fit2->minNll();
373 fit2->printStream( oocoutI(nullptr,Minimization), fit2->defaultPrintContents(0), fit2->defaultPrintStyle(0) );
374
375 if (fit2->status() != 0)
376 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypotest - Conditional fit failed - status = " << fit2->status() << std::endl;
377 }
378
379 }
380 else {
381 // get just the likelihood value (no need to do a fit since the likelihood is a constant function)
382 nLLatCondMLE = nll->getVal();
383 // this value contains the offset
384 if (RooStats::IsNLLOffset() ) nLLatCondMLE -= nlloffset;
385 }
386
387 // Use Wilks' theorem to translate -2 log lambda into a significance/p-value
388 double deltaNLL = std::max( nLLatCondMLE-nLLatMLE, 0.);
389
390 // get number of free parameter of interest
392 int ndf = poiList.getSize();
393
394 double pvalue = ROOT::Math::chisquared_cdf_c( 2* deltaNLL, ndf);
395
396 // in case of one dimension (1 poi) do the one-sided p-value (need to divide by 2)
397 if (ndf == 1) pvalue = 0.5 * pvalue;
398
399 TString name = TString("ProfileLRHypoTestResult_");// + TString(GetName() );
400 HypoTestResult* htr = new HypoTestResult(name, pvalue, 0 );
401
402 // restore previous value of poi
403 for (unsigned int i = 0; i < oldValues.size(); ++i) {
404 RooRealVar * mytarget = (RooRealVar*) constrainedParams->find(poiList[i].GetName());
405 if (mytarget) {
406 mytarget->setVal(oldValues[i] );
407 mytarget->setConstant(false);
408 }
409 }
410
411 delete constrainedParams;
412 delete nll;
413 return htr;
414
415}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
Definition: RooMsgService.h:51
#define oocoutI(o, a)
Definition: RooMsgService.h:49
#define oocoutP(o, a)
Definition: RooMsgService.h:50
#define ClassImp(name)
Definition: Rtypes.h:375
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
static const std::string & DefaultMinimizerAlgo()
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:71
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2185
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:541
bool empty() const
Int_t getSize() const
Return the number of elements in the collection.
const char * GetName() const override
Returns name of object.
virtual RooAbsArg * addClone(const RooAbsArg &var, bool silent=false)
Add a clone of the specified argument to list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:996
void setConstant(bool value=true)
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:480
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
Int_t defaultPrintContents(Option_t *opt) const override
Configure default contents to be printed.
StyleOption defaultPrintStyle(Option_t *opt) const override
Configure mapping of Print() arguments to RooPrintable print styles.
const RooArgList & floatParsFinal() const
Return list of floating parameters after fit.
Definition: RooFitResult.h:111
Int_t status() const
Return MINUIT status code.
Definition: RooFitResult.h:78
double minNll() const
Return minimized -log(L) value.
Definition: RooFitResult.h:99
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:43
void optimizeConst(int flag)
If flag is true, perform constant term optimization on function being minimized.
RooFitResult * save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
void setEvalErrorWall(bool flag)
Definition: RooMinimizer.h:80
void setEps(double eps)
Change MINUIT epsilon.
void setPrintLevel(int newLevel)
Change the MINUIT internal printing level.
std::string const & minimizerType() const
Definition: RooMinimizer.h:117
int minimize(const char *type, const char *alg=nullptr)
Minimise the function passed in the constructor.
void setStrategy(int strat)
Change MINUIT strategy to istrat.
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer,...
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
void setError(double value)
Definition: RooRealVar.h:64
double getError() const
Definition: RooRealVar.h:62
CombinedCalculator is an interface class for a tools which can produce both RooStats HypoTestResults ...
RooArgSet fNullParams
RooArgSet specifying null parameters for hypothesis test.
RooArgSet fPOI
RooArgSet specifying parameters of interest for interval.
RooArgSet fGlobalObs
RooArgSet specifying the global observables.
double fSize
size of the test (eg. specified rate of Type I error)
RooArgSet fConditionalObs
RooArgSet specifying the conditional observables.
HypoTestResult is a base class for results from hypothesis tests.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
void SetConfidenceLevel(double cl) override
set the confidence level for the interval (e.g 0.682 for a 1-sigma interval)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Definition: ModelConfig.h:231
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
static RooFitResult * DoMinimizeNLL(RooAbsReal *nll)
minimize likelihood
HypoTestResult * GetHypoTest() const override
Return the hypothesis test result obtained from the likelihood ratio of the maximum likelihood value ...
RooAbsReal * DoGlobalFit() const
perform a global fit
bool fGlobalFitDone
flag to control if a global fit has been done
~ProfileLikelihoodCalculator() override
destructor cannot delete prod pdf because it will delete all the composing pdf's if (fOwnPdf) delete ...
void DoReset() const
clear internal fit result
LikelihoodInterval * GetInterval() const override
Return a likelihood interval.
ProfileLikelihoodCalculator()
Default constructor (needed for I/O)
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:136
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg Offset(bool flag=true)
RooCmdArg CloneData(bool flag)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
double chisquared_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of the distribution with degrees of freedom (upp...
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ Minimization
Definition: RooGlobalFunc.h:61
Namespace for the RooStats classes.
Definition: Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:67
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
bool IsNLLOffset()
function returning if the flag to check if the flag to use NLLOffset is set