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