Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
15ProfileLikelihoodTestStat is an implementation of the TestStatistic interface
16that calculates the profile likelihood ratio at a particular parameter point
17given a dataset. It does not constitute a statistical test, for that one may
18either 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
42
43using namespace std;
44
46
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->assign(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->assign(*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->assign(*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}
const Bool_t kFALSE
Definition RtypesCore.h:101
const Bool_t kTRUE
Definition RtypesCore.h:100
int type
Definition TGX11.cxx:121
static const std::string & DefaultMinimizerAlgo()
Int_t getSize() const
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.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
RooAbsArg * first() const
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
RooLinkedListIter iterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:82
virtual const RooArgSet * get() const
Definition RooAbsData.h:128
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
void setConstant(Bool_t value=kTRUE)
static void setHideOffset(Bool_t flag)
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:94
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:35
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition RooArgSet.h:158
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Double_t minNll() const
Return minimized -log(L) value.
Int_t status() const
Return MINUIT status code.
A wrapper around TIterator derivatives.
TObject * Next() override
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
void setEvalErrorWall(Bool_t flag)
RooFitResult * save(const char *name=0, const char *title=0)
Save and return a RooFitResult snapshot of current minimizer status.
Int_t minimize(const char *type, const char *alg=0)
Minimise the function passed in the constructor.
void setEps(Double_t eps)
Change MINUIT epsilon.
Int_t setPrintLevel(Int_t newLevel)
Change the MINUIT internal printing level.
void optimizeConst(Int_t flag)
If flag is true, perform constant term optimization on function being minimized.
void setStrategy(Int_t strat)
Change MINUIT strategy to istrat.
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
static RooArgSet * GetAsArgSet(RooFitResult *result, TString prefix="", bool withErrorsAndPulls=false)
static function to translate the given fit result to a RooArgSet in a generic way.
virtual Double_t EvaluateProfileLikelihood(int type, RooAbsData &data, RooArgSet &paramsOfInterest)
internal function to evaluate test statistics can do depending on type:
RooFitResult * GetMinNLL()
find minimum of NLL using RooMinimizer
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
void Stop()
Stop the stopwatch.
Basic string class.
Definition TString.h:136
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg CloneData(Bool_t flag)
RooCmdArg Offset(Bool_t flag=kTRUE)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
void RemoveConstantParameters(RooArgSet *set)
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
void PrintListContent(const RooArgList &l, std::ostream &os=std::cout)
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition TMath.h:858