Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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/*
69// for debugging
70#include "RooNLLVar.h"
71#include "RooDataSet.h"
72#include "RooAbsData.h"
73*/
74
76
77using namespace RooStats;
78using namespace std;
79
80
81////////////////////////////////////////////////////////////////////////////////
82/// Default constructor with name and title
83
85 ConfInterval(name), fBestFitParams(0), fLikelihoodRatio(0), fConfidenceLevel(0.95)
86{
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// Alternate constructor taking a pointer to the profile likelihood ratio, parameter of interest and
91/// optionally a snapshot of best parameter of interest for interval
92
93LikelihoodInterval::LikelihoodInterval(const char* name, RooAbsReal* lr, const RooArgSet* params, RooArgSet * bestParams) :
95 fParameters(*params),
96 fBestFitParams(bestParams),
97 fLikelihoodRatio(lr),
98 fConfidenceLevel(0.95)
99{
100}
101
102
103////////////////////////////////////////////////////////////////////////////////
104/// Destructor
105
107{
108 if (fBestFitParams) delete fBestFitParams;
110}
111
112
113////////////////////////////////////////////////////////////////////////////////
114/// This is the main method to satisfy the RooStats::ConfInterval interface.
115/// It returns true if the parameter point is in the interval.
116
118{
121 // Method to determine if a parameter point is in the interval
122 if( !this->CheckParameters(parameterPoint) ) {
123 std::cout << "parameters don't match" << std::endl;
125 return false;
126 }
127
128 // make sure likelihood ratio is set
129 if(!fLikelihoodRatio) {
130 std::cout << "likelihood ratio not set" << std::endl;
132 return false;
133 }
134
135
136
137 // set parameters
138 SetParameters(&parameterPoint, fLikelihoodRatio->getVariables() );
139
140
141 // evaluate likelihood ratio, see if it's bigger than threshold
142 if (fLikelihoodRatio->getVal()<0){
143 std::cout << "The likelihood ratio is < 0, indicates a bad minimum or numerical precision problems. Will return true" << std::endl;
145 return true;
146 }
147
148
149 // here we use Wilks' theorem.
150 if ( TMath::Prob( 2* fLikelihoodRatio->getVal(), parameterPoint.getSize()) < (1.-fConfidenceLevel) ){
152 return false;
153 }
154
155
157
158 return true;
159
160}
161
162////////////////////////////////////////////////////////////////////////////////
163/// returns list of parameters
164
166{
167 return new RooArgSet(fParameters);
168}
169
170////////////////////////////////////////////////////////////////////////////////
171/// check that the parameters are correct
172
174{
175 if (parameterPoint.getSize() != fParameters.getSize() ) {
176 std::cout << "size is wrong, parameters don't match" << std::endl;
177 return false;
178 }
179 if ( ! parameterPoint.equals( fParameters ) ) {
180 std::cout << "size is ok, but parameters don't match" << std::endl;
181 return false;
182 }
183 return true;
184}
185
186
187
188////////////////////////////////////////////////////////////////////////////////
189/// Compute lower limit, check first if limit has been computed
190/// status is a boolean flag which will b set to false in case of error
191/// and is true if calculation is successful
192/// in case of error return also a lower limit value of zero
193
195{
196 double lower = 0;
197 double upper = 0;
198 status = FindLimits(param, lower, upper);
199 return lower;
200}
201
202////////////////////////////////////////////////////////////////////////////////
203/// Compute upper limit, check first if limit has been computed
204/// status is a boolean flag which will b set to false in case of error
205/// and is true if calculation is successful
206/// in case of error return also a lower limit value of zero
207
209{
210 double lower = 0;
211 double upper = 0;
212 status = FindLimits(param, lower, upper);
213 return upper;
214}
215
216
218 // reset map with cached limits - called every time the test size or CL has been changed
219 fLowerLimits.clear();
220 fUpperLimits.clear();
221}
222
223
225 // internal function to create minimizer object needed to find contours or interval limits
226 // (running MINOS).
227 // Minimizer must be Minuit or Minuit2
228
229 RooProfileLL * profilell = dynamic_cast<RooProfileLL*>(fLikelihoodRatio);
230 if (!profilell) return false;
231
232 RooAbsReal & nll = profilell->nll();
233 // bind the nll function in the right interface for the Minimizer class
234 // as a function of only the parameters (poi + nuisance parameters)
235
236 RooArgSet * partmp = profilell->getVariables();
237 // need to remove constant parameters
239
240 RooArgList params(*partmp);
241 delete partmp;
242
243 // need to restore values and errors for POI
244 if (fBestFitParams) {
245 for (int i = 0; i < params.getSize(); ++i) {
246 RooRealVar & par = (RooRealVar &) params[i];
247 RooRealVar * fitPar = (RooRealVar *) (fBestFitParams->find(par.GetName() ) );
248 if (fitPar) {
249 par.setVal( fitPar->getVal() );
250 par.setError( fitPar->getError() );
251 }
252 }
253 }
254
255 const auto& config = GetGlobalRooStatsConfig();
256
257 // now do binding of NLL with a functor for Minimizer
258 if (config.useLikelihoodOffset) {
259 ccoutI(InputArguments) << "LikelihoodInterval: using nll offset - set all RooAbsReal to hide the offset " << std::endl;
260 RooAbsReal::setHideOffset(kFALSE); // need to keep this false
261 }
262 fFunctor = std::make_shared<RooFunctor>(nll, RooArgSet(), params);
263
265 std::transform(minimType.begin(), minimType.end(), minimType.begin(), (int(*)(int)) tolower );
266 *minimType.begin() = toupper( *minimType.begin());
267
268 if (minimType != "Minuit" && minimType != "Minuit2") {
269 ccoutE(InputArguments) << minimType << " is wrong type of minimizer for getting interval limits or contours - must use Minuit or Minuit2" << std::endl;
270 return false;
271 }
272 // do not use static instance of TMInuit which could interfere with RooFit
273 if (minimType == "Minuit") TMinuitMinimizer::UseStaticMinuit(false);
274 // create minimizer class
275 fMinimizer = std::shared_ptr<ROOT::Math::Minimizer>(ROOT::Math::Factory::CreateMinimizer(minimType, "Migrad"));
276
277 if (!fMinimizer.get()) return false;
278
279 fMinFunc = std::static_pointer_cast<ROOT::Math::IMultiGenFunction>(
281 fMinimizer->SetFunction(*fMinFunc);
282
283 // set minimizer parameters
284 assert( params.getSize() == int(fMinFunc->NDim()) );
285
286 for (unsigned int i = 0; i < fMinFunc->NDim(); ++i) {
287 RooRealVar & v = (RooRealVar &) params[i];
288 fMinimizer->SetLimitedVariable( i, v.GetName(), v.getVal(), v.getError(), v.getMin(), v.getMax() );
289 }
290 // for finding the contour need to find first global minimum
291 bool iret = fMinimizer->Minimize();
292 if (!iret || fMinimizer->X() == 0) {
293 ccoutE(Minimization) << "Error: Minimization failed " << std::endl;
294 return false;
295 }
296
297 //std::cout << "print minimizer result..........." << std::endl;
298 //fMinimizer->PrintResults();
299
300 return true;
301}
302
303bool LikelihoodInterval::FindLimits(const RooRealVar & param, double &lower, double & upper)
304{
305 // Method to find both lower and upper limits using MINOS
306 // If cached values exist (limits have been already found) return them in that case
307 // check first if limit has been computed
308 // otherwise compute limit using MINOS
309 // in case of failure lower and upper will maintain previous value (will not be modified)
310
311 std::map<std::string, double>::const_iterator itrl = fLowerLimits.find(param.GetName());
312 std::map<std::string, double>::const_iterator itru = fUpperLimits.find(param.GetName());
313 if ( itrl != fLowerLimits.end() && itru != fUpperLimits.end() ) {
314 lower = itrl->second;
315 upper = itru->second;
316 return true;
317 }
318
319
322 RooArgList params(*partmp);
323 delete partmp;
324 int ix = params.index(&param);
325 if (ix < 0 ) {
326 ccoutE(InputArguments) << "Error - invalid parameter " << param.GetName() << " specified for finding the interval limits " << std::endl;
327 return false;
328 }
329
330 bool ret = true;
331 if (!fMinimizer.get()) ret = CreateMinimizer();
332 if (!ret) {
333 ccoutE(Eval) << "Error returned from minimization of likelihood function - cannot find interval limits " << std::endl;
334 return false;
335 }
336
337 assert(fMinimizer.get());
338
339 // getting a 1D interval so ndf = 1
340 double err_level = TMath::ChisquareQuantile(ConfidenceLevel(),1); // level for -2log LR
341 err_level = err_level/2; // since we are using -log LR
342 fMinimizer->SetErrorDef(err_level);
343
344 unsigned int ivarX = ix;
345
346 double elow = 0;
347 double eup = 0;
348 ret = fMinimizer->GetMinosError(ivarX, elow, eup );
349 if (!ret) {
350 ccoutE(Minimization) << "Error running Minos for parameter " << param.GetName() << std::endl;
351 return false;
352 }
353
354 // WHEN error is zero normally is at limit
355 if (elow == 0) {
356 lower = param.getMin();
357 ccoutW(Minimization) << "Warning: lower value for " << param.GetName() << " is at limit " << lower << std::endl;
358 }
359 else
360 lower = fMinimizer->X()[ivarX] + elow; // elow is negative
361
362 if (eup == 0) {
363 ccoutW(Minimization) << "Warning: upper value for " << param.GetName() << " is at limit " << upper << std::endl;
364 upper = param.getMax();
365 }
366 else
367 upper = fMinimizer->X()[ivarX] + eup;
368
369 // store limits in the map
370 // minos return error limit = minValue +/- error
371 fLowerLimits[param.GetName()] = lower;
372 fUpperLimits[param.GetName()] = upper;
373
374 return true;
375}
376
377
379 // use Minuit to find the contour of the likelihood function at the desired CL
380
381 // check the parameters
382 // variable index in minimizer
383 // is index in the RooArgList obtained from the profileLL variables
386 RooArgList params(*partmp);
387 delete partmp;
388 int ix = params.index(&paramX);
389 int iy = params.index(&paramY);
390 if (ix < 0 || iy < 0) {
391 coutE(InputArguments) << "LikelihoodInterval - Error - invalid parameters specified for finding the contours; parX = " << paramX.GetName()
392 << " parY = " << paramY.GetName() << std::endl;
393 return 0;
394 }
395
396 bool ret = true;
397 if (!fMinimizer.get()) ret = CreateMinimizer();
398 if (!ret) {
399 coutE(Eval) << "LikelihoodInterval - Error returned creating minimizer for likelihood function - cannot find contour points " << std::endl;
400 return 0;
401 }
402
403 assert(fMinimizer.get());
404
405 // getting a 2D contour so ndf = 2
406 double cont_level = TMath::ChisquareQuantile(ConfidenceLevel(),2); // level for -2log LR
407 cont_level = cont_level/2; // since we are using -log LR
408 fMinimizer->SetErrorDef(cont_level);
409
410 unsigned int ncp = npoints;
411 unsigned int ivarX = ix;
412 unsigned int ivarY = iy;
413 coutI(Minimization) << "LikelihoodInterval - Finding the contour of " << paramX.GetName() << " ( " << ivarX << " ) and " << paramY.GetName() << " ( " << ivarY << " ) " << std::endl;
414 ret = fMinimizer->Contour(ivarX, ivarY, ncp, x, y );
415 if (!ret) {
416 coutE(Minimization) << "LikelihoodInterval - Error finding contour for parameters " << paramX.GetName() << " and " << paramY.GetName() << std::endl;
417 return 0;
418 }
419 if (int(ncp) < npoints) {
420 coutW(Minimization) << "LikelihoodInterval -Warning - Less points calculated in contours np = " << ncp << " / " << npoints << std::endl;
421 }
422
423 return ncp;
424 }
#define coutI(a)
#define ccoutE(a)
#define coutW(a)
#define coutE(a)
#define ccoutW(a)
#define ccoutI(a)
const Bool_t kFALSE
Definition RtypesCore.h:101
#define ClassImp(name)
Definition Rtypes.h:364
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 corrisponding 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_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Int_t getSize() const
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:64
static void setHideOffset(Bool_t flag)
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:94
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:35
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...
RooAbsReal & nll()
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
void setError(Double_t value)
Definition RooRealVar.h:64
Double_t getError() const
Definition RooRealVar.h:62
virtual void setVal(Double_t value)
Set value of variable to 'value'.
ConfInterval is an interface class for a generic interval in the RooStats framework.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual RooArgSet * GetParameters() const
return a cloned list of parameters of interest. User manages the return object
Double_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
void ResetLimits()
reset the cached limit values
bool CreateMinimizer()
internal function to create the minimizer for finding the contours
Int_t GetContourPoints(const RooRealVar &paramX, const RooRealVar &paramY, Double_t *x, Double_t *y, Int_t npoints=30)
return the 2D-contour points for the given subset of parameters by default make the contour using 30 ...
Bool_t CheckParameters(const RooArgSet &) const
check if parameters are correct (i.e. they are the POI of this interval)
Double_t fConfidenceLevel
likelihood ratio function used to make contours (managed internally)
RooArgSet * fBestFitParams
parameters of interest for this interval
virtual Double_t ConfidenceLevel() const
return confidence level
std::shared_ptr< RooFunctor > fFunctor
transient pointer to minimizer class used to find limits and contour
Bool_t 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_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
std::map< std::string, double > fLowerLimits
Requested confidence level (eg. 0.95 for 95% CL)
LikelihoodInterval(const char *name=0)
default constructor
RooAbsReal * fLikelihoodRatio
snapshot of the model parameters with best fit value (managed internally)
std::map< std::string, double > fUpperLimits
map with cached lower bound values
std::shared_ptr< ROOT::Math::Minimizer > fMinimizer
map with cached upper bound values
virtual ~LikelihoodInterval()
destructor
virtual Bool_t IsInInterval(const RooArgSet &) const
check if given point is in the interval
std::shared_ptr< ROOT::Math::IMultiGenFunction > fMinFunc
transient pointer to functor class used by the minimizer
static bool UseStaticMinuit(bool on=true)
static function to switch on/off usage of static global TMinuit instance (gMinuit) By default it is u...
virtual const char * GetName() const
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.
Namespace for the RooStats classes.
Definition Asimov.h:19
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
void RemoveConstantParameters(RooArgSet *set)
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:614
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition TMath.cxx:2157