Logo ROOT   6.10/09
Reference Guide
RooProfileLL.cxx
Go to the documentation of this file.
1  /*****************************************************************************
2  * Project: RooFit *
3  * *
4  * Copyright (c) 2000-2005, Regents of the University of California *
5  * and Stanford University. All rights reserved. *
6  * *
7  * Redistribution and use in source and binary forms, *
8  * with or without modification, are permitted according to the terms *
9  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
10  *****************************************************************************/
11 
12 /**
13 \file RooProfileLL.cxx
14 \class RooProfileLL
15 \ingroup Roofitcore
16 
17 Class RooProfileLL implements the profile likelihood estimator for
18 a given likelihood and set of parameters of interest. The value return by
19 RooProfileLL is the input likelihood nll minimized w.r.t all nuisance parameters
20 (which are all parameters except for those listed in the constructor) minus
21 the -log(L) of the best fit. Note that this function is slow to evaluate
22 as a MIGRAD minimization step is executed for each function evaluation
23 **/
24 
25 #include "Riostream.h"
26 
27 #include "RooFit.h"
28 #include "RooProfileLL.h"
29 #include "RooAbsReal.h"
30 #include "RooMinuit.h"
31 #include "RooMinimizer.h"
32 #include "RooMsgService.h"
33 #include "RooRealVar.h"
34 #include "RooMsgService.h"
35 
36 using namespace std ;
37 
39 
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 /// Default constructor
43 /// Should only be used by proof.
44 
46  RooAbsReal("RooProfileLL","RooProfileLL"),
47  _nll(),
48  _obs("paramOfInterest","Parameters of interest",this),
49  _par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE),
50  _startFromMin(kTRUE),
51  _minimizer(0),
52  _absMinValid(kFALSE),
53  _absMin(0),
54  _neval(0)
55 {
56  _piter = _par.createIterator() ;
57  _oiter = _obs.createIterator() ;
58 }
59 
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Constructor of profile likelihood given input likelihood nll w.r.t
63 /// the given set of variables. The input log likelihood is minimized w.r.t
64 /// to all other variables of the likelihood at each evaluation and the
65 /// value of the global log likelihood minimum is always subtracted.
66 
67 RooProfileLL::RooProfileLL(const char *name, const char *title,
68  RooAbsReal& nllIn, const RooArgSet& observables) :
69  RooAbsReal(name,title),
70  _nll("input","-log(L) function",this,nllIn),
71  _obs("paramOfInterest","Parameters of interest",this),
72  _par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE),
73  _startFromMin(kTRUE),
74  _minimizer(0),
75  _absMinValid(kFALSE),
76  _absMin(0),
77  _neval(0)
78 {
79  // Determine actual parameters and observables
80  RooArgSet* actualObs = nllIn.getObservables(observables) ;
81  RooArgSet* actualPars = nllIn.getParameters(observables) ;
82 
83  _obs.add(*actualObs) ;
84  _par.add(*actualPars) ;
85 
86  delete actualObs ;
87  delete actualPars ;
88 
91 }
92 
93 
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Copy constructor
97 
98 RooProfileLL::RooProfileLL(const RooProfileLL& other, const char* name) :
99  RooAbsReal(other,name),
100  _nll("nll",this,other._nll),
101  _obs("obs",this,other._obs),
102  _par("par",this,other._par),
104  _minimizer(0),
106  _absMin(0),
107  _paramFixed(other._paramFixed),
108  _neval(0)
109 {
112 
115 
116 }
117 
118 
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 /// Destructor
122 
124 {
125  // Delete instance of minuit if it was ever instantiated
126  if (_minimizer) {
127  delete _minimizer ;
128  }
129 
130  delete _piter ;
131  delete _oiter ;
132 }
133 
134 
135 
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 
140 {
141  validateAbsMin() ;
142  return _paramAbsMin ;
143 }
144 
145 
146 ////////////////////////////////////////////////////////////////////////////////
147 
149 {
150  validateAbsMin() ;
151  return _obsAbsMin ;
152 }
153 
154 
155 
156 
157 ////////////////////////////////////////////////////////////////////////////////
158 /// Optimized implementation of createProfile for profile likelihoods.
159 /// Return profile of original function in terms of stated parameters
160 /// of interest rather than profiling recursively.
161 
163 {
164  return nll().createProfile(paramsOfInterest) ;
165 }
166 
167 
168 
169 
170 ////////////////////////////////////////////////////////////////////////////////
171 
173 {
174  coutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") Creating instance of MINUIT" << endl ;
175 
178  _minimizer = new MINIMIZER(const_cast<RooAbsReal&>(_nll.arg())) ;
180 
181  //_minimizer->setPrintLevel(-999) ;
182  //_minimizer->setNoWarn() ;
183  //_minimizer->setVerbose(1) ;
184 }
185 
186 
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// Evaluate profile likelihood by minimizing likelihood w.r.t. all
190 /// parameters that are not considered observables of this profile
191 /// likelihood object.
192 
194 {
195  // Instantiate minimizer if we haven't done that already
196  if (!_minimizer) {
198  }
199 
200  // Save current value of observables
201  RooArgSet* obsSetOrig = (RooArgSet*) _obs.snapshot() ;
202 
203  validateAbsMin() ;
204 
205 
206  // Set all observables constant in the minimization
207  const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",kTRUE) ;
208  ccoutP(Eval) << "." ; ccoutP(Eval).flush() ;
209 
210  // If requested set initial parameters to those corresponding to absolute minimum
211  if (_startFromMin) {
212  const_cast<RooProfileLL&>(*this)._par = _paramAbsMin ;
213  }
214 
215  _minimizer->zeroEvalCount() ;
216 
217  //TString minim=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str();
218  //TString algorithm = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str();
219  //if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
220  //_minimizer->minimize(minim.Data(),algorithm.Data());
221 
222  _minimizer->migrad() ;
223  _neval = _minimizer->evalCounter() ;
224 
225  // Restore original values and constant status of observables
226  TIterator* iter = obsSetOrig->createIterator() ;
227  RooRealVar* var ;
228  while((var=dynamic_cast<RooRealVar*>(iter->Next()) ) ) {
229  RooRealVar* target = (RooRealVar*) _obs.find(var->GetName()) ;
230  target->setVal(var->getVal()) ;
231  target->setConstant(var->isConstant()) ;
232  }
233  delete iter ;
234  delete obsSetOrig ;
235 
236  return _nll - _absMin ;
237 }
238 
239 
240 
241 ////////////////////////////////////////////////////////////////////////////////
242 /// Check that parameters and likelihood value for 'best fit' are still valid. If not,
243 /// because the best fit has never been calculated, or because constant parameters have
244 /// changed value or parameters have changed const/float status, the minimum is recalculated
245 
247 {
248  // Check if constant status of any of the parameters have changed
249  if (_absMinValid) {
250  _piter->Reset() ;
251  RooAbsArg* par ;
252  while((par=(RooAbsArg*)_piter->Next())) {
253  if (_paramFixed[par->GetName()] != par->isConstant()) {
254  cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") constant status of parameter " << par->GetName() << " has changed from "
255  << (_paramFixed[par->GetName()]?"fixed":"floating") << " to " << (par->isConstant()?"fixed":"floating")
256  << ", recalculating absolute minimum" << endl ;
257  _absMinValid = kFALSE ;
258  break ;
259  }
260  }
261  }
262 
263 
264  // If we don't have the absolute minimum w.r.t all observables, calculate that first
265  if (!_absMinValid) {
266 
267  cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") determining minimum likelihood for current configurations w.r.t all observable" << endl ;
268 
269 
270  if (!_minimizer) {
272  }
273 
274  // Save current values of non-marginalized parameters
275  RooArgSet* obsStart = (RooArgSet*) _obs.snapshot(kFALSE) ;
276 
277  // Start from previous global minimum
278  if (_paramAbsMin.getSize()>0) {
279  const_cast<RooSetProxy&>(_par).assignValueOnly(_paramAbsMin) ;
280  }
281  if (_obsAbsMin.getSize()>0) {
282  const_cast<RooSetProxy&>(_obs).assignValueOnly(_obsAbsMin) ;
283  }
284 
285  // Find minimum with all observables floating
286  const_cast<RooSetProxy&>(_obs).setAttribAll("Constant",kFALSE) ;
287 
288  //TString minim=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str();
289  //TString algorithm = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str();
290  //if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
291  //_minimizer->minimize(minim.Data(),algorithm.Data());
292  _minimizer->migrad() ;
293 
294  // Save value and remember
295  _absMin = _nll ;
296  _absMinValid = kTRUE ;
297 
298  // Save parameter values at abs minimum as well
300 
301  // Only store non-constant parameters here!
302  RooArgSet* tmp = (RooArgSet*) _par.selectByAttrib("Constant",kFALSE) ;
303  _paramAbsMin.addClone(*tmp) ;
304  delete tmp ;
305 
307 
308  // Save constant status of all parameters
309  _piter->Reset() ;
310  RooAbsArg* par ;
311  while((par=(RooAbsArg*)_piter->Next())) {
312  _paramFixed[par->GetName()] = par->isConstant() ;
313  }
314 
315  if (dologI(Minimization)) {
316  cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") minimum found at (" ;
317 
318  RooAbsReal* arg ;
319  Bool_t first=kTRUE ;
320  _oiter->Reset() ;
321  while ((arg=(RooAbsReal*)_oiter->Next())) {
322  ccxcoutI(Minimization) << (first?"":", ") << arg->GetName() << "=" << arg->getVal() ;
323  first=kFALSE ;
324  }
325  ccxcoutI(Minimization) << ")" << endl ;
326  }
327 
328  // Restore original parameter values
329  const_cast<RooSetProxy&>(_obs) = *obsStart ;
330  delete obsStart ;
331 
332  }
333 }
334 
335 
336 
337 ////////////////////////////////////////////////////////////////////////////////
338 
339 Bool_t RooProfileLL::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
340  Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
341 {
342  if (_minimizer) {
343  delete _minimizer ;
344  _minimizer = 0 ;
345  }
346  return kFALSE ;
347 }
348 
349 
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
double par[1]
Definition: unuranDistr.cxx:38
TIterator * createIterator(Bool_t dir=kIterForward) const
#define cxcoutI(a)
Definition: RooMsgService.h:83
virtual void Reset()=0
RooArgSet _paramAbsMin
Definition: RooProfileLL.h:69
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:26
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Definition: RooAbsArg.h:194
#define coutI(a)
Definition: RooMsgService.h:31
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
RooProfileLL()
Default constructor Should only be used by proof.
bool Bool_t
Definition: RtypesCore.h:59
const RooArgSet & bestFitObs() const
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
Iterator abstract base class.
Definition: TIterator.h:30
const RooArgSet & bestFitParams() const
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Optimized implementation of createProfile for profile likelihoods.
Bool_t _absMinValid
Internal minuit instance.
Definition: RooProfileLL.h:67
Bool_t _startFromMin
Definition: RooProfileLL.h:60
Bool_t silentMode() const
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
void validateAbsMin() const
Check that parameters and likelihood value for &#39;best fit&#39; are still valid.
RooAbsCollection * selectByAttrib(const char *name, Bool_t value) const
Create a subset of the current collection, consisting only of those elements with the specified attri...
#define dologI(a)
Definition: RooMsgService.h:64
std::map< std::string, bool > _paramFixed
Definition: RooProfileLL.h:71
RooRealProxy _nll
Definition: RooProfileLL.h:57
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
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
RooCmdArg Parameters(const RooArgSet &params)
void setConstant(Bool_t value=kTRUE)
RooArgSet _obsAbsMin
Definition: RooProfileLL.h:70
RooAbsReal & nll()
Definition: RooProfileLL.h:39
RooSetProxy _obs
Definition: RooProfileLL.h:58
#define ccoutP(a)
Definition: RooMsgService.h:39
TIterator * _piter
Definition: RooProfileLL.h:62
void setSilentMode(Bool_t flag)
virtual ~RooProfileLL()
Destructor.
const Bool_t kFALSE
Definition: RtypesCore.h:92
virtual Bool_t redirectServersHook(const RooAbsCollection &, Bool_t, Bool_t, Bool_t)
RooSetProxy _par
Definition: RooProfileLL.h:59
#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
Double_t evaluate() const
Evaluate profile likelihood by minimizing likelihood w.r.t.
#define MINIMIZER
Definition: RooProfileLL.h:24
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects...
MINIMIZER * _minimizer
Iterator of profile likelihood output parameter(s)
Definition: RooProfileLL.h:65
#define ccxcoutI(a)
Definition: RooMsgService.h:84
void initializeMinimizer() const
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:461
virtual TObject * Next()=0
RooSetProxy is the concrete proxy for RooArgSet objects.
Definition: RooSetProxy.h:24
Double_t _absMin
Definition: RooProfileLL.h:68
Definition: first.py:1
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
const Bool_t kTRUE
Definition: RtypesCore.h:91
Bool_t isConstant() const
Definition: RooAbsArg.h:266
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Overloaded RooArgSet::add() method inserts &#39;var&#39; into set and registers &#39;var&#39; as server to owner with...
TIterator * _oiter
Iterator over profile likelihood parameters to be minimized.
Definition: RooProfileLL.h:63