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 <cctype> // 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
79 ConfInterval(name), fBestFitParams(nullptr), fLikelihoodRatio(nullptr), 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, std::unique_ptr<RooArgSet>{fLikelihoodRatio->getVariables()}.get());
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.size()) < (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.size() != fParameters.size() ) {
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 std::unique_ptr<RooArgSet> partmp{profilell->getVariables()};
231 // need to remove constant parameters
232 RemoveConstantParameters(&*partmp);
233
234 RooArgList params(*partmp);
235
236 // need to restore values and errors for POI
237 if (fBestFitParams) {
238 for (std::size_t i = 0; i < params.size(); ++i) {
239 RooRealVar & par = static_cast<RooRealVar &>( params[i]);
240 RooRealVar * fitPar = static_cast<RooRealVar *> (fBestFitParams->find(par.GetName() ) );
241 if (fitPar) {
242 par.setVal( fitPar->getVal() );
243 par.setError( fitPar->getError() );
244 }
245 }
246 }
247
248 const auto& config = GetGlobalRooStatsConfig();
249
250 // now do binding of NLL with a functor for Minimizer
251 if (config.useLikelihoodOffset) {
252 ccoutI(InputArguments) << "LikelihoodInterval: using nll offset - set all RooAbsReal to hide the offset " << std::endl;
253 RooAbsReal::setHideOffset(false); // need to keep this false
254 }
255 fFunctor = std::make_shared<RooFunctor>(nll, RooArgSet(), params);
256
258 std::transform(minimType.begin(), minimType.end(), minimType.begin(), (int(*)(int)) tolower );
259 *minimType.begin() = toupper( *minimType.begin());
260
261 if (minimType != "Minuit" && minimType != "Minuit2") {
262 ccoutE(InputArguments) << minimType << " is wrong type of minimizer for getting interval limits or contours - must use Minuit or Minuit2" << std::endl;
263 return false;
264 }
265 // do not use static instance of TMInuit which could interfere with RooFit
266 if (minimType == "Minuit") TMinuitMinimizer::UseStaticMinuit(false);
267 // create minimizer class
268 fMinimizer = std::shared_ptr<ROOT::Math::Minimizer>(ROOT::Math::Factory::CreateMinimizer(minimType, "Migrad"));
269
270 if (!fMinimizer.get()) return false;
271
272 fMinFunc = std::static_pointer_cast<ROOT::Math::IMultiGenFunction>(
274 fMinimizer->SetFunction(*fMinFunc);
275
276 // set minimizer parameters
277 assert(params.size() == fMinFunc->NDim());
278
279 for (unsigned int i = 0; i < fMinFunc->NDim(); ++i) {
280 RooRealVar & v = static_cast<RooRealVar &>( params[i]);
281 fMinimizer->SetLimitedVariable( i, v.GetName(), v.getVal(), v.getError(), v.getMin(), v.getMax() );
282 }
283 // for finding the contour need to find first global minimum
284 bool iret = fMinimizer->Minimize();
285 if (!iret || fMinimizer->X() == nullptr) {
286 ccoutE(Minimization) << "Error: Minimization failed " << std::endl;
287 return false;
288 }
289
290 //std::cout << "print minimizer result..........." << std::endl;
291 //fMinimizer->PrintResults();
292
293 return true;
294}
295
296bool LikelihoodInterval::FindLimits(const RooRealVar & param, double &lower, double & upper)
297{
298 // Method to find both lower and upper limits using MINOS
299 // If cached values exist (limits have been already found) return them in that case
300 // check first if limit has been computed
301 // otherwise compute limit using MINOS
302 // in case of failure lower and upper will maintain previous value (will not be modified)
303
304 std::map<std::string, double>::const_iterator itrl = fLowerLimits.find(param.GetName());
305 std::map<std::string, double>::const_iterator itru = fUpperLimits.find(param.GetName());
306 if ( itrl != fLowerLimits.end() && itru != fUpperLimits.end() ) {
307 lower = itrl->second;
308 upper = itru->second;
309 return true;
310 }
311
312
313 std::unique_ptr<RooArgSet> partmp{fLikelihoodRatio->getVariables()};
314 RemoveConstantParameters(&*partmp);
315 RooArgList params(*partmp);
316 int ix = params.index(&param);
317 if (ix < 0 ) {
318 ccoutE(InputArguments) << "Error - invalid parameter " << param.GetName() << " specified for finding the interval limits " << std::endl;
319 return false;
320 }
321
322 bool ret = true;
323 if (!fMinimizer.get()) ret = CreateMinimizer();
324 if (!ret) {
325 ccoutE(Eval) << "Error returned from minimization of likelihood function - cannot find interval limits " << std::endl;
326 return false;
327 }
328
329 assert(fMinimizer.get());
330
331 // getting a 1D interval so ndf = 1
332 double err_level = TMath::ChisquareQuantile(ConfidenceLevel(),1); // level for -2log LR
333 err_level = err_level/2; // since we are using -log LR
334 fMinimizer->SetErrorDef(err_level);
335
336 unsigned int ivarX = ix;
337
338 double elow = 0;
339 double eup = 0;
340 ret = fMinimizer->GetMinosError(ivarX, elow, eup );
341
342 // WHEN error is zero normally is at limit
343 if (elow == 0) {
344 lower = param.getMin();
345 ccoutW(Minimization) << "Warning: lower value for " << param.GetName() << " is at limit " << lower << std::endl;
346 }
347 else
348 lower = fMinimizer->X()[ivarX] + elow; // elow is negative
349
350 if (eup == 0) {
351 ccoutW(Minimization) << "Warning: upper value for " << param.GetName() << " is at limit " << upper << std::endl;
352 upper = param.getMax();
353 }
354 else
355 upper = fMinimizer->X()[ivarX] + eup;
356
357 // store limits in the map
358 // minos return error limit = minValue +/- error
359 fLowerLimits[param.GetName()] = lower;
360 fUpperLimits[param.GetName()] = upper;
361
362 return true;
363}
364
365
366Int_t LikelihoodInterval::GetContourPoints(const RooRealVar & paramX, const RooRealVar & paramY, double * x, double *y, Int_t npoints ) {
367 // use Minuit to find the contour of the likelihood function at the desired CL
368
369 // check the parameters
370 // variable index in minimizer
371 // is index in the RooArgList obtained from the profileLL variables
372 std::unique_ptr<RooArgSet> partmp{fLikelihoodRatio->getVariables()};
373 RemoveConstantParameters(&*partmp);
374 RooArgList params(*partmp);
375 int ix = params.index(&paramX);
376 int iy = params.index(&paramY);
377 if (ix < 0 || iy < 0) {
378 coutE(InputArguments) << "LikelihoodInterval - Error - invalid parameters specified for finding the contours; parX = " << paramX.GetName()
379 << " parY = " << paramY.GetName() << std::endl;
380 return 0;
381 }
382
383 bool ret = true;
384 if (!fMinimizer.get()) ret = CreateMinimizer();
385 if (!ret) {
386 coutE(Eval) << "LikelihoodInterval - Error returned creating minimizer for likelihood function - cannot find contour points " << std::endl;
387 return 0;
388 }
389
390 assert(fMinimizer.get());
391
392 // getting a 2D contour so ndf = 2
393 double cont_level = TMath::ChisquareQuantile(ConfidenceLevel(),2); // level for -2log LR
394 cont_level = cont_level/2; // since we are using -log LR
395 fMinimizer->SetErrorDef(cont_level);
396
397 unsigned int ncp = npoints;
398 unsigned int ivarX = ix;
399 unsigned int ivarY = iy;
400 coutI(Minimization) << "LikelihoodInterval - Finding the contour of " << paramX.GetName() << " ( " << ivarX << " ) and " << paramY.GetName() << " ( " << ivarY << " ) " << std::endl;
401 ret = fMinimizer->Contour(ivarX, ivarY, ncp, x, y );
402 if (!ret) {
403 coutE(Minimization) << "LikelihoodInterval - Error finding contour for parameters " << paramX.GetName() << " and " << paramY.GetName() << std::endl;
404 return 0;
405 }
406 if (int(ncp) < npoints) {
407 coutW(Minimization) << "LikelihoodInterval -Warning - Less points calculated in contours np = " << ncp << " / " << npoints << std::endl;
408 }
409
410 return ncp;
411 }
#define coutI(a)
#define ccoutE(a)
#define coutW(a)
#define coutE(a)
#define ccoutW(a)
#define ccoutI(a)
#define ClassImp(name)
Definition Rtypes.h:377
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-...
RooFit::OwningPtr< RooArgSet > getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expression tree)
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
Storage_t::size_type size() const
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.
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
static void setHideOffset(bool flag)
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
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Implements the profile likelihood estimator for a given likelihood and set of parameters of interest.
RooAbsReal & nll()
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
ConfInterval is an interface class for a generic interval in the RooStats framework.
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.
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:637
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition TMath.cxx:2189