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
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 "RooMsgService.h"
35#include "RooMinimizer.h"
36#include "RooArgSet.h"
37#include "RooDataSet.h"
38#include "TStopwatch.h"
39
41
42using namespace std;
43
45
47
48////////////////////////////////////////////////////////////////////////////////
49/// internal function to evaluate test statistics
50/// can do depending on type:
51/// - type = 0 standard evaluation,
52/// - type = 1 find only unconditional NLL minimum,
53/// - type = 2 conditional MLL
54
56
57 if( fDetailedOutputEnabled && fDetailedOutput ) {
58 delete fDetailedOutput;
59 fDetailedOutput = 0;
60 }
61 if( fDetailedOutputEnabled && !fDetailedOutput ) {
62 fDetailedOutput = new RooArgSet();
63 }
64
65 //data.Print("V");
66
67 TStopwatch tsw;
68 tsw.Start();
69
70 double initial_mu_value = 0;
71 RooRealVar* firstPOI = dynamic_cast<RooRealVar*>( paramsOfInterest.first());
72 if (firstPOI) initial_mu_value = firstPOI->getVal();
73 //paramsOfInterest.getRealValue(firstPOI->GetName());
74 if (fPrintLevel > 1) {
75 cout << "POIs: " << endl;
76 paramsOfInterest.Print("v");
77 }
78
81
82 // simple
83 bool reuse=(fReuseNll || fgAlwaysReuseNll) ;
84
85 bool created(false) ;
86 if (!reuse || fNll==0) {
87 RooArgSet* allParams = fPdf->getParameters(data);
89
90 // need to call constrain for RooSimultaneous until stripDisconnected problem fixed
91 fNll = fPdf->createNLL(data, RooFit::CloneData(false),RooFit::Constrain(*allParams),
92 RooFit::GlobalObservables(fGlobalObs), RooFit::ConditionalObservables(fConditionalObs), RooFit::Offset(fLOffset));
93
94 if (fPrintLevel > 0 && fLOffset) cout << "ProfileLikelihoodTestStat::Evaluate - Use Offset in creating NLL " << endl ;
95
96 created = true ;
97 delete allParams;
98 if (fPrintLevel > 1) cout << "creating NLL " << fNll << " with data = " << &data << endl ;
99 }
100 if (reuse && !created) {
101 if (fPrintLevel > 1) cout << "reusing NLL " << fNll << " new data = " << &data << endl ;
102 fNll->setData(data,false) ;
103 }
104 // print data in case of number counting (simple data sets)
105 if (fPrintLevel > 1 && data.numEntries() == 1) {
106 std::cout << "Data set used is: ";
107 RooStats::PrintListContent(*data.get(0), std::cout);
108 }
109
110
111 // make sure we set the variables attached to this nll
112 RooArgSet* attachedSet = fNll->getVariables();
113
114 attachedSet->assign(paramsOfInterest);
115 RooArgSet* origAttachedSet = (RooArgSet*) attachedSet->snapshot();
116
117 ///////////////////////////////////////////////////////////////////////
118 // New profiling based on RooMinimizer (allows for Minuit2)
119 // based on major speed increases seen by CMS for complex problems
120
121
122 // other order
123 // get the numerator
124 RooArgSet* snap = (RooArgSet*)paramsOfInterest.snapshot();
125
126 tsw.Stop();
127 double createTime = tsw.CpuTime();
128 tsw.Start();
129
130 // get the denominator
131 double uncondML = 0;
132 double fit_favored_mu = 0;
133 int statusD = 0;
134 RooArgSet * detOutput = 0;
135 if (type != 2) {
136 // minimize and count eval errors
137 fNll->clearEvalErrorLog();
138 if (fPrintLevel>1) std::cout << "Do unconditional fit" << std::endl;
139 RooFitResult* result = GetMinNLL();
140 if (result) {
141 uncondML = result->minNll();
142 statusD = result->status();
143
144 // get best fit value for one-sided interval
145 if (firstPOI) fit_favored_mu = attachedSet->getRealValue(firstPOI->GetName()) ;
146
147 // save this snapshot
148 if( fDetailedOutputEnabled ) {
149 detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitUncond_", fDetailedOutputWithErrorsAndPulls);
150 fDetailedOutput->addOwned(*detOutput);
151 delete detOutput;
152 }
153 delete result;
154 }
155 else {
156 return TMath::SignalingNaN(); // this should not really happen
157 }
158 }
159 tsw.Stop();
160 double fitTime1 = tsw.CpuTime();
161
162 //double ret = 0;
163 int statusN = 0;
164 tsw.Start();
165
166 double condML = 0;
167
168 bool doConditionalFit = (type != 1);
169
170 // skip the conditional ML (the numerator) only when fit value is smaller than test value
171 if (!fSigned && type==0 &&
172 ((fLimitType==oneSided && fit_favored_mu >= initial_mu_value) ||
173 (fLimitType==oneSidedDiscovery && fit_favored_mu <= initial_mu_value))) {
174 doConditionalFit = false;
175 condML = uncondML;
176 }
177
178 if (doConditionalFit) {
179
180 if (fPrintLevel>1) std::cout << "Do conditional fit " << std::endl;
181
182
183 // cout <<" reestablish snapshot"<<endl;
184 attachedSet->assign(*snap);
185
186
187 // set the POI to constant
188 for (auto *tmpPar : paramsOfInterest) {
189 RooRealVar *tmpParA = dynamic_cast<RooRealVar *>(attachedSet->find(tmpPar->GetName()));
190 if (tmpParA) tmpParA->setConstant();
191 }
192
193
194 // check if there are non-const parameters so it is worth to do the minimization
195 RooArgSet allParams(*attachedSet);
197
198 // in case no nuisance parameters are present
199 // no need to minimize just evaluate the nll
200 if (allParams.empty() ) {
201 // be sure to evaluate with offsets
202 if (fLOffset) RooAbsReal::setHideOffset(false);
203 condML = fNll->getVal();
204 if (fLOffset) RooAbsReal::setHideOffset(true);
205 }
206 else {
207 fNll->clearEvalErrorLog();
208 RooFitResult* result = GetMinNLL();
209 if (result) {
210 condML = result->minNll();
211 statusN = result->status();
212 if( fDetailedOutputEnabled ) {
213 detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitCond_", fDetailedOutputWithErrorsAndPulls);
214 fDetailedOutput->addOwned(*detOutput);
215 delete detOutput;
216 }
217 delete result;
218 }
219 else {
220 return TMath::SignalingNaN(); // this should not really happen
221 }
222 }
223
224 }
225
226 tsw.Stop();
227 double fitTime2 = tsw.CpuTime();
228
229 double pll = 0;
230 if (type != 0) {
231 // for conditional only or unconditional fits
232 // need to compute nll value without the offset
233 if (fLOffset) {
235 pll = fNll->getVal();
236 }
237 else {
238 if (type == 1)
239 pll = uncondML;
240 else if (type == 2)
241 pll = condML;
242 }
243 }
244 else { // type == 0
245 // for standard profile likelihood evaluations
246 pll = condML-uncondML;
247
248 if (fSigned) {
249 if (pll<0.0) {
250 if (fPrintLevel > 0) std::cout << "pll is negative - setting it to zero " << std::endl;
251 pll = 0.0; // bad fit
252 }
253 if (fLimitType==oneSidedDiscovery ? (fit_favored_mu < initial_mu_value)
254 : (fit_favored_mu > initial_mu_value))
255 pll = -pll;
256 }
257 }
258
259 if (fPrintLevel > 0) {
260 std::cout << "EvaluateProfileLikelihood - ";
261 if (type <= 1)
262 std::cout << "mu hat = " << fit_favored_mu << ", uncond ML = " << uncondML;
263 if (type != 1)
264 std::cout << ", cond ML = " << condML;
265 if (type == 0)
266 std::cout << " pll = " << pll;
267 std::cout << " time (create/fit1/2) " << createTime << " , " << fitTime1 << " , " << fitTime2
268 << std::endl;
269 }
270
271
272 // need to restore the values ?
273 attachedSet->assign(*origAttachedSet);
274
275 delete attachedSet;
276 delete origAttachedSet;
277 delete snap;
278
279 if (!reuse) {
280 delete fNll;
281 fNll = 0;
282 }
283
285
286 if(statusN!=0 || statusD!=0) {
287 return -1; // indicate failed fit (WVE is not used anywhere yet)
288 }
289
290 return pll;
291
292 }
293
294////////////////////////////////////////////////////////////////////////////////
295/// find minimum of NLL using RooMinimizer
296
298
299 const auto& config = GetGlobalRooStatsConfig();
300 RooMinimizer minim(*fNll);
301 minim.setStrategy(fStrategy);
302 minim.setEvalErrorWall(config.useEvalErrorWall);
303 //LM: RooMinimizer.setPrintLevel has +1 offset - so subtract here -1 + an extra -1
304 int level = (fPrintLevel == 0) ? -1 : fPrintLevel -2;
305 minim.setPrintLevel(level);
306 minim.setEps(fTolerance);
307 // this causes a memory leak
308 minim.optimizeConst(2);
309 TString minimizer = fMinimizer;
311 if (algorithm == "Migrad") algorithm = "Minimize"; // prefer to use Minimize instead of Migrad
312 int status;
313 for (int tries = 1, maxtries = 4; tries <= maxtries; ++tries) {
314 status = minim.minimize(minimizer,algorithm);
315 if (status%1000 == 0) { // ignore erros from Improve
316 break;
317 } else if (tries < maxtries) {
318 cout << " ----> Doing a re-scan first" << endl;
319 minim.minimize(minimizer,"Scan");
320 if (tries == 2) {
321 if (fStrategy == 0 ) {
322 cout << " ----> trying with strategy = 1" << endl;;
323 minim.setStrategy(1);
324 }
325 else
326 tries++; // skip this trial if strategy is already 1
327 }
328 if (tries == 3) {
329 cout << " ----> trying with improve" << endl;;
330 minimizer = "Minuit";
331 algorithm = "migradimproved";
332 }
333 }
334 }
335
336 //how to get cov quality faster?
337 return minim.save();
338 //minim.optimizeConst(false);
339}
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
static const std::string & DefaultMinimizerAlgo()
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
bool empty() const
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
RooAbsArg * find(const char *name) const
Find object with given name in list.
void Print(Option_t *options=nullptr) const override
This method must be overridden when a class wants to print itself.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:61
void setConstant(bool value=true)
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:104
static void setHideOffset(bool flag)
Definition: RooAbsReal.cxx:115
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:179
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooMinimizer is a wrapper class around ROOT::Fit:Fitter that provides a seamless interface between th...
Definition: RooMinimizer.h:43
void optimizeConst(int flag)
If flag is true, perform constant term optimization on function being minimized.
RooFitResult * save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
void setEvalErrorWall(bool flag)
Definition: RooMinimizer.h:92
void setEps(double eps)
Change MINUIT epsilon.
void setPrintLevel(int newLevel)
Change the MINUIT internal printing level.
int minimize(const char *type, const char *alg=nullptr)
Minimise the function passed in the constructor.
void setStrategy(int 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:40
static RooArgSet * GetAsArgSet(RooFitResult *result, TString prefix="", bool withErrorsAndPulls=false)
Translate the given fit result to a RooArgSet in a generic way.
virtual double EvaluateProfileLikelihood(int type, RooAbsData &data, RooArgSet &paramsOfInterest)
evaluate the profile likelihood ratio (type = 0) or the minimum of likelihood (type=1) or the conditi...
RooFitResult * GetMinNLL()
find minimum of NLL using RooMinimizer
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Stopwatch class.
Definition: TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
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
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Basic string class.
Definition: TString.h:136
RooCmdArg Offset(std::string const &mode)
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg CloneData(bool flag)
RooCmdArg ConditionalObservables(Args_t &&... argsOrArgSet)
Create a RooCmdArg to declare conditional observables.
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:60
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:67
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)
useful function to print in one line the content of a set with their values
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition: TMath.h:908