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