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;
72
73
74////////////////////////////////////////////////////////////////////////////////
75/// Default constructor with name and title
76
78 ConfInterval(name), fBestFitParams(nullptr), fLikelihoodRatio(nullptr), fConfidenceLevel(0.95)
79{
80}
81
82////////////////////////////////////////////////////////////////////////////////
83/// Alternate constructor taking a pointer to the profile likelihood ratio, parameter of interest and
84/// optionally a snapshot of best parameter of interest for interval
85
86LikelihoodInterval::LikelihoodInterval(const char* name, RooAbsReal* lr, const RooArgSet* params, RooArgSet * bestParams) :
88 fParameters(*params),
89 fBestFitParams(bestParams),
90 fLikelihoodRatio(lr),
91 fConfidenceLevel(0.95)
92{
93}
94
95
96////////////////////////////////////////////////////////////////////////////////
97/// Destructor
98
100{
101 if (fBestFitParams) delete fBestFitParams;
103}
104
105
106////////////////////////////////////////////////////////////////////////////////
107/// This is the main method to satisfy the RooStats::ConfInterval interface.
108/// It returns true if the parameter point is in the interval.
109
110bool LikelihoodInterval::IsInInterval(const RooArgSet &parameterPoint) const
111{
114 // Method to determine if a parameter point is in the interval
115 if( !this->CheckParameters(parameterPoint) ) {
116 std::cout << "parameters don't match" << std::endl;
118 return false;
119 }
120
121 // make sure likelihood ratio is set
122 if(!fLikelihoodRatio) {
123 std::cout << "likelihood ratio not set" << std::endl;
125 return false;
126 }
127
128
129
130 // set parameters
131 SetParameters(&parameterPoint, std::unique_ptr<RooArgSet>{fLikelihoodRatio->getVariables()}.get());
132
133
134 // evaluate likelihood ratio, see if it's bigger than threshold
135 if (fLikelihoodRatio->getVal()<0){
136 std::cout << "The likelihood ratio is < 0, indicates a bad minimum or numerical precision problems. Will return true" << std::endl;
138 return true;
139 }
140
141
142 // here we use Wilks' theorem.
143 if ( TMath::Prob( 2* fLikelihoodRatio->getVal(), parameterPoint.size()) < (1.-fConfidenceLevel) ){
145 return false;
146 }
147
148
150
151 return true;
152
153}
154
155////////////////////////////////////////////////////////////////////////////////
156/// returns list of parameters
157
159{
160 return new RooArgSet(fParameters);
161}
162
163////////////////////////////////////////////////////////////////////////////////
164/// check that the parameters are correct
165
166bool LikelihoodInterval::CheckParameters(const RooArgSet &parameterPoint) const
167{
168 if (parameterPoint.size() != fParameters.size() ) {
169 std::cout << "size is wrong, parameters don't match" << std::endl;
170 return false;
171 }
172 if ( ! parameterPoint.equals( fParameters ) ) {
173 std::cout << "size is ok, but parameters don't match" << std::endl;
174 return false;
175 }
176 return true;
177}
178
179
180
181////////////////////////////////////////////////////////////////////////////////
182/// Compute lower limit, check first if limit has been computed
183/// status is a boolean flag which will b set to false in case of error
184/// and is true if calculation is successful
185/// in case of error return also a lower limit value of zero
186
187double LikelihoodInterval::LowerLimit(const RooRealVar& param, bool & status)
188{
189 double lower = 0;
190 double upper = 0;
191 status = FindLimits(param, lower, upper);
192 return lower;
193}
194
195////////////////////////////////////////////////////////////////////////////////
196/// Compute upper limit, check first if limit has been computed
197/// status is a boolean flag which will b set to false in case of error
198/// and is true if calculation is successful
199/// in case of error return also a lower limit value of zero
200
201double LikelihoodInterval::UpperLimit(const RooRealVar& param, bool & status)
202{
203 double lower = 0;
204 double upper = 0;
205 status = FindLimits(param, lower, upper);
206 return upper;
207}
208
209
211 // reset map with cached limits - called every time the test size or CL has been changed
212 fLowerLimits.clear();
213 fUpperLimits.clear();
214}
215
216
218 // internal function to create minimizer object needed to find contours or interval limits
219 // (running MINOS).
220 // Minimizer must be Minuit or Minuit2
221
222 RooProfileLL * profilell = dynamic_cast<RooProfileLL*>(fLikelihoodRatio);
223 if (!profilell) return false;
224
225 RooAbsReal & nll = profilell->nll();
226 // bind the nll function in the right interface for the Minimizer class
227 // as a function of only the parameters (poi + nuisance parameters)
228
229 std::unique_ptr<RooArgSet> partmp{profilell->getVariables()};
230 // need to remove constant parameters
231 RemoveConstantParameters(&*partmp);
232
233 RooArgList params(*partmp);
234
235 // need to restore values and errors for POI
236 if (fBestFitParams) {
237 for (std::size_t i = 0; i < params.size(); ++i) {
238 RooRealVar & par = static_cast<RooRealVar &>( params[i]);
239 RooRealVar * fitPar = static_cast<RooRealVar *> (fBestFitParams->find(par.GetName() ) );
240 if (fitPar) {
241 par.setVal( fitPar->getVal() );
242 par.setError( fitPar->getError() );
243 }
244 }
245 }
246
247 const auto& config = GetGlobalRooStatsConfig();
248
249 // now do binding of NLL with a functor for Minimizer
250 if (config.useLikelihoodOffset) {
251 ccoutI(InputArguments) << "LikelihoodInterval: using nll offset - set all RooAbsReal to hide the offset " << std::endl;
252 RooAbsReal::setHideOffset(false); // need to keep this false
253 }
254 fFunctor = std::make_shared<RooFunctor>(nll, RooArgSet(), params);
255
257 std::transform(minimType.begin(), minimType.end(), minimType.begin(), (int(*)(int)) tolower );
258 *minimType.begin() = toupper( *minimType.begin());
259
260 if (minimType != "Minuit" && minimType != "Minuit2") {
261 ccoutE(InputArguments) << minimType << " is wrong type of minimizer for getting interval limits or contours - must use Minuit or Minuit2" << std::endl;
262 return false;
263 }
264 // do not use static instance of TMInuit which could interfere with RooFit
265 if (minimType == "Minuit") TMinuitMinimizer::UseStaticMinuit(false);
266 // create minimizer class
267 fMinimizer = std::shared_ptr<ROOT::Math::Minimizer>(ROOT::Math::Factory::CreateMinimizer(minimType, "Migrad"));
268
269 if (!fMinimizer.get()) return false;
270
271 fMinFunc = std::static_pointer_cast<ROOT::Math::IMultiGenFunction>(
273 fMinimizer->SetFunction(*fMinFunc);
274
275 // set minimizer parameters
276 assert(params.size() == fMinFunc->NDim());
277
278 for (unsigned int i = 0; i < fMinFunc->NDim(); ++i) {
279 RooRealVar & v = static_cast<RooRealVar &>( params[i]);
280 fMinimizer->SetLimitedVariable( i, v.GetName(), v.getVal(), v.getError(), v.getMin(), v.getMax() );
281 }
282 // for finding the contour need to find first global minimum
283 bool iret = fMinimizer->Minimize();
284 if (!iret || fMinimizer->X() == nullptr) {
285 ccoutE(Minimization) << "Error: Minimization failed " << std::endl;
286 return false;
287 }
288
289 //std::cout << "print minimizer result..........." << std::endl;
290 //fMinimizer->PrintResults();
291
292 return true;
293}
294
295bool LikelihoodInterval::FindLimits(const RooRealVar & param, double &lower, double & upper)
296{
297 // Method to find both lower and upper limits using MINOS
298 // If cached values exist (limits have been already found) return them in that case
299 // check first if limit has been computed
300 // otherwise compute limit using MINOS
301 // in case of failure lower and upper will maintain previous value (will not be modified)
302
303 std::map<std::string, double>::const_iterator itrl = fLowerLimits.find(param.GetName());
304 std::map<std::string, double>::const_iterator itru = fUpperLimits.find(param.GetName());
305 if ( itrl != fLowerLimits.end() && itru != fUpperLimits.end() ) {
306 lower = itrl->second;
307 upper = itru->second;
308 return true;
309 }
310
311
312 std::unique_ptr<RooArgSet> partmp{fLikelihoodRatio->getVariables()};
313 RemoveConstantParameters(&*partmp);
314 RooArgList params(*partmp);
315 int iX = params.index(&param);
316 if (iX < 0 ) {
317 ccoutE(InputArguments) << "Error - invalid parameter " << param.GetName() << " specified for finding the interval limits " << std::endl;
318 return false;
319 }
320
321 bool ret = true;
322 if (!fMinimizer.get()) ret = CreateMinimizer();
323 if (!ret) {
324 ccoutE(Eval) << "Error returned from minimization of likelihood function - cannot find interval limits " << std::endl;
325 return false;
326 }
327
328 assert(fMinimizer.get());
329
330 // getting a 1D interval so ndf = 1
331 double err_level = TMath::ChisquareQuantile(ConfidenceLevel(),1); // level for -2log LR
332 err_level = err_level/2; // since we are using -log LR
333 fMinimizer->SetErrorDef(err_level);
334
335 unsigned int ivarX = iX;
336
337 double elow = 0;
338 double eup = 0;
339 ret = fMinimizer->GetMinosError(ivarX, elow, eup );
340
341 // WHEN error is zero normally is at limit
342 if (elow == 0) {
343 lower = param.getMin();
344 ccoutW(Minimization) << "Warning: lower value for " << param.GetName() << " is at limit " << lower << std::endl;
345 }
346 else
347 lower = fMinimizer->X()[ivarX] + elow; // elow is negative
348
349 if (eup == 0) {
350 ccoutW(Minimization) << "Warning: upper value for " << param.GetName() << " is at limit " << upper << std::endl;
351 upper = param.getMax();
352 }
353 else
354 upper = fMinimizer->X()[ivarX] + eup;
355
356 // store limits in the map
357 // minos return error limit = minValue +/- error
358 fLowerLimits[param.GetName()] = lower;
359 fUpperLimits[param.GetName()] = upper;
360
361 return true;
362}
363
364
365Int_t LikelihoodInterval::GetContourPoints(const RooRealVar & paramX, const RooRealVar & paramY, double * x, double *y, Int_t npoints ) {
366 // use Minuit to find the contour of the likelihood function at the desired CL
367
368 // check the parameters
369 // variable index in minimizer
370 // is index in the RooArgList obtained from the profileLL variables
371 std::unique_ptr<RooArgSet> partmp{fLikelihoodRatio->getVariables()};
372 RemoveConstantParameters(&*partmp);
373 RooArgList params(*partmp);
374 int iX = params.index(&paramX);
375 int iY = params.index(&paramY);
376 if (iX < 0 || iY < 0) {
377 coutE(InputArguments) << "LikelihoodInterval - Error - invalid parameters specified for finding the contours; parX = " << paramX.GetName()
378 << " parY = " << paramY.GetName() << std::endl;
379 return 0;
380 }
381
382 bool ret = true;
383 if (!fMinimizer.get()) ret = CreateMinimizer();
384 if (!ret) {
385 coutE(Eval) << "LikelihoodInterval - Error returned creating minimizer for likelihood function - cannot find contour points " << std::endl;
386 return 0;
387 }
388
389 assert(fMinimizer.get());
390
391 // getting a 2D contour so ndf = 2
392 double cont_level = TMath::ChisquareQuantile(ConfidenceLevel(),2); // level for -2log LR
393 cont_level = cont_level/2; // since we are using -log LR
394 fMinimizer->SetErrorDef(cont_level);
395
396 unsigned int ncp = npoints;
397 unsigned int ivarX = iX;
398 unsigned int ivarY = iY;
399 coutI(Minimization) << "LikelihoodInterval - Finding the contour of " << paramX.GetName() << " ( " << ivarX << " ) and " << paramY.GetName() << " ( " << ivarY << " ) " << std::endl;
400 ret = fMinimizer->Contour(ivarX, ivarY, ncp, x, y );
401 if (!ret) {
402 coutE(Minimization) << "LikelihoodInterval - Error finding contour for parameters " << paramX.GetName() << " and " << paramY.GetName() << std::endl;
403 return 0;
404 }
405 if (int(ncp) < npoints) {
406 coutW(Minimization) << "LikelihoodInterval -Warning - Less points calculated in contours np = " << ncp << " / " << npoints << std::endl;
407 }
408
409 return ncp;
410 }
#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:2193