/***************************************************************************** * Project: RooFit * * * * Copyright (c) 2000-2005, Regents of the University of California * * and Stanford University. All rights reserved. * * * * Redistribution and use in source and binary forms, * * with or without modification, are permitted according to the terms * * listed in LICENSE (http://roofit.sourceforge.net/license.txt) * *****************************************************************************/ ////////////////////////////////////////////////////////////////////////////// // // BEGIN_HTML // Class RooProfileLL implements the profile likelihood estimator for // a given likelihood and set of parameters of interest. The value return by // RooProfileLL is the input likelihood nll minimized w.r.t all nuisance parameters // (which are all parameters except for those listed in the constructor) minus // the -log(L) of the best fit. Note that this function is slow to evaluate // as a MIGRAD minimization step is executed for each function evaluation // END_HTML // #include "Riostream.h" #include "RooFit.h" #include "RooProfileLL.h" #include "RooAbsReal.h" #include "RooMinuit.h" #include "RooMinimizer.h" #include "RooMsgService.h" #include "RooRealVar.h" #include "RooMsgService.h" using namespace std ; ClassImp(RooProfileLL) //_____________________________________________________________________________ RooProfileLL::RooProfileLL() : RooAbsReal("RooProfileLL","RooProfileLL"), _nll(), _obs("paramOfInterest","Parameters of interest",this), _par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE), _startFromMin(kTRUE), _minimizer(0), _absMinValid(kFALSE), _absMin(0), _neval(0) { // Default constructor // Should only be used by proof. _piter = _par.createIterator() ; _oiter = _obs.createIterator() ; } //_____________________________________________________________________________ RooProfileLL::RooProfileLL(const char *name, const char *title, RooAbsReal& nllIn, const RooArgSet& observables) : RooAbsReal(name,title), _nll("input","-log(L) function",this,nllIn), _obs("paramOfInterest","Parameters of interest",this), _par("nuisanceParam","Nuisance parameters",this,kFALSE,kFALSE), _startFromMin(kTRUE), _minimizer(0), _absMinValid(kFALSE), _absMin(0), _neval(0) { // Constructor of profile likelihood given input likelihood nll w.r.t // the given set of variables. The input log likelihood is minimized w.r.t // to all other variables of the likelihood at each evaluation and the // value of the global log likelihood minimum is always subtracted. // Determine actual parameters and observables RooArgSet* actualObs = nllIn.getObservables(observables) ; RooArgSet* actualPars = nllIn.getParameters(observables) ; _obs.add(*actualObs) ; _par.add(*actualPars) ; delete actualObs ; delete actualPars ; _piter = _par.createIterator() ; _oiter = _obs.createIterator() ; } //_____________________________________________________________________________ RooProfileLL::RooProfileLL(const RooProfileLL& other, const char* name) : RooAbsReal(other,name), _nll("nll",this,other._nll), _obs("obs",this,other._obs), _par("par",this,other._par), _startFromMin(other._startFromMin), _minimizer(0), _absMinValid(kFALSE), _absMin(0), _paramFixed(other._paramFixed), _neval(0) { // Copy constructor _piter = _par.createIterator() ; _oiter = _obs.createIterator() ; _paramAbsMin.addClone(other._paramAbsMin) ; _obsAbsMin.addClone(other._obsAbsMin) ; } //_____________________________________________________________________________ RooProfileLL::~RooProfileLL() { // Destructor // Delete instance of minuit if it was ever instantiated if (_minimizer) { delete _minimizer ; } delete _piter ; delete _oiter ; } //_____________________________________________________________________________ const RooArgSet& RooProfileLL::bestFitParams() const { validateAbsMin() ; return _paramAbsMin ; } //_____________________________________________________________________________ const RooArgSet& RooProfileLL::bestFitObs() const { validateAbsMin() ; return _obsAbsMin ; } //_____________________________________________________________________________ RooAbsReal* RooProfileLL::createProfile(const RooArgSet& paramsOfInterest) { // Optimized implementation of createProfile for profile likelihoods. // Return profile of original function in terms of stated parameters // of interest rather than profiling recursively. return nll().createProfile(paramsOfInterest) ; } //_____________________________________________________________________________ void RooProfileLL::initializeMinimizer() const { coutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") Creating instance of MINUIT" << endl ; Bool_t smode = RooMsgService::instance().silentMode() ; RooMsgService::instance().setSilentMode(kTRUE) ; _minimizer = new MINIMIZER(const_cast(_nll.arg())) ; if (!smode) RooMsgService::instance().setSilentMode(kFALSE) ; //_minimizer->setPrintLevel(-999) ; //_minimizer->setNoWarn() ; //_minimizer->setVerbose(1) ; } //_____________________________________________________________________________ Double_t RooProfileLL::evaluate() const { // Evaluate profile likelihood by minimizing likelihood w.r.t. all // parameters that are not considered observables of this profile // likelihood object. // Instantiate minimizer if we haven't done that already if (!_minimizer) { initializeMinimizer() ; } // Save current value of observables RooArgSet* obsSetOrig = (RooArgSet*) _obs.snapshot() ; validateAbsMin() ; // Set all observables constant in the minimization const_cast(_obs).setAttribAll("Constant",kTRUE) ; ccoutP(Eval) << "." ; ccoutP(Eval).flush() ; // If requested set initial parameters to those corresponding to absolute minimum if (_startFromMin) { const_cast(*this)._par = _paramAbsMin ; } _minimizer->zeroEvalCount() ; //TString minim=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(); //TString algorithm = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str(); //if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad //_minimizer->minimize(minim.Data(),algorithm.Data()); _minimizer->migrad() ; _neval = _minimizer->evalCounter() ; // Restore original values and constant status of observables TIterator* iter = obsSetOrig->createIterator() ; RooRealVar* var ; while((var=(RooRealVar*)iter->Next())) { RooRealVar* target = (RooRealVar*) _obs.find(var->GetName()) ; target->setVal(var->getVal()) ; target->setConstant(var->isConstant()) ; } delete iter ; delete obsSetOrig ; return _nll - _absMin ; } //_____________________________________________________________________________ void RooProfileLL::validateAbsMin() const { // Check that parameters and likelihood value for 'best fit' are still valid. If not, // because the best fit has never been calculated, or because constant parameters have // changed value or parameters have changed const/float status, the minimum is recalculated // Check if constant status of any of the parameters have changed if (_absMinValid) { _piter->Reset() ; RooAbsArg* par ; while((par=(RooAbsArg*)_piter->Next())) { if (_paramFixed[par->GetName()] != par->isConstant()) { cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") constant status of parameter " << par->GetName() << " has changed from " << (_paramFixed[par->GetName()]?"fixed":"floating") << " to " << (par->isConstant()?"fixed":"floating") << ", recalculating absolute minimum" << endl ; _absMinValid = kFALSE ; break ; } } } // If we don't have the absolute minimum w.r.t all observables, calculate that first if (!_absMinValid) { cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") determining minimum likelihood for current configurations w.r.t all observable" << endl ; if (!_minimizer) { initializeMinimizer() ; } // Save current values of non-marginalized parameters RooArgSet* obsStart = (RooArgSet*) _obs.snapshot(kFALSE) ; // Start from previous global minimum if (_paramAbsMin.getSize()>0) { const_cast(_par).assignValueOnly(_paramAbsMin) ; } if (_obsAbsMin.getSize()>0) { const_cast(_obs).assignValueOnly(_obsAbsMin) ; } // Find minimum with all observables floating const_cast(_obs).setAttribAll("Constant",kFALSE) ; //TString minim=::ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(); //TString algorithm = ROOT::Math::MinimizerOptions::DefaultMinimizerAlgo().c_str(); //if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad //_minimizer->minimize(minim.Data(),algorithm.Data()); _minimizer->migrad() ; // Save value and remember _absMin = _nll ; _absMinValid = kTRUE ; // Save parameter values at abs minimum as well _paramAbsMin.removeAll() ; // Only store non-constant parameters here! RooArgSet* tmp = (RooArgSet*) _par.selectByAttrib("Constant",kFALSE) ; _paramAbsMin.addClone(*tmp) ; delete tmp ; _obsAbsMin.addClone(_obs) ; // Save constant status of all parameters _piter->Reset() ; RooAbsArg* par ; while((par=(RooAbsArg*)_piter->Next())) { _paramFixed[par->GetName()] = par->isConstant() ; } if (dologI(Minimization)) { cxcoutI(Minimization) << "RooProfileLL::evaluate(" << GetName() << ") minimum found at (" ; RooAbsReal* arg ; Bool_t first=kTRUE ; _oiter->Reset() ; while ((arg=(RooAbsReal*)_oiter->Next())) { ccxcoutI(Minimization) << (first?"":", ") << arg->GetName() << "=" << arg->getVal() ; first=kFALSE ; } ccxcoutI(Minimization) << ")" << endl ; } // Restore original parameter values const_cast(_obs) = *obsStart ; delete obsStart ; } } //_____________________________________________________________________________ Bool_t RooProfileLL::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/, Bool_t /*nameChange*/, Bool_t /*isRecursive*/) { if (_minimizer) { delete _minimizer ; _minimizer = 0 ; } return kFALSE ; }