Logo ROOT  
Reference Guide
LikelihoodInterval.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id$
2// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11/*****************************************************************************
12 * Project: RooStats
13 * Package: RooFit/RooStats
14 * @(#)root/roofit/roostats:$Id$
15 * Authors:
16 * Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
17 *
18 *****************************************************************************/
19
20
21/** \class RooStats::LikelihoodInterval
22 \ingroup Roostats
23
24 LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
25 It implements a connected N-dimensional intervals based on the contour of a likelihood ratio.
26 The boundary of the interval is equivalent to a MINUIT/MINOS contour about the maximum likelihood estimator
27
28 The interval does not need to be an ellipse (eg. it is not the HESSE error matrix).
29 The level used to make the contour is the same as that used in MINOS, eg. it uses Wilks' theorem,
30 which states that under certain regularity conditions the function -2* log (profile likelihood ratio) is asymptotically distributed as a chi^2 with N-dof, where
31 N is the number of parameters of interest.
32
33
34 Note, a boundary on the parameter space (eg. s>= 0) or a degeneracy (eg. mass of signal if Nsig = 0) can lead to violations of the conditions necessary for Wilks'
35 theorem to be true.
36
37 Also note, one can use any RooAbsReal as the function that will be used in the contour; however, the level of the contour
38 is based on Wilks' theorem as stated above.
39
40
41#### References
42
43* 1. F. James., Minuit.Long writeup D506, CERN, 1998.
44
45*/
46
47
50
51#include "RooAbsReal.h"
52#include "RooMsgService.h"
53
55#include "Math/Minimizer.h"
56#include "Math/Factory.h"
58#include "RooFunctor.h"
59#include "RooProfileLL.h"
60
61#include "TMinuitMinimizer.h"
62
63#include <string>
64#include <algorithm>
65#include <functional>
66#include <ctype.h> // need to use c version of toupper defined here
67
68
70
71using namespace RooStats;
72using namespace std;
73
74
75////////////////////////////////////////////////////////////////////////////////
76/// Default constructor with name and title
77
78LikelihoodInterval::LikelihoodInterval(const char* name) :
79 ConfInterval(name), fBestFitParams(0), fLikelihoodRatio(0), fConfidenceLevel(0.95)
80{
81}
82
83////////////////////////////////////////////////////////////////////////////////
84/// Alternate constructor taking a pointer to the profile likelihood ratio, parameter of interest and
85/// optionally a snapshot of best parameter of interest for interval
86
87LikelihoodInterval::LikelihoodInterval(const char* name, RooAbsReal* lr, const RooArgSet* params, RooArgSet * bestParams) :
89 fParameters(*params),
90 fBestFitParams(bestParams),
91 fLikelihoodRatio(lr),
92 fConfidenceLevel(0.95)
93{
94}
95
96
97////////////////////////////////////////////////////////////////////////////////
98/// Destructor
99
101{
102 if (fBestFitParams) delete fBestFitParams;
104}
105
106
107////////////////////////////////////////////////////////////////////////////////
108/// This is the main method to satisfy the RooStats::ConfInterval interface.
109/// It returns true if the parameter point is in the interval.
110
111bool LikelihoodInterval::IsInInterval(const RooArgSet &parameterPoint) const
112{
115 // Method to determine if a parameter point is in the interval
116 if( !this->CheckParameters(parameterPoint) ) {
117 std::cout << "parameters don't match" << std::endl;
119 return false;
120 }
121
122 // make sure likelihood ratio is set
123 if(!fLikelihoodRatio) {
124 std::cout << "likelihood ratio not set" << std::endl;
126 return false;
127 }
128
129
130
131 // set parameters
132 SetParameters(&parameterPoint, fLikelihoodRatio->getVariables() );
133
134
135 // evaluate likelihood ratio, see if it's bigger than threshold
136 if (fLikelihoodRatio->getVal()<0){
137 std::cout << "The likelihood ratio is < 0, indicates a bad minimum or numerical precision problems. Will return true" << std::endl;
139 return true;
140 }
141
142
143 // here we use Wilks' theorem.
144 if ( TMath::Prob( 2* fLikelihoodRatio->getVal(), parameterPoint.getSize()) < (1.-fConfidenceLevel) ){
146 return false;
147 }
148
149
151
152 return true;
153
154}
155
156////////////////////////////////////////////////////////////////////////////////
157/// returns list of parameters
158
160{
161 return new RooArgSet(fParameters);
162}
163
164////////////////////////////////////////////////////////////////////////////////
165/// check that the parameters are correct
166
167bool LikelihoodInterval::CheckParameters(const RooArgSet &parameterPoint) const
168{
169 if (parameterPoint.getSize() != fParameters.getSize() ) {
170 std::cout << "size is wrong, parameters don't match" << std::endl;
171 return false;
172 }
173 if ( ! parameterPoint.equals( fParameters ) ) {
174 std::cout << "size is ok, but parameters don't match" << std::endl;
175 return false;
176 }
177 return true;
178}
179
180
181
182////////////////////////////////////////////////////////////////////////////////
183/// Compute lower limit, check first if limit has been computed
184/// status is a boolean flag which will b set to false in case of error
185/// and is true if calculation is successful
186/// in case of error return also a lower limit value of zero
187
188double LikelihoodInterval::LowerLimit(const RooRealVar& param, bool & status)
189{
190 double lower = 0;
191 double upper = 0;
192 status = FindLimits(param, lower, upper);
193 return lower;
194}
195
196////////////////////////////////////////////////////////////////////////////////
197/// Compute upper limit, check first if limit has been computed
198/// status is a boolean flag which will b set to false in case of error
199/// and is true if calculation is successful
200/// in case of error return also a lower limit value of zero
201
202double LikelihoodInterval::UpperLimit(const RooRealVar& param, bool & status)
203{
204 double lower = 0;
205 double upper = 0;
206 status = FindLimits(param, lower, upper);
207 return upper;
208}
209
210
212 // reset map with cached limits - called every time the test size or CL has been changed
213 fLowerLimits.clear();
214 fUpperLimits.clear();
215}
216
217
219 // internal function to create minimizer object needed to find contours or interval limits
220 // (running MINOS).
221 // Minimizer must be Minuit or Minuit2
222
223 RooProfileLL * profilell = dynamic_cast<RooProfileLL*>(fLikelihoodRatio);
224 if (!profilell) return false;
225
226 RooAbsReal & nll = profilell->nll();
227 // bind the nll function in the right interface for the Minimizer class
228 // as a function of only the parameters (poi + nuisance parameters)
229
230 RooArgSet * partmp = profilell->getVariables();
231 // need to remove constant parameters
233
234 RooArgList params(*partmp);
235 delete partmp;
236
237 // need to restore values and errors for POI
238 if (fBestFitParams) {
239 for (int i = 0; i < params.getSize(); ++i) {
240 RooRealVar & par = (RooRealVar &) params[i];
241 RooRealVar * fitPar = (RooRealVar *) (fBestFitParams->find(par.GetName() ) );
242 if (fitPar) {
243 par.setVal( fitPar->getVal() );
244 par.setError( fitPar->getError() );
245 }
246 }
247 }
248
249 const auto& config = GetGlobalRooStatsConfig();
250
251 // now do binding of NLL with a functor for Minimizer
252 if (config.useLikelihoodOffset) {
253 ccoutI(InputArguments) << "LikelihoodInterval: using nll offset - set all RooAbsReal to hide the offset " << std::endl;
254 RooAbsReal::setHideOffset(false); // need to keep this false
255 }
256 fFunctor = std::make_shared<RooFunctor>(nll, RooArgSet(), params);
257
259 std::transform(minimType.begin(), minimType.end(), minimType.begin(), (int(*)(int)) tolower );
260 *minimType.begin() = toupper( *minimType.begin());
261
262 if (minimType != "Minuit" && minimType != "Minuit2") {
263 ccoutE(InputArguments) << minimType << " is wrong type of minimizer for getting interval limits or contours - must use Minuit or Minuit2" << std::endl;
264 return false;
265 }
266 // do not use static instance of TMInuit which could interfere with RooFit
267 if (minimType == "Minuit") TMinuitMinimizer::UseStaticMinuit(false);
268 // create minimizer class
269 fMinimizer = std::shared_ptr<ROOT::Math::Minimizer>(ROOT::Math::Factory::CreateMinimizer(minimType, "Migrad"));
270
271 if (!fMinimizer.get()) return false;
272
273 fMinFunc = std::static_pointer_cast<ROOT::Math::IMultiGenFunction>(
275 fMinimizer->SetFunction(*fMinFunc);
276
277 // set minimizer parameters
278 assert( params.getSize() == int(fMinFunc->NDim()) );
279
280 for (unsigned int i = 0; i < fMinFunc->NDim(); ++i) {
281 RooRealVar & v = (RooRealVar &) params[i];
282 fMinimizer->SetLimitedVariable( i, v.GetName(), v.getVal(), v.getError(), v.getMin(), v.getMax() );
283 }
284 // for finding the contour need to find first global minimum
285 bool iret = fMinimizer->Minimize();
286 if (!iret || fMinimizer->X() == 0) {
287 ccoutE(Minimization) << "Error: Minimization failed " << std::endl;
288 return false;
289 }
290
291 //std::cout << "print minimizer result..........." << std::endl;
292 //fMinimizer->PrintResults();
293
294 return true;
295}
296
297bool LikelihoodInterval::FindLimits(const RooRealVar & param, double &lower, double & upper)
298{
299 // Method to find both lower and upper limits using MINOS
300 // If cached values exist (limits have been already found) return them in that case
301 // check first if limit has been computed
302 // otherwise compute limit using MINOS
303 // in case of failure lower and upper will maintain previous value (will not be modified)
304
305 std::map<std::string, double>::const_iterator itrl = fLowerLimits.find(param.GetName());
306 std::map<std::string, double>::const_iterator itru = fUpperLimits.find(param.GetName());
307 if ( itrl != fLowerLimits.end() && itru != fUpperLimits.end() ) {
308 lower = itrl->second;
309 upper = itru->second;
310 return true;
311 }
312
313
316 RooArgList params(*partmp);
317 delete partmp;
318 int ix = params.index(&param);
319 if (ix < 0 ) {
320 ccoutE(InputArguments) << "Error - invalid parameter " << param.GetName() << " specified for finding the interval limits " << std::endl;
321 return false;
322 }
323
324 bool ret = true;
325 if (!fMinimizer.get()) ret = CreateMinimizer();
326 if (!ret) {
327 ccoutE(Eval) << "Error returned from minimization of likelihood function - cannot find interval limits " << std::endl;
328 return false;
329 }
330
331 assert(fMinimizer.get());
332
333 // getting a 1D interval so ndf = 1
334 double err_level = TMath::ChisquareQuantile(ConfidenceLevel(),1); // level for -2log LR
335 err_level = err_level/2; // since we are using -log LR
336 fMinimizer->SetErrorDef(err_level);
337
338 unsigned int ivarX = ix;
339
340 double elow = 0;
341 double eup = 0;
342 ret = fMinimizer->GetMinosError(ivarX, elow, eup );
343 if (!ret) {
344 ccoutE(Minimization) << "Error running Minos for parameter " << param.GetName() << std::endl;
345 return false;
346 }
347
348 // WHEN error is zero normally is at limit
349 if (elow == 0) {
350 lower = param.getMin();
351 ccoutW(Minimization) << "Warning: lower value for " << param.GetName() << " is at limit " << lower << std::endl;
352 }
353 else
354 lower = fMinimizer->X()[ivarX] + elow; // elow is negative
355
356 if (eup == 0) {
357 ccoutW(Minimization) << "Warning: upper value for " << param.GetName() << " is at limit " << upper << std::endl;
358 upper = param.getMax();
359 }
360 else
361 upper = fMinimizer->X()[ivarX] + eup;
362
363 // store limits in the map
364 // minos return error limit = minValue +/- error
365 fLowerLimits[param.GetName()] = lower;
366 fUpperLimits[param.GetName()] = upper;
367
368 return true;
369}
370
371
372Int_t LikelihoodInterval::GetContourPoints(const RooRealVar & paramX, const RooRealVar & paramY, double * x, double *y, Int_t npoints ) {
373 // use Minuit to find the contour of the likelihood function at the desired CL
374
375 // check the parameters
376 // variable index in minimizer
377 // is index in the RooArgList obtained from the profileLL variables
380 RooArgList params(*partmp);
381 delete partmp;
382 int ix = params.index(&paramX);
383 int iy = params.index(&paramY);
384 if (ix < 0 || iy < 0) {
385 coutE(InputArguments) << "LikelihoodInterval - Error - invalid parameters specified for finding the contours; parX = " << paramX.GetName()
386 << " parY = " << paramY.GetName() << std::endl;
387 return 0;
388 }
389
390 bool ret = true;
391 if (!fMinimizer.get()) ret = CreateMinimizer();
392 if (!ret) {
393 coutE(Eval) << "LikelihoodInterval - Error returned creating minimizer for likelihood function - cannot find contour points " << std::endl;
394 return 0;
395 }
396
397 assert(fMinimizer.get());
398
399 // getting a 2D contour so ndf = 2
400 double cont_level = TMath::ChisquareQuantile(ConfidenceLevel(),2); // level for -2log LR
401 cont_level = cont_level/2; // since we are using -log LR
402 fMinimizer->SetErrorDef(cont_level);
403
404 unsigned int ncp = npoints;
405 unsigned int ivarX = ix;
406 unsigned int ivarY = iy;
407 coutI(Minimization) << "LikelihoodInterval - Finding the contour of " << paramX.GetName() << " ( " << ivarX << " ) and " << paramY.GetName() << " ( " << ivarY << " ) " << std::endl;
408 ret = fMinimizer->Contour(ivarX, ivarY, ncp, x, y );
409 if (!ret) {
410 coutE(Minimization) << "LikelihoodInterval - Error finding contour for parameters " << paramX.GetName() << " and " << paramY.GetName() << std::endl;
411 return 0;
412 }
413 if (int(ncp) < npoints) {
414 coutW(Minimization) << "LikelihoodInterval -Warning - Less points calculated in contours np = " << ncp << " / " << npoints << std::endl;
415 }
416
417 return ncp;
418 }
#define coutI(a)
Definition: RooMsgService.h:34
#define ccoutE(a)
Definition: RooMsgService.h:45
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define ccoutW(a)
Definition: RooMsgService.h:44
#define ccoutI(a)
Definition: RooMsgService.h:42
#define ClassImp(name)
Definition: Rtypes.h:375
char name[80]
Definition: TGX11.cxx:110
static ROOT::Math::Minimizer * CreateMinimizer(const std::string &minimizerType="", const std::string &algoType="")
static method to create the corresponding Minimizer given the string Supported Minimizers types are: ...
Definition: Factory.cxx:63
static const std::string & DefaultMinimizerType()
Template class to wrap any C++ callable object implementing operator() (const double * x) in a multi-...
RooArgSet * getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2083
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
Int_t getSize() const
Return the number of elements in the collection.
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:104
static void setHideOffset(bool flag)
Definition: RooAbsReal.cxx:115
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
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:22
RooAbsReal & nll()
Definition: RooProfileLL.h:34
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
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
double ConfidenceLevel() const override
return confidence level
double UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
Int_t GetContourPoints(const RooRealVar &paramX, const RooRealVar &paramY, double *x, double *y, Int_t npoints=30)
return the 2D-contour points for the given subset of parameters by default make the contour using 30 ...
void ResetLimits()
reset the cached limit values
bool CreateMinimizer()
internal function to create the minimizer for finding the contours
RooArgSet * GetParameters() const override
return a cloned list of parameters of interest. User manages the return object
RooArgSet * fBestFitParams
snapshot of the model parameters with best fit value (managed internally)
RooArgSet fParameters
parameters of interest for this interval
~LikelihoodInterval() override
destructor
double LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
std::shared_ptr< RooFunctor > fFunctor
! transient pointer to functor class used by the minimizer
bool FindLimits(const RooRealVar &param, double &lower, double &upper)
find both lower and upper interval boundaries for a given parameter return false if the bounds have n...
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
std::map< std::string, double > fLowerLimits
map with cached lower bound values
bool IsInInterval(const RooArgSet &) const override
check if given point is in the interval
bool CheckParameters(const RooArgSet &) const override
check if parameters are correct (i.e. they are the POI of this interval)
RooAbsReal * fLikelihoodRatio
likelihood ratio function used to make contours (managed internally)
std::map< std::string, double > fUpperLimits
map with cached upper bound values
std::shared_ptr< ROOT::Math::Minimizer > fMinimizer
! transient pointer to minimizer class used to find limits and contour
LikelihoodInterval(const char *name=nullptr)
default constructor
std::shared_ptr< ROOT::Math::IMultiGenFunction > fMinFunc
! transient pointer to the minimization function
static bool UseStaticMinuit(bool on=true)
static function to switch on/off usage of static global TMinuit instance (gMinuit) By default it is u...
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:60
@ Minimization
Definition: RooGlobalFunc.h:62
@ InputArguments
Definition: RooGlobalFunc.h:63
Namespace for the RooStats classes.
Definition: Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:63
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...
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition: TMath.cxx:637
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition: TMath.cxx:2189