ROOT   Reference Guide
Searching...
No Matches
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 std::cout, std::endl;
72
74
75using namespace RooFit;
76using namespace RooStats;
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// default constructor
81
83
85 double size, const RooArgSet* nullParams ) :
86 CombinedCalculator(data,pdf, paramsOfInterest, size, nullParams ),
87 fGlobalFitDone(false)
88{
89 // constructor from pdf and parameters
90 // the pdf must contain eventually the nuisance parameters
91}
92
95 fGlobalFitDone(false)
96{
97 // construct from a ModelConfig. Assume data model.GetPdf() will provide full description of model including
98 // constraint term on the nuisances parameters
99 assert(model.GetPdf() );
100}
101
102
103////////////////////////////////////////////////////////////////////////////////
104/// destructor
105/// cannot delete prod pdf because it will delete all the composing pdf's
106/// if (fOwnPdf) delete fPdf;
107/// fPdf = 0;
108
110
112 // reset and clear fit result
113 // to be called when a new model or data are set in the calculator
114 fFitResult.reset();
115}
116
118 // perform a global fit of the likelihood letting with all parameter of interest and
119 // nuisance parameters
120 // keep the list of fitted parameters
121
122 DoReset();
123 RooAbsPdf * pdf = GetPdf();
125 if (!data || !pdf ) return nullptr;
126
127 // get all non-const parameters
128 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
129 if (!constrainedParams) return nullptr;
130 RemoveConstantParameters(&*constrainedParams);
131
132 const auto& config = GetGlobalRooStatsConfig();
133 std::unique_ptr<RooAbsReal> nll{pdf->createNLL(*data, CloneData(true), Constrain(*constrainedParams),ConditionalObservables(fConditionalObs), GlobalObservables(fGlobalObs),
134 RooFit::Offset(config.useLikelihoodOffset) )};
135
136 // check if global fit has been already done
137 if (fFitResult && fGlobalFitDone) {
138 return RooFit::makeOwningPtr(std::move(nll));
139 }
140
141 // calculate MLE
142 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGLobalFit - find MLE " << std::endl;
143
144 fFitResult = std::unique_ptr<RooFitResult>{DoMinimizeNLL(nll.get())};
145
146 // print fit result
147 if (fFitResult) {
148 fFitResult->printStream( oocoutI(nullptr,Minimization), fFitResult->defaultPrintContents(nullptr), fFitResult->defaultPrintStyle(nullptr) );
149
150 if (fFitResult->status() != 0) {
151 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGlobalFit - Global fit failed - status = " << fFitResult->status() << std::endl;
152 } else {
153 fGlobalFitDone = true;
154 }
155 }
156
157 return RooFit::makeOwningPtr(std::move(nll));
158}
159
161 // Minimizer the given NLL using the default options
162
163 const char * minimType = ""; // empty string to select RooMinimizer default
164 const char * minimAlgo = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str();
166 int level = ROOT::Math::MinimizerOptions::DefaultPrintLevel() -1;// RooFit level starts from -1
168 // do global fit and store fit result for further use
169
170 const auto& config = GetGlobalRooStatsConfig();
171
172 RooMinimizer minim(*nll);
173 minim.setStrategy(strategy);
174 minim.setEps(tolerance);
175 minim.setPrintLevel(level);
176 minim.optimizeConst(2); // to optimize likelihood calculations
177 minim.setEvalErrorWall(config.useEvalErrorWall);
178
179 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoMinimizeNLL - using " << minim.minimizerType()
180 << " / " << minimAlgo << " with strategy " << strategy << std::endl;
181
182 int status = -1;
183 for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
184 status = minim.minimize(minimType,minimAlgo);
185 if (status%1000 == 0) { // ignore errors from Improve
186 break;
187 } else if (tries < maxtries) {
188 cout << " ----> Doing a re-scan first" << endl;
189 minim.minimize(minimType,"Scan");
190 if (tries == 2) {
191 if (strategy == 0 ) {
192 cout << " ----> trying with strategy = 1" << endl;;
193 minim.setStrategy(1);
194 }
195 else
196 tries++; // skip this trial if strategy is already 1
197 }
198 if (tries == 3) {
199 cout << " ----> trying with improve" << endl;;
200 minimType = "Minuit";
202 }
203 }
204 }
205
206 return minim.save();
207}
208
209////////////////////////////////////////////////////////////////////////////////
210/// Main interface to get a RooStats::ConfInterval.
211/// It constructs a profile likelihood ratio, and uses that to construct a RooStats::LikelihoodInterval.
212
214// RooAbsPdf* pdf = fWS->pdf(fPdfName);
215// RooAbsData* data = fWS->data(fDataName);
216 RooAbsPdf * pdf = GetPdf();
218 if (!data || !pdf || fPOI.empty()) return nullptr;
219
220 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
221 RemoveConstantParameters(&*constrainedParams);
222
223
224 /*
225 RooNLLVar* nll = new RooNLLVar("nll","",*pdf,*data, Extended(),Constrain(*constrainedParams));
226 RooProfileLL* profile = new RooProfileLL("pll","",*nll, *fPOI);
227 profile->addOwnedComponents(*nll) ; // to avoid memory leak
228 */
229
230 // do a global fit cloning the data
231 std::unique_ptr<RooAbsReal> nll{DoGlobalFit()};
232 if (!nll) return nullptr;
233
234 if (!fFitResult) {
235 return nullptr;
236 }
237
238 std::unique_ptr<RooAbsReal> profile{nll->createProfile(fPOI)};
239 profile->addOwnedComponents(std::move(nll)) ; // to avoid memory leak
240
241 // t.b.f. " RooProfileLL should keep and provide possibility to query on global minimum
242 // set POI to fit value (this will speed up profileLL calculation of global minimum)
243 const RooArgList & fitParams = fFitResult->floatParsFinal();
244 for (std::size_t i = 0; i < fitParams.size(); ++i) {
245 RooRealVar & fitPar = static_cast<RooRealVar &>( fitParams[i]);
246 RooRealVar * par = static_cast<RooRealVar*>(fPOI.find( fitPar.GetName() ));
247 if (par) {
248 par->setVal( fitPar.getVal() );
249 par->setError( fitPar.getError() );
250 }
251 }
252
253 // do this so profile will cache inside the absolute minimum and
254 // minimum values of nuisance parameters
255 // (no need to this here)
256 // profile->getVal();
257 //RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
258 // profile->Print();
259
260 TString name = TString("LikelihoodInterval_");// + TString(GetName() );
261
262 // make a list of fPOI with fit result values and pass to LikelihoodInterval class
263 // bestPOI is a cloned list of POI only with their best fit values
264 RooArgSet fitParSet(fitParams);
265 RooArgSet * bestPOI = new RooArgSet();
266 for (auto const *arg : fPOI){
267 RooAbsArg * p = fitParSet.find( arg->GetName() );
270 }
271 // fPOI contains the parameter of interest of the PL object
272 // and bestPOI contains a snapshot with the best fit values
273 LikelihoodInterval* interval = new LikelihoodInterval(name, profile.release(), &fPOI, bestPOI);
274 interval->SetConfidenceLevel(1.-fSize);
275 return interval;
276}
277
278////////////////////////////////////////////////////////////////////////////////
279/// Main interface to get a HypoTestResult.
280/// It does two fits:
281/// 1. The first lets the null parameters float, so it's a maximum likelihood estimate.
282/// 2. The second is to the null model (fixing null parameters to their specified values): *e.g.* a conditional maximum likelihood.
283/// Since not all parameters are floating, this likelihood will be lower than the unconditional model.
284///
285/// The ratio of the likelihood obtained from the conditional MLE to the MLE is the profile likelihood ratio.
286/// Wilks' theorem is used to get \f$p\f$-values.
287
289// RooAbsPdf* pdf = fWS->pdf(fPdfName);
290// RooAbsData* data = fWS->data(fDataName);
291 RooAbsPdf * pdf = GetPdf();
293
294
295 if (!data || !pdf) return nullptr;
296
297 if (fNullParams.empty()) return nullptr;
298
299 // make a clone and ordered list since a vector will be associated to keep parameter values
300 // clone the list since first fit will changes the fNullParams values
301 RooArgList poiList;
302 poiList.addClone(fNullParams); // make a clone list
303
304
305 // do a global fit
306 std::unique_ptr<RooAbsReal> nll{DoGlobalFit()};
307 if (!nll) return nullptr;
308
309 if (!fFitResult) {
310 return nullptr;
311 }
312
313 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
314 RemoveConstantParameters(&*constrainedParams);
315
316 double nLLatMLE = fFitResult->minNll();
317 // in case of using offset need to save offset value
318 double nlloffset = (RooStats::IsNLLOffset() ) ? nll->getVal() - nLLatMLE : 0;
319
320 // set POI to given values, set constant, calculate conditional MLE
321 std::vector<double> oldValues(poiList.size() );
322 for (unsigned int i = 0; i < oldValues.size(); ++i) {
323 RooRealVar * mytarget = static_cast<RooRealVar*>(constrainedParams->find(poiList[i].GetName()));
324 if (mytarget) {
325 oldValues[i] = mytarget->getVal();
326 mytarget->setVal( ( static_cast<RooRealVar&>( poiList[i]) ).getVal() );
327 mytarget->setConstant(true);
328 }
329 }
330
331
332
333 // perform the fit only if nuisance parameters are available
334 // get nuisance parameters
335 // nuisance parameters are the non const parameters from the likelihood parameters
336 RooArgSet nuisParams(*constrainedParams);
337
338 // need to remove the parameter of interest
339 RemoveConstantParameters(&nuisParams);
340
341 // check there are variable parameter in order to do a fit
342 bool existVarParams = false;
343 for (auto const *myarg : static_range_cast<RooRealVar *> (nuisParams)) {
344 if ( !myarg->isConstant() ) {
345 existVarParams = true;
346 break;
347 }
348 }
349
350 double nLLatCondMLE = nLLatMLE;
351 if (existVarParams) {
352 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit " << std::endl;
353
354 std::unique_ptr<RooFitResult> fit2{DoMinimizeNLL(&*nll)};
355
356 // print fit result
357 if (fit2) {
358 nLLatCondMLE = fit2->minNll();
359 fit2->printStream( oocoutI(nullptr,Minimization), fit2->defaultPrintContents(nullptr), fit2->defaultPrintStyle(nullptr) );
360
361 if (fit2->status() != 0)
362 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypotest - Conditional fit failed - status = " << fit2->status() << std::endl;
363 }
364
365 }
366 else {
367 // get just the likelihood value (no need to do a fit since the likelihood is a constant function)
368 nLLatCondMLE = nll->getVal();
369 // this value contains the offset
370 if (RooStats::IsNLLOffset() ) nLLatCondMLE -= nlloffset;
371 }
372
373 // Use Wilks' theorem to translate -2 log lambda into a significance/p-value
374 double deltaNLL = std::max( nLLatCondMLE-nLLatMLE, 0.);
375
376 // get number of free parameter of interest
378 int ndf = poiList.size();
379
380 double pvalue = ROOT::Math::chisquared_cdf_c( 2* deltaNLL, ndf);
381
382 // in case of one dimension (1 poi) do the one-sided p-value (need to divide by 2)
383 if (ndf == 1) pvalue = 0.5 * pvalue;
384
385 TString name = TString("ProfileLRHypoTestResult_");// + TString(GetName() );
386 HypoTestResult* htr = new HypoTestResult(name, pvalue, 0 );
387
388 // restore previous value of poi
389 for (unsigned int i = 0; i < oldValues.size(); ++i) {
390 RooRealVar * mytarget = static_cast<RooRealVar*>(constrainedParams->find(poiList[i].GetName()));
391 if (mytarget) {
392 mytarget->setVal(oldValues[i] );
393 mytarget->setConstant(false);
394 }
395 }
396
397 return htr;
398
399}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define oocoutW(o, a)
#define oocoutI(o, a)
#define oocoutP(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition TGX11.cxx:110
static const std::string & DefaultMinimizerAlgo()
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
RooFit::OwningPtr< 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...
Storage_t::size_type size() const
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.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:163
void setConstant(bool value=true)
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:55
Wrapper class around ROOT::Fit:Fitter that provides a seamless interface between the minimizer functi...
void optimizeConst(int flag)
If flag is true, perform constant term optimization on function being minimized.
RooFit::OwningPtr< RooFitResult > save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
void setEvalErrorWall(bool flag)
void setEps(double eps)
Change MINUIT epsilon.
void setPrintLevel(int newLevel)
Change the MINUIT internal printing level.
std::string const & minimizerType() const
int minimize(const char *type, const char *alg=nullptr)
Minimise the function passed in the constructor.
void setStrategy(int istrat)
Change MINUIT strategy to istrat.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
void setError(double value)
Definition RooRealVar.h:60
double getError() const
Definition RooRealVar.h:58
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:35
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
static RooFit::OwningPtr< 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 ...
RooFit::OwningPtr< 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:139
RooCmdArg Offset(std::string const &mode)
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
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 JSONIO.h:26
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35
Namespace for the RooStats classes.
Definition Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
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