Logo ROOT   6.10/09
Reference Guide
ProfileLikelihoodCalculator.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer 28/07/2008
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class RooStats::ProfileLikelihoodCalculator
13  \ingroup Roostats
14 
15 ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator
16 (the interface class for a tools which can produce both RooStats HypoTestResults
17 and ConfIntervals). The tool uses the profile likelihood ratio as a test statistic,
18 and assumes that Wilks' theorem is valid. Wilks' theorem states that -2* log
19 (profile likelihood ratio) is asymptotically distributed as a chi^2 distribution
20 with N-dof, where N is the number of degrees of freedom. Thus, p-values can be
21 constructed and the profile likelihood ratio can be used to construct a
22 LikelihoodInterval. (In the future, this class could be extended to use toy
23 Monte Carlo to calibrate the distribution of the test statistic).
24 
25 Usage: It uses the interface of the CombinedCalculator, so that it can be
26 configured by specifying:
27 
28  - a model common model (eg. a family of specific models which includes both
29  the null and alternate),
30  - a data set,
31  - a set of parameters of interest. The nuisance parameters will be all other
32  parameters of the model
33  - a set of parameters of which specify the null hypothesis (including values
34  and const/non-const status)
35 
36 The interface allows one to pass the model, data, and parameters either directly
37 or via a ModelConfig class. The alternate hypothesis leaves the parameter free
38 to take any value other than those specified by the null hypotesis. There is
39 therefore no need to specify the alternate parameters.
40 
41 After configuring the calculator, one only needs to ask GetHypoTest() (which
42 will return a HypoTestResult pointer) or GetInterval() (which will return an
43 ConfInterval pointer).
44 
45 This calculator can work with both one-dimensional intervals or multi-
46 dimensional ones (contours)
47 
48 Note that for hypothesis tests, it is often better to use the
49 RooStats::AsymptoricCalculator class, which can compute in addition the expected
50 p-value using an Asimov data set.
51 
52 */
53 
55 
56 #include "RooStats/RooStatsUtils.h"
57 
60 
61 #include "RooFitResult.h"
62 #include "RooRealVar.h"
63 #include "RooProfileLL.h"
64 #include "RooNLLVar.h"
65 #include "RooGlobalFunc.h"
66 #include "RooMsgService.h"
67 
68 #include "Math/MinimizerOptions.h"
69 #include "RooMinimizer.h"
70 //#include "RooProdPdf.h"
71 
72 using namespace std;
73 
75 
76 using namespace RooFit;
77 using namespace RooStats;
78 
79 
80 ////////////////////////////////////////////////////////////////////////////////
81 /// default constructor
82 
83 ProfileLikelihoodCalculator::ProfileLikelihoodCalculator() :
84  CombinedCalculator(), fFitResult(0), fGlobalFitDone(false)
85 {
86 }
87 
89  Double_t size, const RooArgSet* nullParams ) :
90  CombinedCalculator(data,pdf, paramsOfInterest, size, nullParams ),
91  fFitResult(0), fGlobalFitDone(false)
92 {
93  // constructor from pdf and parameters
94  // the pdf must contain eventually the nuisance parameters
95 }
96 
98  CombinedCalculator(data, model, size),
99  fFitResult(0), fGlobalFitDone(false)
100 {
101  // construct from a ModelConfig. Assume data model.GetPdf() will provide full description of model including
102  // constraint term on the nuisances parameters
103  assert(model.GetPdf() );
104 }
105 
106 
107 ////////////////////////////////////////////////////////////////////////////////
108 /// destructor
109 /// cannot delete prod pdf because it will delete all the composing pdf's
110 /// if (fOwnPdf) delete fPdf;
111 /// fPdf = 0;
112 
114  if (fFitResult) delete fFitResult;
115 }
116 
118  // reset and clear fit result
119  // to be called when a new model or data are set in the calculator
120  if (fFitResult) delete fFitResult;
121  fFitResult = 0;
122 }
123 
125  // perform a global fit of the likelihood letting with all parameter of interest and
126  // nuisance parameters
127  // keep the list of fitted parameters
128 
129  DoReset();
130  RooAbsPdf * pdf = GetPdf();
131  RooAbsData* data = GetData();
132  if (!data || !pdf ) return 0;
133 
134  // get all non-const parameters
135  RooArgSet* constrainedParams = pdf->getParameters(*data);
136  if (!constrainedParams) return 0;
137  RemoveConstantParameters(constrainedParams);
138 
139 
140  RooAbsReal * nll = pdf->createNLL(*data, CloneData(true), Constrain(*constrainedParams),ConditionalObservables(fConditionalObs), Offset(RooStats::IsNLLOffset() ) );
141 
142  // check if global fit has been already done
143  if (fFitResult && fGlobalFitDone) {
144  delete constrainedParams;
145  return nll;
146  }
147 
148  // calculate MLE
149  oocoutP((TObject*)0,Minimization) << "ProfileLikelihoodCalcultor::DoGLobalFit - find MLE " << std::endl;
150 
151  if (fFitResult) delete fFitResult;
152  fFitResult = DoMinimizeNLL(nll);
153 
154  // print fit result
155  if (fFitResult) {
157 
158  if (fFitResult->status() != 0)
159  oocoutW((TObject*)0,Minimization) << "ProfileLikelihoodCalcultor::DoGlobalFit - Global fit failed - status = " << fFitResult->status() << std::endl;
160  else
161  fGlobalFitDone = true;
162  }
163 
164  delete constrainedParams;
165  return nll;
166 }
167 
169  // Minimizer the given NLL using the default options
170 
171  const char * minimType = ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str();
172  const char * minimAlgo = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str();
174  int level = ROOT::Math::MinimizerOptions::DefaultPrintLevel() -1;// RooFit level starts from -1
176  oocoutP((TObject*)0,Minimization) << "ProfileLikelihoodCalcultor::DoMinimizeNLL - using " << minimType << " / " << minimAlgo << " with strategy " << strategy << std::endl;
177  // do global fit and store fit result for further use
178 
179  RooMinimizer minim(*nll);
180  minim.setStrategy(strategy);
181  minim.setEps(tolerance);
182  minim.setPrintLevel(level);
183  minim.optimizeConst(2); // to optimize likelihood calculations
184 
185  int status = -1;
186  for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
187  status = minim.minimize(minimType,minimAlgo);
188  if (status%1000 == 0) { // ignore erros from Improve
189  break;
190  } else if (tries < maxtries) {
191  cout << " ----> Doing a re-scan first" << endl;
192  minim.minimize(minimType,"Scan");
193  if (tries == 2) {
194  if (strategy == 0 ) {
195  cout << " ----> trying with strategy = 1" << endl;;
196  minim.setStrategy(1);
197  }
198  else
199  tries++; // skip this trial if strategy is already 1
200  }
201  if (tries == 3) {
202  cout << " ----> trying with improve" << endl;;
203  minimType = "Minuit";
204  minimAlgo = "migradimproved";
205  }
206  }
207  }
208 
209  RooFitResult * result = minim.save();
210 
211 
212  return result;
213 }
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 /// Main interface to get a RooStats::ConfInterval.
217 /// It constructs a profile likelihood ratio and uses that to construct a RooStats::LikelihoodInterval.
218 
220 // RooAbsPdf* pdf = fWS->pdf(fPdfName);
221 // RooAbsData* data = fWS->data(fDataName);
222  RooAbsPdf * pdf = GetPdf();
223  RooAbsData* data = GetData();
224  if (!data || !pdf || fPOI.getSize() == 0) return 0;
225 
226  RooArgSet* constrainedParams = pdf->getParameters(*data);
227  RemoveConstantParameters(constrainedParams);
228 
229 
230  /*
231  RooNLLVar* nll = new RooNLLVar("nll","",*pdf,*data, Extended(),Constrain(*constrainedParams));
232  RooProfileLL* profile = new RooProfileLL("pll","",*nll, *fPOI);
233  profile->addOwnedComponents(*nll) ; // to avoid memory leak
234  */
235 
236  // do a global fit cloning the data
237  RooAbsReal * nll = DoGlobalFit();
238  if (!nll) return 0;
239 
240  if (!fFitResult) {
241  delete nll;
242  return 0;
243  }
244 
245  RooAbsReal* profile = nll->createProfile(fPOI);
246  profile->addOwnedComponents(*nll) ; // to avoid memory leak
247 
248  // t.b.f. " RooProfileLL should keep and provide possibility to query on global minimum
249  // set POI to fit value (this will speed up profileLL calculation of global minimum)
250  const RooArgList & fitParams = fFitResult->floatParsFinal();
251  for (int i = 0; i < fitParams.getSize(); ++i) {
252  RooRealVar & fitPar = (RooRealVar &) fitParams[i];
253  RooRealVar * par = (RooRealVar*) fPOI.find( fitPar.GetName() );
254  if (par) {
255  par->setVal( fitPar.getVal() );
256  par->setError( fitPar.getError() );
257  }
258  }
259 
260  // do this so profile will cache inside the absolute minimum and
261  // minimum values of nuisance parameters
262  // (no need to this here)
263  // profile->getVal();
264  //RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
265  // profile->Print();
266 
267  TString name = TString("LikelihoodInterval_");// + TString(GetName() );
268 
269  // make a list of fPOI with fit result values and pass to LikelihoodInterval class
270  // bestPOI is a cloned list of POI only with their best fit values
271  TIter iter = fPOI.createIterator();
272  RooArgSet fitParSet(fitParams);
273  RooArgSet * bestPOI = new RooArgSet();
274  while (RooAbsArg * arg = (RooAbsArg*) iter.Next() ) {
275  RooAbsArg * p = fitParSet.find( arg->GetName() );
276  if (p) bestPOI->addClone(*p);
277  else bestPOI->addClone(*arg);
278  }
279  // fPOI contains the parameter of interest of the PL object
280  // and bestPOI contains a snapshot with the best fit values
281  LikelihoodInterval* interval = new LikelihoodInterval(name, profile, &fPOI, bestPOI);
282  interval->SetConfidenceLevel(1.-fSize);
283  delete constrainedParams;
284  return interval;
285 }
286 
287 ////////////////////////////////////////////////////////////////////////////////
288 /// Main interface to get a HypoTestResult.
289 /// It does two fits:
290 /// the first lets the null parameters float, so it's a maximum likelihood estimate
291 /// the second is to the null (fixing null parameters to their specified values): eg. a conditional maximum likelihood
292 /// the ratio of the likelihood at the conditional MLE to the MLE is the profile likelihood ratio.
293 /// Wilks' theorem is used to get p-values
294 
296 // RooAbsPdf* pdf = fWS->pdf(fPdfName);
297 // RooAbsData* data = fWS->data(fDataName);
298  RooAbsPdf * pdf = GetPdf();
299  RooAbsData* data = GetData();
300 
301 
302  if (!data || !pdf) return 0;
303 
304  if (fNullParams.getSize() == 0) return 0;
305 
306  // make a clone and ordered list since a vector will be associated to keep parameter values
307  // clone the list since first fit will changes the fNullParams values
308  RooArgList poiList;
309  poiList.addClone(fNullParams); // make a clone list
310 
311 
312  // do a global fit
313  RooAbsReal * nll = DoGlobalFit();
314  if (!nll) return 0;
315 
316  if (!fFitResult) {
317  delete nll;
318  return 0;
319  }
320 
321  RooArgSet* constrainedParams = pdf->getParameters(*data);
322  RemoveConstantParameters(constrainedParams);
323 
324  Double_t nLLatMLE = fFitResult->minNll();
325  // in case of using offset need to save offset value
326  Double_t nlloffset = (RooStats::IsNLLOffset() ) ? nll->getVal() - nLLatMLE : 0;
327 
328  // set POI to given values, set constant, calculate conditional MLE
329  std::vector<double> oldValues(poiList.getSize() );
330  for (unsigned int i = 0; i < oldValues.size(); ++i) {
331  RooRealVar * mytarget = (RooRealVar*) constrainedParams->find(poiList[i].GetName());
332  if (mytarget) {
333  oldValues[i] = mytarget->getVal();
334  mytarget->setVal( ( (RooRealVar&) poiList[i] ).getVal() );
335  mytarget->setConstant(kTRUE);
336  }
337  }
338 
339 
340 
341  // perform the fit only if nuisance parameters are available
342  // get nuisance parameters
343  // nuisance parameters are the non const parameters from the likelihood parameters
344  RooArgSet nuisParams(*constrainedParams);
345 
346  // need to remove the parameter of interest
347  RemoveConstantParameters(&nuisParams);
348 
349  // check there are variable parameter in order to do a fit
350  bool existVarParams = false;
351  TIter it = nuisParams.createIterator();
352  RooRealVar * myarg = 0;
353  while ((myarg = (RooRealVar *)it.Next())) {
354  if ( !myarg->isConstant() ) {
355  existVarParams = true;
356  break;
357  }
358  }
359 
360  Double_t nLLatCondMLE = nLLatMLE;
361  if (existVarParams) {
362  oocoutP((TObject*)0,Minimization) << "ProfileLikelihoodCalcultor::GetHypoTest - do conditional fit " << std::endl;
363 
364  RooFitResult * fit2 = DoMinimizeNLL(nll);
365 
366  // print fit result
367  if (fit2) {
368  nLLatCondMLE = fit2->minNll();
370 
371  if (fit2->status() != 0)
372  oocoutW((TObject*)0,Minimization) << "ProfileLikelihoodCalcultor::GetHypotest - Conditional fit failed - status = " << fit2->status() << std::endl;
373  }
374 
375  }
376  else {
377  // get just the likelihood value (no need to do a fit since the likelihood is a constant function)
378  nLLatCondMLE = nll->getVal();
379  // this value contains the offset
380  if (RooStats::IsNLLOffset() ) nLLatCondMLE -= nlloffset;
381  }
382 
383  // Use Wilks' theorem to translate -2 log lambda into a significance/p-value
384  Double_t deltaNLL = std::max( nLLatCondMLE-nLLatMLE, 0.);
385 
386  // get number of free parameter of interest
387  RemoveConstantParameters(poiList);
388  int ndf = poiList.getSize();
389 
390  Double_t pvalue = ROOT::Math::chisquared_cdf_c( 2* deltaNLL, ndf);
391 
392  // in case of one dimension (1 poi) do the one-sided p-value (need to divide by 2)
393  if (ndf == 1) pvalue = 0.5 * pvalue;
394 
395  TString name = TString("ProfileLRHypoTestResult_");// + TString(GetName() );
396  HypoTestResult* htr = new HypoTestResult(name, pvalue, 0 );
397 
398  // restore previous value of poi
399  for (unsigned int i = 0; i < oldValues.size(); ++i) {
400  RooRealVar * mytarget = (RooRealVar*) constrainedParams->find(poiList[i].GetName());
401  if (mytarget) {
402  mytarget->setVal(oldValues[i] );
403  mytarget->setConstant(false);
404  }
405  }
406 
407  delete constrainedParams;
408  delete nll;
409  return htr;
410 
411 }
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:777
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
double par[1]
Definition: unuranDistr.cxx:38
RooCmdArg Offset(Bool_t flag=kTRUE)
TIterator * createIterator(Bool_t dir=kIterForward) const
virtual void printStream(std::ostream &os, Int_t contents, StyleOption style, TString indent="") const
Print description of object on ostream, printing contents set by contents integer, which is interpreted as an OR of &#39;enum ContentsOptions&#39; values and in the style given by &#39;enum StyleOption&#39;.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:108
void optimizeConst(Int_t flag)
If flag is true, perform constant term optimization on function being minimized.
virtual RooAbsArg * addClone(const RooAbsArg &var, Bool_t silent=kFALSE)
Add a clone of the specified argument to list.
virtual LikelihoodInterval * GetInterval() const
Return a likelihood interval.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
RooCmdArg CloneData(Bool_t flag)
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
static RooFitResult * DoMinimizeNLL(RooAbsReal *nll)
virtual HypoTestResult * GetHypoTest() const
Return the hypothesis test result obtained from the likelihood ratio of the maximum likelihood value ...
#define oocoutI(o, a)
Definition: RooMsgService.h:44
HypoTestResult is a base class for results from hypothesis tests.
ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface class f...
virtual Int_t defaultPrintContents(Option_t *opt) const
Configure default contents to be printed.
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of &#39;comps&#39;.
Definition: RooAbsArg.cxx:2282
STL namespace.
ProfileLikelihoodCalculator()
Default constructor (needed for I/O)
void setEps(Double_t eps)
Change MINUIT epsilon.
void setStrategy(Int_t strat)
Change MINUIT strategy to istrat.
virtual StyleOption defaultPrintStyle(Option_t *opt) const
Configure mapping of Print() arguments to RooPrintable print styles.
#define oocoutP(o, a)
Definition: RooMsgService.h:45
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void addClone(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling addOwned() for each element in the source...
Definition: RooArgSet.h:94
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
Int_t getSize() const
static const std::string & DefaultMinimizerType()
void setConstant(Bool_t value=kTRUE)
RooFitResult * save(const char *name=0, const char *title=0)
Save and return a RooFitResult snaphot of current minimizer status.
TObject * Next()
Definition: TCollection.h:153
Double_t minNll() const
Definition: RooFitResult.h:96
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
CombinedCalculator is an interface class for a tools which can produce both RooStats HypoTestResults ...
static const std::string & DefaultMinimizerAlgo()
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:222
#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
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:560
Int_t minimize(const char *type, const char *alg=0)
virtual ~ProfileLikelihoodCalculator()
destructor cannot delete prod pdf because it will delete all the composing pdf&#39;s if (fOwnPdf) delete ...
virtual void SetConfidenceLevel(Double_t cl)
set the confidence level for the interval (e.g 0.682 for a 1-sigma interval)
bool IsNLLOffset()
const char * GetName() const
Returns name of object.
#define oocoutW(o, a)
Definition: RooMsgService.h:46
Mother of all ROOT objects.
Definition: TObject.h:37
Int_t setPrintLevel(Int_t newLevel)
Change the MINUIT internal printing level.
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:62
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
double chisquared_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of the distribution with degrees of freedom (upp...
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:38
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:461
double result[121]
const int strategy
Definition: testNdimFit.cxx:46
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
RooCmdArg ConditionalObservables(const RooArgSet &set)
Double_t getError() const
Definition: RooRealVar.h:53
const Bool_t kTRUE
Definition: RtypesCore.h:91
Bool_t isConstant() const
Definition: RooAbsArg.h:266
RooCmdArg Constrain(const RooArgSet &params)
void setError(Double_t value)
Definition: RooRealVar.h:55
Int_t status() const
Definition: RooFitResult.h:75