Logo ROOT  
Reference Guide
ProfileLikelihoodTestStat.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
3 // Additional Contributions: Giovanni Petrucciani
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::ProfileLikelihoodTestStat
13  \ingroup Roostats
14 
15 ProfileLikelihoodTestStat is an implementation of the TestStatistic interface
16 that calculates the profile likelihood ratio at a particular parameter point
17 given a dataset. It does not constitute a statistical test, for that one may
18 either use:
19 
20  - the ProfileLikelihoodCalculator that relies on asymptotic properties of the
21  Profile Likelihood Ratio
22  - the NeymanConstruction class with this class as a test statistic
23  - the HybridCalculator class with this class as a test statistic
24 
25 
26 */
27 
29 #include "RooFitResult.h"
30 #include "RooPullVar.h"
32 
33 #include "RooProfileLL.h"
34 #include "RooNLLVar.h"
35 #include "RooMsgService.h"
36 #include "RooMinimizer.h"
37 #include "RooArgSet.h"
38 #include "RooDataSet.h"
39 #include "TStopwatch.h"
40 
41 #include "RooStats/RooStatsUtils.h"
42 
43 using namespace std;
44 
46 
47 void RooStats::ProfileLikelihoodTestStat::SetAlwaysReuseNLL(Bool_t flag) { fgAlwaysReuseNll = flag ; }
48 
49 ////////////////////////////////////////////////////////////////////////////////
50 /// internal function to evaluate test statistics
51 /// can do depending on type:
52 /// - type = 0 standard evaluation,
53 /// - type = 1 find only unconditional NLL minimum,
54 /// - type = 2 conditional MLL
55 
57 
58  if( fDetailedOutputEnabled && fDetailedOutput ) {
59  delete fDetailedOutput;
60  fDetailedOutput = 0;
61  }
62  if( fDetailedOutputEnabled && !fDetailedOutput ) {
63  fDetailedOutput = new RooArgSet();
64  }
65 
66  //data.Print("V");
67 
68  TStopwatch tsw;
69  tsw.Start();
70 
71  double initial_mu_value = 0;
72  RooRealVar* firstPOI = dynamic_cast<RooRealVar*>( paramsOfInterest.first());
73  if (firstPOI) initial_mu_value = firstPOI->getVal();
74  //paramsOfInterest.getRealValue(firstPOI->GetName());
75  if (fPrintLevel > 1) {
76  cout << "POIs: " << endl;
77  paramsOfInterest.Print("v");
78  }
79 
82 
83  // simple
84  Bool_t reuse=(fReuseNll || fgAlwaysReuseNll) ;
85 
86  Bool_t created(kFALSE) ;
87  if (!reuse || fNll==0) {
88  RooArgSet* allParams = fPdf->getParameters(data);
90 
91  // need to call constrain for RooSimultaneous until stripDisconnected problem fixed
92  fNll = fPdf->createNLL(data, RooFit::CloneData(kFALSE),RooFit::Constrain(*allParams),
93  RooFit::GlobalObservables(fGlobalObs), RooFit::ConditionalObservables(fConditionalObs), RooFit::Offset(fLOffset));
94 
95  if (fPrintLevel > 0 && fLOffset) cout << "ProfileLikelihoodTestStat::Evaluate - Use Offset in creating NLL " << endl ;
96 
97  created = kTRUE ;
98  delete allParams;
99  if (fPrintLevel > 1) cout << "creating NLL " << fNll << " with data = " << &data << endl ;
100  }
101  if (reuse && !created) {
102  if (fPrintLevel > 1) cout << "reusing NLL " << fNll << " new data = " << &data << endl ;
103  fNll->setData(data,kFALSE) ;
104  }
105  // print data in case of number counting (simple data sets)
106  if (fPrintLevel > 1 && data.numEntries() == 1) {
107  std::cout << "Data set used is: ";
108  RooStats::PrintListContent(*data.get(0), std::cout);
109  }
110 
111 
112  // make sure we set the variables attached to this nll
113  RooArgSet* attachedSet = fNll->getVariables();
114 
115  *attachedSet = paramsOfInterest;
116  RooArgSet* origAttachedSet = (RooArgSet*) attachedSet->snapshot();
117 
118  ///////////////////////////////////////////////////////////////////////
119  // New profiling based on RooMinimizer (allows for Minuit2)
120  // based on major speed increases seen by CMS for complex problems
121 
122 
123  // other order
124  // get the numerator
125  RooArgSet* snap = (RooArgSet*)paramsOfInterest.snapshot();
126 
127  tsw.Stop();
128  double createTime = tsw.CpuTime();
129  tsw.Start();
130 
131  // get the denominator
132  double uncondML = 0;
133  double fit_favored_mu = 0;
134  int statusD = 0;
135  RooArgSet * detOutput = 0;
136  if (type != 2) {
137  // minimize and count eval errors
138  fNll->clearEvalErrorLog();
139  if (fPrintLevel>1) std::cout << "Do unconditional fit" << std::endl;
140  RooFitResult* result = GetMinNLL();
141  if (result) {
142  uncondML = result->minNll();
143  statusD = result->status();
144 
145  // get best fit value for one-sided interval
146  if (firstPOI) fit_favored_mu = attachedSet->getRealValue(firstPOI->GetName()) ;
147 
148  // save this snapshot
149  if( fDetailedOutputEnabled ) {
150  detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitUncond_", fDetailedOutputWithErrorsAndPulls);
151  fDetailedOutput->addOwned(*detOutput);
152  delete detOutput;
153  }
154  delete result;
155  }
156  else {
157  return TMath::SignalingNaN(); // this should not really happen
158  }
159  }
160  tsw.Stop();
161  double fitTime1 = tsw.CpuTime();
162 
163  //double ret = 0;
164  int statusN = 0;
165  tsw.Start();
166 
167  double condML = 0;
168 
169  bool doConditionalFit = (type != 1);
170 
171  // skip the conditional ML (the numerator) only when fit value is smaller than test value
172  if (!fSigned && type==0 &&
173  ((fLimitType==oneSided && fit_favored_mu >= initial_mu_value) ||
174  (fLimitType==oneSidedDiscovery && fit_favored_mu <= initial_mu_value))) {
175  doConditionalFit = false;
176  condML = uncondML;
177  }
178 
179  if (doConditionalFit) {
180 
181  if (fPrintLevel>1) std::cout << "Do conditional fit " << std::endl;
182 
183 
184  // cout <<" reestablish snapshot"<<endl;
185  *attachedSet = *snap;
186 
187 
188  // set the POI to constant
189  RooLinkedListIter it = paramsOfInterest.iterator();
190  RooRealVar* tmpPar = NULL, *tmpParA=NULL;
191  while((tmpPar = (RooRealVar*)it.Next())){
192  tmpParA = dynamic_cast<RooRealVar*>( attachedSet->find(tmpPar->GetName()));
193  if (tmpParA) tmpParA->setConstant();
194  }
195 
196 
197  // check if there are non-const parameters so it is worth to do the minimization
198  RooArgSet allParams(*attachedSet);
200 
201  // in case no nuisance parameters are present
202  // no need to minimize just evaluate the nll
203  if (allParams.getSize() == 0 ) {
204  // be sure to evaluate with offsets
205  if (fLOffset) RooAbsReal::setHideOffset(false);
206  condML = fNll->getVal();
207  if (fLOffset) RooAbsReal::setHideOffset(true);
208  }
209  else {
210  fNll->clearEvalErrorLog();
211  RooFitResult* result = GetMinNLL();
212  if (result) {
213  condML = result->minNll();
214  statusN = result->status();
215  if( fDetailedOutputEnabled ) {
216  detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitCond_", fDetailedOutputWithErrorsAndPulls);
217  fDetailedOutput->addOwned(*detOutput);
218  delete detOutput;
219  }
220  delete result;
221  }
222  else {
223  return TMath::SignalingNaN(); // this should not really happen
224  }
225  }
226 
227  }
228 
229  tsw.Stop();
230  double fitTime2 = tsw.CpuTime();
231 
232  double pll = 0;
233  if (type != 0) {
234  // for conditional only or unconditional fits
235  // need to compute nll value without the offset
236  if (fLOffset) {
238  pll = fNll->getVal();
239  }
240  else {
241  if (type == 1)
242  pll = uncondML;
243  else if (type == 2)
244  pll = condML;
245  }
246  }
247  else { // type == 0
248  // for standard profile likelihood evaluations
249  pll = condML-uncondML;
250 
251  if (fSigned) {
252  if (pll<0.0) {
253  if (fPrintLevel > 0) std::cout << "pll is negative - setting it to zero " << std::endl;
254  pll = 0.0; // bad fit
255  }
256  if (fLimitType==oneSidedDiscovery ? (fit_favored_mu < initial_mu_value)
257  : (fit_favored_mu > initial_mu_value))
258  pll = -pll;
259  }
260  }
261 
262  if (fPrintLevel > 0) {
263  std::cout << "EvaluateProfileLikelihood - ";
264  if (type <= 1)
265  std::cout << "mu hat = " << fit_favored_mu << ", uncond ML = " << uncondML;
266  if (type != 1)
267  std::cout << ", cond ML = " << condML;
268  if (type == 0)
269  std::cout << " pll = " << pll;
270  std::cout << " time (create/fit1/2) " << createTime << " , " << fitTime1 << " , " << fitTime2
271  << std::endl;
272  }
273 
274 
275  // need to restore the values ?
276  *attachedSet = *origAttachedSet;
277 
278  delete attachedSet;
279  delete origAttachedSet;
280  delete snap;
281 
282  if (!reuse) {
283  delete fNll;
284  fNll = 0;
285  }
286 
288 
289  if(statusN!=0 || statusD!=0) {
290  return -1; // indicate failed fit (WVE is not used anywhere yet)
291  }
292 
293  return pll;
294 
295  }
296 
297 ////////////////////////////////////////////////////////////////////////////////
298 /// find minimum of NLL using RooMinimizer
299 
301 
302  const auto& config = GetGlobalRooStatsConfig();
303  RooMinimizer minim(*fNll);
304  minim.setStrategy(fStrategy);
305  minim.setEvalErrorWall(config.useEvalErrorWall);
306  //LM: RooMinimizer.setPrintLevel has +1 offset - so subtract here -1 + an extra -1
307  int level = (fPrintLevel == 0) ? -1 : fPrintLevel -2;
308  minim.setPrintLevel(level);
309  minim.setEps(fTolerance);
310  // this causes a memory leak
311  minim.optimizeConst(2);
312  TString minimizer = fMinimizer;
314  if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
315  int status;
316  for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
317  status = minim.minimize(minimizer,algorithm);
318  if (status%1000 == 0) { // ignore erros from Improve
319  break;
320  } else if (tries < maxtries) {
321  cout << " ----> Doing a re-scan first" << endl;
322  minim.minimize(minimizer,"Scan");
323  if (tries == 2) {
324  if (fStrategy == 0 ) {
325  cout << " ----> trying with strategy = 1" << endl;;
326  minim.setStrategy(1);
327  }
328  else
329  tries++; // skip this trial if strategy is already 1
330  }
331  if (tries == 3) {
332  cout << " ----> trying with improve" << endl;;
333  minimizer = "Minuit";
334  algorithm = "migradimproved";
335  }
336  }
337  }
338 
339  //how to get cov quality faster?
340  return minim.save();
341  //minim.optimizeConst(false);
342 }
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooMinimizer::setEps
void setEps(Double_t eps)
Change MINUIT epsilon.
Definition: RooMinimizer.cxx:225
RooMsgService.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooMinimizer::save
RooFitResult * save(const char *name=0, const char *title=0)
Save and return a RooFitResult snaphot of current minimizer status.
Definition: RooMinimizer.cxx:635
RooAbsData
Definition: RooAbsData.h:46
RooArgSet::getRealValue
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:474
RooArgSet.h
DetailedOutputAggregator.h
TStopwatch.h
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
RooMinimizer::setStrategy
void setStrategy(Int_t strat)
Change MINUIT strategy to istrat.
Definition: RooMinimizer.cxx:176
RooFit::Offset
RooCmdArg Offset(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:212
RooAbsCollection::find
RooAbsArg * find(const char *name) const
Find object with given name in list.
Definition: RooAbsCollection.cxx:813
RooStats::ProfileLikelihoodTestStat::fgAlwaysReuseNll
static Bool_t fgAlwaysReuseNll
Definition: ProfileLikelihoodTestStat.h:181
TString
Definition: TString.h:136
RooFit::MsgLevel
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
RooDataSet.h
RooNLLVar.h
RooLinkedListIter::Next
TObject * Next() override
Definition: RooLinkedListIter.h:220
bool
RooFitResult
Definition: RooFitResult.h:40
RooAbsCollection::iterator
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
Definition: RooAbsCollection.h:126
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
RooFit::GlobalObservables
RooCmdArg GlobalObservables(const RooArgSet &globs)
Definition: RooGlobalFunc.cxx:201
RooStats::ProfileLikelihoodTestStat::EvaluateProfileLikelihood
virtual Double_t EvaluateProfileLikelihood(int type, RooAbsData &data, RooArgSet &paramsOfInterest)
internal function to evaluate test statistics can do depending on type:
Definition: ProfileLikelihoodTestStat.cxx:56
RooStats::ProfileLikelihoodTestStat::GetMinNLL
RooFitResult * GetMinNLL()
find minimum of NLL using RooMinimizer
Definition: ProfileLikelihoodTestStat.cxx:300
RooLinkedListIter
A wrapper around TIterator derivatives.
Definition: RooLinkedListIter.h:204
RooFit::FATAL
@ FATAL
Definition: RooGlobalFunc.h:65
RooFitResult::status
Int_t status() const
Definition: RooFitResult.h:77
RooAbsReal::setHideOffset
static void setHideOffset(Bool_t flag)
Definition: RooAbsReal.cxx:117
RooPullVar.h
RooFit::Constrain
RooCmdArg Constrain(const RooArgSet &params)
Definition: RooGlobalFunc.cxx:200
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooAbsData::get
virtual const RooArgSet * get() const
Definition: RooAbsData.h:89
RooMsgService::setGlobalKillBelow
void setGlobalKillBelow(RooFit::MsgLevel level)
Definition: RooMsgService.h:160
RooStats::PrintListContent
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
Definition: RooStatsUtils.cxx:317
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
RooStats::ProfileLikelihoodTestStat::SetAlwaysReuseNLL
static void SetAlwaysReuseNLL(Bool_t flag)
Definition: ProfileLikelihoodTestStat.cxx:47
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
RooFit::ConditionalObservables
RooCmdArg ConditionalObservables(const RooArgSet &set)
Definition: RooGlobalFunc.cxx:196
TMath::SignalingNaN
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition: TMath.h:908
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
RooFitResult.h
RooStatsUtils.h
RooMinimizer::minimize
Int_t minimize(const char *type, const char *alg=0)
Minimise the function passed in the constructor.
Definition: RooMinimizer.cxx:312
Double_t
double Double_t
Definition: RtypesCore.h:59
ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo
static const std::string & DefaultMinimizerAlgo()
Definition: MinimizerOptions.cxx:92
RooMinimizer.h
TStopwatch
Definition: TStopwatch.h:28
RooStats::RemoveConstantParameters
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooAbsRealLValue::setConstant
void setConstant(Bool_t value=kTRUE)
Definition: RooAbsRealLValue.h:115
RooFit::CloneData
RooCmdArg CloneData(Bool_t flag)
Definition: RooGlobalFunc.cxx:209
RooAbsCollection::Print
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
Definition: RooAbsCollection.h:199
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooStats::GetGlobalRooStatsConfig
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
Definition: RooStatsUtils.cxx:33
type
int type
Definition: TGX11.cxx:121
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
RooRealVar
Definition: RooRealVar.h:36
RooMsgService::instance
static RooMsgService & instance()
Return reference to singleton instance.
Definition: RooMsgService.cxx:363
RooMinimizer
Definition: RooMinimizer.h:40
RooMinimizer::setEvalErrorWall
void setEvalErrorWall(Bool_t flag)
Definition: RooMinimizer.h:52
RooMinimizer::optimizeConst
void optimizeConst(Int_t flag)
If flag is true, perform constant term optimization on function being minimized.
Definition: RooMinimizer.cxx:603
RooMsgService::globalKillBelow
RooFit::MsgLevel globalKillBelow() const
Definition: RooMsgService.h:161
RooMinimizer::setPrintLevel
Int_t setPrintLevel(Int_t newLevel)
Change the MINUIT internal printing level.
Definition: RooMinimizer.cxx:591
RooFitResult::minNll
Double_t minNll() const
Definition: RooFitResult.h:98
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooProfileLL.h
RooArgSet
Definition: RooArgSet.h:28
ProfileLikelihoodTestStat.h