ROOT  6.06/09
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 
13 #include "RooFitResult.h"
14 #include "RooPullVar.h"
16 
17 #include "RooProfileLL.h"
18 #include "RooNLLVar.h"
19 #include "RooMsgService.h"
20 #include "RooMinimizer.h"
21 #include "RooArgSet.h"
22 #include "RooDataSet.h"
23 #include "TStopwatch.h"
24 
25 #include "RooStats/RooStatsUtils.h"
26 
27 using namespace std;
28 
30 
31 void RooStats::ProfileLikelihoodTestStat::SetAlwaysReuseNLL(Bool_t flag) { fgAlwaysReuseNll = flag ; }
32 
34  // interna function to evaluate test statistics
35  // can do depending on type:
36  // type = 0 standard evaluation, type = 1 find only unconditional NLL minimum, type = 2 conditional MLL
37 
38  if( fDetailedOutputEnabled && fDetailedOutput ) {
39  delete fDetailedOutput;
40  fDetailedOutput = 0;
41  }
42  if( fDetailedOutputEnabled && !fDetailedOutput ) {
43  fDetailedOutput = new RooArgSet();
44  }
45 
46  //data.Print("V");
47 
48  TStopwatch tsw;
49  tsw.Start();
50 
51  double initial_mu_value = 0;
52  RooRealVar* firstPOI = dynamic_cast<RooRealVar*>( paramsOfInterest.first());
53  if (firstPOI) initial_mu_value = firstPOI->getVal();
54  //paramsOfInterest.getRealValue(firstPOI->GetName());
55  if (fPrintLevel > 1) {
56  cout << "POIs: " << endl;
57  paramsOfInterest.Print("v");
58  }
59 
62 
63  // simple
64  Bool_t reuse=(fReuseNll || fgAlwaysReuseNll) ;
65 
66  Bool_t created(kFALSE) ;
67  if (!reuse || fNll==0) {
68  RooArgSet* allParams = fPdf->getParameters(data);
70 
71  // need to call constrain for RooSimultaneous until stripDisconnected problem fixed
72  fNll = fPdf->createNLL(data, RooFit::CloneData(kFALSE),RooFit::Constrain(*allParams),RooFit::ConditionalObservables(fConditionalObs), RooFit::Offset(fLOffset));
73 
74  if (fPrintLevel > 0 && fLOffset) cout << "ProfileLikelihoodTestStat::Evaluate - Use Offset in creating NLL " << endl ;
75 
76  created = kTRUE ;
77  delete allParams;
78  if (fPrintLevel > 1) cout << "creating NLL " << fNll << " with data = " << &data << endl ;
79  }
80  if (reuse && !created) {
81  if (fPrintLevel > 1) cout << "reusing NLL " << fNll << " new data = " << &data << endl ;
82  fNll->setData(data,kFALSE) ;
83  }
84  // print data in case of number counting (simple data sets)
85  if (fPrintLevel > 1 && data.numEntries() == 1) {
86  std::cout << "Data set used is: ";
87  RooStats::PrintListContent(*data.get(0), std::cout);
88  }
89 
90 
91  // make sure we set the variables attached to this nll
92  RooArgSet* attachedSet = fNll->getVariables();
93 
94  *attachedSet = paramsOfInterest;
95  RooArgSet* origAttachedSet = (RooArgSet*) attachedSet->snapshot();
96 
97  ///////////////////////////////////////////////////////////////////////
98  // New profiling based on RooMinimizer (allows for Minuit2)
99  // based on major speed increases seen by CMS for complex problems
100 
101 
102  // other order
103  // get the numerator
104  RooArgSet* snap = (RooArgSet*)paramsOfInterest.snapshot();
105 
106  tsw.Stop();
107  double createTime = tsw.CpuTime();
108  tsw.Start();
109 
110  // get the denominator
111  double uncondML = 0;
112  double fit_favored_mu = 0;
113  int statusD = 0;
114  RooArgSet * detOutput = 0;
115  if (type != 2) {
116  // minimize and count eval errors
117  fNll->clearEvalErrorLog();
118  if (fPrintLevel>1) std::cout << "Do unconditional fit" << std::endl;
119  RooFitResult* result = GetMinNLL();
120  if (result) {
121  uncondML = result->minNll();
122  statusD = result->status();
123 
124  // get best fit value for one-sided interval
125  if (firstPOI) fit_favored_mu = attachedSet->getRealValue(firstPOI->GetName()) ;
126 
127  // save this snapshot
128  if( fDetailedOutputEnabled ) {
129  detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitUncond_", fDetailedOutputWithErrorsAndPulls);
130  fDetailedOutput->addOwned(*detOutput);
131  delete detOutput;
132  }
133  delete result;
134  }
135  else {
136  return TMath::SignalingNaN(); // this should not really happen
137  }
138  }
139  tsw.Stop();
140  double fitTime1 = tsw.CpuTime();
141 
142  //double ret = 0;
143  int statusN = 0;
144  tsw.Start();
145 
146  double condML = 0;
147 
148  bool doConditionalFit = (type != 1);
149 
150  // skip the conditional ML (the numerator) only when fit value is smaller than test value
151  if (!fSigned && type==0 &&
152  ((fLimitType==oneSided && fit_favored_mu >= initial_mu_value) ||
153  (fLimitType==oneSidedDiscovery && fit_favored_mu <= initial_mu_value))) {
154  doConditionalFit = false;
155  condML = uncondML;
156  }
157 
158  if (doConditionalFit) {
159 
160  if (fPrintLevel>1) std::cout << "Do conditional fit " << std::endl;
161 
162 
163  // cout <<" reestablish snapshot"<<endl;
164  *attachedSet = *snap;
165 
166 
167  // set the POI to constant
168  RooLinkedListIter it = paramsOfInterest.iterator();
169  RooRealVar* tmpPar = NULL, *tmpParA=NULL;
170  while((tmpPar = (RooRealVar*)it.Next())){
171  tmpParA = dynamic_cast<RooRealVar*>( attachedSet->find(tmpPar->GetName()));
172  if (tmpParA) tmpParA->setConstant();
173  }
174 
175 
176  // check if there are non-const parameters so it is worth to do the minimization
177  RooArgSet allParams(*attachedSet);
179 
180  // in case no nuisance parameters are present
181  // no need to minimize just evaluate the nll
182  if (allParams.getSize() == 0 ) {
183  // be sure to evaluate with offsets
184  if (fLOffset) RooAbsReal::setHideOffset(false);
185  condML = fNll->getVal();
186  if (fLOffset) RooAbsReal::setHideOffset(true);
187  }
188  else {
189  fNll->clearEvalErrorLog();
190  RooFitResult* result = GetMinNLL();
191  if (result) {
192  condML = result->minNll();
193  statusN = result->status();
194  if( fDetailedOutputEnabled ) {
195  detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitCond_", fDetailedOutputWithErrorsAndPulls);
196  fDetailedOutput->addOwned(*detOutput);
197  delete detOutput;
198  }
199  delete result;
200  }
201  else {
202  return TMath::SignalingNaN(); // this should not really happen
203  }
204  }
205 
206  }
207 
208  tsw.Stop();
209  double fitTime2 = tsw.CpuTime();
210 
211  double pll = 0;
212  if (type != 0) {
213  // for conditional only or unconditional fits
214  // need to compute nll value without the offset
215  if (fLOffset) {
217  pll = fNll->getVal();
218  }
219  else {
220  if (type == 1)
221  pll = uncondML;
222  else if (type == 2)
223  pll = condML;
224  }
225  }
226  else { // type == 0
227  // for standard profile likelihood evaluations
228  pll = condML-uncondML;
229 
230  if (fSigned) {
231  if (pll<0.0) {
232  if (fPrintLevel > 0) std::cout << "pll is negative - setting it to zero " << std::endl;
233  pll = 0.0; // bad fit
234  }
235  if (fLimitType==oneSidedDiscovery ? (fit_favored_mu < initial_mu_value)
236  : (fit_favored_mu > initial_mu_value))
237  pll = -pll;
238  }
239  }
240 
241  if (fPrintLevel > 0) {
242  std::cout << "EvaluateProfileLikelihood - ";
243  if (type <= 1)
244  std::cout << "mu hat = " << fit_favored_mu << ", uncond ML = " << uncondML;
245  if (type != 1)
246  std::cout << ", cond ML = " << condML;
247  if (type == 0)
248  std::cout << " pll = " << pll;
249  std::cout << " time (create/fit1/2) " << createTime << " , " << fitTime1 << " , " << fitTime2
250  << std::endl;
251  }
252 
253 
254  // need to restore the values ?
255  *attachedSet = *origAttachedSet;
256 
257  delete attachedSet;
258  delete origAttachedSet;
259  delete snap;
260 
261  if (!reuse) {
262  delete fNll;
263  fNll = 0;
264  }
265 
267 
268  if(statusN!=0 || statusD!=0) {
269  return -1; // indicate failed fit (WVE is not used anywhere yet)
270  }
271 
272  return pll;
273 
274  }
275 
277  //find minimum of NLL using RooMinimizer
278 
279  RooMinimizer minim(*fNll);
280  minim.setStrategy(fStrategy);
281  //LM: RooMinimizer.setPrintLevel has +1 offset - so subtruct here -1 + an extra -1
282  int level = (fPrintLevel == 0) ? -1 : fPrintLevel -2;
283  minim.setPrintLevel(level);
284  minim.setEps(fTolerance);
285  // this cayses a memory leak
286  minim.optimizeConst(2);
287  TString minimizer = fMinimizer;
289  if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
290  int status;
291  for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
292  status = minim.minimize(minimizer,algorithm);
293  if (status%1000 == 0) { // ignore erros from Improve
294  break;
295  } else if (tries < maxtries) {
296  cout << " ----> Doing a re-scan first" << endl;
297  minim.minimize(minimizer,"Scan");
298  if (tries == 2) {
299  if (fStrategy == 0 ) {
300  cout << " ----> trying with strategy = 1" << endl;;
301  minim.setStrategy(1);
302  }
303  else
304  tries++; // skip this trial if stratehy is already 1
305  }
306  if (tries == 3) {
307  cout << " ----> trying with improve" << endl;;
308  minimizer = "Minuit";
309  algorithm = "migradimproved";
310  }
311  }
312  }
313 
314  //how to get cov quality faster?
315  return minim.save();
316  //minim.optimizeConst(false);
317 }
RooCmdArg Offset(Bool_t flag=kTRUE)
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
void optimizeConst(Int_t flag)
If flag is true, perform constant term optimization on function being minimized.
RooCmdArg CloneData(Bool_t flag)
Basic string class.
Definition: TString.h:137
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual Double_t EvaluateProfileLikelihood(int type, RooAbsData &data, RooArgSet &paramsOfInterest)
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
static void setHideOffset(Bool_t flag)
Definition: RooAbsReal.cxx:118
void setEps(Double_t eps)
Change MINUIT epsilon.
void setStrategy(Int_t strat)
Change MINUIT strategy to istrat.
RooAbsArg * first() const
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
virtual TObject * Next()
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
virtual const RooArgSet * get() const
Definition: RooAbsData.h:77
RooFit::MsgLevel globalKillBelow() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
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.
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:291
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
void setGlobalKillBelow(RooFit::MsgLevel level)
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
Double_t SignalingNaN()
Definition: TMath.h:642
static const std::string & DefaultMinimizerAlgo()
Int_t status() const
Definition: RooFitResult.h:76
double Double_t
Definition: RtypesCore.h:55
Double_t minNll() const
Definition: RooFitResult.h:97
int type
Definition: TGX11.cxx:120
Int_t minimize(const char *type, const char *alg=0)
Int_t setPrintLevel(Int_t newLevel)
Change the MINUIT internal printing level.
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:73
#define NULL
Definition: Rtypes.h:82
double result[121]
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
Int_t getSize() const
RooCmdArg ConditionalObservables(const RooArgSet &set)
const Bool_t kTRUE
Definition: Rtypes.h:91
RooCmdArg Constrain(const RooArgSet &params)
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:527
Stopwatch class.
Definition: TStopwatch.h:30