Logo ROOT  
Reference Guide
 
Loading...
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/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
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
70
71using namespace RooFit;
72using namespace RooStats;
73
74
75////////////////////////////////////////////////////////////////////////////////
76/// default constructor
77
79
81 double size, const RooArgSet* nullParams ) :
83 fGlobalFitDone(false)
84{
85 // constructor from pdf and parameters
86 // the pdf must contain eventually the nuisance parameters
87}
88
91 fGlobalFitDone(false)
92{
93 // construct from a ModelConfig. Assume data model.GetPdf() will provide full description of model including
94 // constraint term on the nuisances parameters
95 assert(model.GetPdf() );
96}
97
98
99////////////////////////////////////////////////////////////////////////////////
100/// destructor
101/// cannot delete prod pdf because it will delete all the composing pdf's
102/// if (fOwnPdf) delete fPdf;
103/// fPdf = 0;
104
106
108 // reset and clear fit result
109 // to be called when a new model or data are set in the calculator
110 fFitResult.reset();
111}
112
114 // perform a global fit of the likelihood letting with all parameter of interest and
115 // nuisance parameters
116 // keep the list of fitted parameters
117
118 DoReset();
119 RooAbsPdf * pdf = GetPdf();
121 if (!data || !pdf ) return nullptr;
122
123 // get all non-const parameters
124 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
125 if (!constrainedParams) return nullptr;
127
128 const auto& config = GetGlobalRooStatsConfig();
130 RooFit::Offset(config.useLikelihoodOffset) )};
131
132 // check if global fit has been already done
133 if (fFitResult && fGlobalFitDone) {
134 return RooFit::makeOwningPtr(std::move(nll));
135 }
136
137 // calculate MLE
138 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGLobalFit - find MLE " << std::endl;
139
140 fFitResult = std::unique_ptr<RooFitResult>{DoMinimizeNLL(nll.get())};
141
142 // print fit result
143 if (fFitResult) {
144 fFitResult->printStream( oocoutI(nullptr,Minimization), fFitResult->defaultPrintContents(nullptr), fFitResult->defaultPrintStyle(nullptr) );
145
146 if (fFitResult->status() != 0) {
147 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoGlobalFit - Global fit failed - status = " << fFitResult->status() << std::endl;
148 } else {
149 fGlobalFitDone = true;
150 }
151 }
152
153 return RooFit::makeOwningPtr(std::move(nll));
154}
155
157 // Minimizer the given NLL using the default options
158
159 const char * minimType = ""; // empty string to select RooMinimizer default
162 int level = ROOT::Math::MinimizerOptions::DefaultPrintLevel() -1;// RooFit level starts from -1
164 // do global fit and store fit result for further use
165
166 const auto& config = GetGlobalRooStatsConfig();
167
168 RooMinimizer minim(*nll);
169 minim.setStrategy(strategy);
170 minim.setEps(tolerance);
171 minim.setPrintLevel(level);
172 minim.optimizeConst(2); // to optimize likelihood calculations
173 minim.setEvalErrorWall(config.useEvalErrorWall);
174
175 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::DoMinimizeNLL - using " << minim.minimizerType()
176 << " / " << minimAlgo << " with strategy " << strategy << std::endl;
177
178 int status = -1;
179 for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
180 status = minim.minimize(minimType,minimAlgo);
181 if (status%1000 == 0) { // ignore errors from Improve
182 break;
183 } else if (tries < maxtries) {
184 std::cout << " ----> Doing a re-scan first" << std::endl;
185 minim.minimize(minimType,"Scan");
186 if (tries == 2) {
187 if (strategy == 0 ) {
188 std::cout << " ----> trying with strategy = 1" << std::endl;
189 minim.setStrategy(1);
190 }
191 else
192 tries++; // skip this trial if strategy is already 1
193 }
194 if (tries == 3) {
195 std::cout << " ----> trying with improve" << std::endl;
196 minimType = "Minuit";
197 minimAlgo = "migradimproved";
198 }
199 }
200 }
201
202 return minim.save();
203}
204
205////////////////////////////////////////////////////////////////////////////////
206/// Main interface to get a RooStats::ConfInterval.
207/// It constructs a profile likelihood ratio, and uses that to construct a RooStats::LikelihoodInterval.
208
210// RooAbsPdf* pdf = fWS->pdf(fPdfName);
211// RooAbsData* data = fWS->data(fDataName);
212 RooAbsPdf * pdf = GetPdf();
214 if (!data || !pdf || fPOI.empty()) return nullptr;
215
216 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
218
219
220 /*
221 RooNLLVar* nll = new RooNLLVar("nll","",*pdf,*data, Extended(),Constrain(*constrainedParams));
222 RooProfileLL* profile = new RooProfileLL("pll","",*nll, *fPOI);
223 profile->addOwnedComponents(*nll) ; // to avoid memory leak
224 */
225
226 // do a global fit cloning the data
227 std::unique_ptr<RooAbsReal> nll{DoGlobalFit()};
228 if (!nll) return nullptr;
229
230 if (!fFitResult) {
231 return nullptr;
232 }
233
234 std::unique_ptr<RooAbsReal> profile{nll->createProfile(fPOI)};
235 profile->addOwnedComponents(std::move(nll)) ; // to avoid memory leak
236
237 // t.b.f. " RooProfileLL should keep and provide possibility to query on global minimum
238 // set POI to fit value (this will speed up profileLL calculation of global minimum)
239 const RooArgList & fitParams = fFitResult->floatParsFinal();
240 for (std::size_t i = 0; i < fitParams.size(); ++i) {
241 RooRealVar & fitPar = static_cast<RooRealVar &>( fitParams[i]);
242 RooRealVar * par = static_cast<RooRealVar*>(fPOI.find( fitPar.GetName() ));
243 if (par) {
244 par->setVal( fitPar.getVal() );
245 par->setError( fitPar.getError() );
246 }
247 }
248
249 // do this so profile will cache inside the absolute minimum and
250 // minimum values of nuisance parameters
251 // (no need to this here)
252 // profile->getVal();
253 //RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
254 // profile->Print();
255
256 TString name = TString("LikelihoodInterval_");// + TString(GetName() );
257
258 // make a list of fPOI with fit result values and pass to LikelihoodInterval class
259 // bestPOI is a cloned list of POI only with their best fit values
260 RooArgSet fitParSet(fitParams);
261 RooArgSet * bestPOI = new RooArgSet();
262 for (auto const *arg : fPOI){
263 RooAbsArg * p = fitParSet.find( arg->GetName() );
264 if (p) bestPOI->addClone(*p);
265 else bestPOI->addClone(*arg);
266 }
267 // fPOI contains the parameter of interest of the PL object
268 // and bestPOI contains a snapshot with the best fit values
269 LikelihoodInterval* interval = new LikelihoodInterval(name, profile.release(), &fPOI, bestPOI);
270 interval->SetConfidenceLevel(1.-fSize);
271 return interval;
272}
273
274////////////////////////////////////////////////////////////////////////////////
275/// Main interface to get a HypoTestResult.
276/// It does two fits:
277/// 1. The first lets the null parameters float, so it's a maximum likelihood estimate.
278/// 2. The second is to the null model (fixing null parameters to their specified values): *e.g.* a conditional maximum likelihood.
279/// Since not all parameters are floating, this likelihood will be lower than the unconditional model.
280///
281/// The ratio of the likelihood obtained from the conditional MLE to the MLE is the profile likelihood ratio.
282/// Wilks' theorem is used to get \f$p\f$-values.
283
285// RooAbsPdf* pdf = fWS->pdf(fPdfName);
286// RooAbsData* data = fWS->data(fDataName);
287 RooAbsPdf * pdf = GetPdf();
289
290
291 if (!data || !pdf) return nullptr;
292
293 if (fNullParams.empty()) return nullptr;
294
295 // make a clone and ordered list since a vector will be associated to keep parameter values
296 // clone the list since first fit will changes the fNullParams values
298 poiList.addClone(fNullParams); // make a clone list
299
300
301 // do a global fit
302 std::unique_ptr<RooAbsReal> nll{DoGlobalFit()};
303 if (!nll) return nullptr;
304
305 if (!fFitResult) {
306 return nullptr;
307 }
308
309 std::unique_ptr<RooArgSet> constrainedParams{pdf->getParameters(*data)};
311
312 double nLLatMLE = fFitResult->minNll();
313 // in case of using offset need to save offset value
314 double nlloffset = RooStats::NLLOffsetMode() == "initial" ? nll->getVal() - nLLatMLE : 0;
315
316 // set POI to given values, set constant, calculate conditional MLE
317 std::vector<double> oldValues(poiList.size() );
318 for (unsigned int i = 0; i < oldValues.size(); ++i) {
319 RooRealVar * mytarget = static_cast<RooRealVar*>(constrainedParams->find(poiList[i].GetName()));
320 if (mytarget) {
321 oldValues[i] = mytarget->getVal();
322 mytarget->setVal( ( static_cast<RooRealVar&>( poiList[i]) ).getVal() );
323 mytarget->setConstant(true);
324 }
325 }
326
327
328
329 // perform the fit only if nuisance parameters are available
330 // get nuisance parameters
331 // nuisance parameters are the non const parameters from the likelihood parameters
333
334 // need to remove the parameter of interest
336
337 // check there are variable parameter in order to do a fit
338 bool existVarParams = false;
339 for (auto const *myarg : static_range_cast<RooRealVar *> (nuisParams)) {
340 if ( !myarg->isConstant() ) {
341 existVarParams = true;
342 break;
343 }
344 }
345
346 double nLLatCondMLE = nLLatMLE;
347 if (existVarParams) {
348 oocoutP(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit " << std::endl;
349
350 std::unique_ptr<RooFitResult> fit2{DoMinimizeNLL(&*nll)};
351
352 // print fit result
353 if (fit2) {
354 nLLatCondMLE = fit2->minNll();
355 fit2->printStream( oocoutI(nullptr,Minimization), fit2->defaultPrintContents(nullptr), fit2->defaultPrintStyle(nullptr) );
356
357 if (fit2->status() != 0)
358 oocoutW(nullptr,Minimization) << "ProfileLikelihoodCalcultor::GetHypotest - Conditional fit failed - status = " << fit2->status() << std::endl;
359 }
360
361 }
362 else {
363 // get just the likelihood value (no need to do a fit since the likelihood is a constant function)
364 nLLatCondMLE = nll->getVal();
365 // this value contains the offset
366 if (RooStats::NLLOffsetMode() == "initial") nLLatCondMLE -= nlloffset;
367 }
368
369 // Use Wilks' theorem to translate -2 log lambda into a significance/p-value
370 double deltaNLL = std::max( nLLatCondMLE-nLLatMLE, 0.);
371
372 // get number of free parameter of interest
374 int ndf = poiList.size();
375
377
378 // in case of one dimension (1 poi) do the one-sided p-value (need to divide by 2)
379 if (ndf == 1) pvalue = 0.5 * pvalue;
380
381 TString name = TString("ProfileLRHypoTestResult_");// + TString(GetName() );
383
384 // restore previous value of poi
385 for (unsigned int i = 0; i < oldValues.size(); ++i) {
386 RooRealVar * mytarget = static_cast<RooRealVar*>(constrainedParams->find(poiList[i].GetName()));
387 if (mytarget) {
388 mytarget->setVal(oldValues[i] );
389 mytarget->setConstant(false);
390 }
391 }
392
393 return htr;
394
395}
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)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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
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
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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:24
Wrapper class around ROOT::Math::Minimizer that provides a seamless interface between the minimizer f...
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
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.
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)
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)
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 CodegenImpl.h:64
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
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
Namespace for the RooStats classes.
Definition CodegenImpl.h:58
void RemoveConstantParameters(RooArgSet *set)
std::string const & NLLOffsetMode()
Test what offsetting mode RooStats should use by default.
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...