Logo ROOT   master
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 
49 #include "RooStats/RooStatsUtils.h"
50 
51 #include "RooAbsReal.h"
52 #include "RooMsgService.h"
53 
54 #include "Math/WrappedFunction.h"
55 #include "Math/Minimizer.h"
56 #include "Math/Factory.h"
57 #include "Math/MinimizerOptions.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 
77 using namespace RooStats;
78 using 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 
93 LikelihoodInterval::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
238  RemoveConstantParameters(partmp);
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 
264  std::string minimType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();
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>(
280  std::make_shared<ROOT::Math::WrappedMultiFunction<RooFunctor &>>(*fFunctor, fFunctor->nPar()) );
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 
303 bool 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 
321  RemoveConstantParameters(partmp);
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 
378 Int_t LikelihoodInterval::GetContourPoints(const RooRealVar & paramX, const RooRealVar & paramY, Double_t * x, Double_t *y, Int_t npoints ) {
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
385  RemoveConstantParameters(partmp);
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  }
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1910
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooAbsReal * fLikelihoodRatio
snapshot of the model parameters with best fit value (managed internally)
#define coutE(a)
Definition: RooMsgService.h:33
std::map< std::string, double > fLowerLimits
Requested confidence level (eg. 0.95 for 95% CL)
std::shared_ptr< ROOT::Math::Minimizer > fMinimizer
map with cached upper bound values
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
std::map< std::string, double > fUpperLimits
map with cached lower bound values
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:26
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
std::shared_ptr< ROOT::Math::IMultiGenFunction > fMinFunc
transient pointer to functor class used by the minimizer
#define coutI(a)
Definition: RooMsgService.h:30
LikelihoodInterval(const char *name=0)
default constructor
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
Definition: RooStatsUtils.h:65
RooFit::MsgLevel globalKillBelow() const
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
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
virtual RooArgSet * GetParameters() const
return a cloned list of parameters of interest. User manages the return object
virtual Double_t ConfidenceLevel() const
return confidence level
Bool_t equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
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...
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
#define coutW(a)
Definition: RooMsgService.h:32
RooArgSet * fBestFitParams
parameters of interest for this interval
static void setHideOffset(Bool_t flag)
Definition: RooAbsReal.cxx:116
Double_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
#define ccoutI(a)
Definition: RooMsgService.h:38
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
Double_t x[n]
Definition: legend1.C:17
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
virtual Bool_t IsInInterval(const RooArgSet &) const
check if given point is in the interval
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:225
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
Int_t getSize() const
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
#define ccoutE(a)
Definition: RooMsgService.h:41
bool CreateMinimizer()
internal function to create the minimizer for finding the contours
void ResetLimits()
reset the cached limit values
static const std::string & DefaultMinimizerType()
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooAbsReal & nll()
Definition: RooProfileLL.h:39
void setGlobalKillBelow(RooFit::MsgLevel level)
const Bool_t kFALSE
Definition: RtypesCore.h:90
std::shared_ptr< RooFunctor > fFunctor
transient pointer to minimizer class used to find limits and contour
ConfInterval is an interface class for a generic interval in the RooStats framework.
Definition: ConfInterval.h:35
Bool_t CheckParameters(const RooArgSet &) const
check if parameters are correct (i.e. they are the POI of this interval)
Namespace for the RooStats classes.
Definition: Asimov.h:19
virtual ~LikelihoodInterval()
destructor
#define ClassImp(name)
Definition: Rtypes.h:361
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
Double_t y[n]
Definition: legend1.C:17
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 ...
#define ccoutW(a)
Definition: RooMsgService.h:40
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Double_t fConfidenceLevel
likelihood ratio function used to make contours (managed internally)
static bool UseStaticMinuit(bool on=true)
static function to switch on/off usage of static global TMinuit instance (gMinuit) By default it is u...
Double_t getError() const
Definition: RooRealVar.h:62
char name[80]
Definition: TGX11.cxx:109
void setError(Double_t value)
Definition: RooRealVar.h:64