Logo ROOT  
Reference Guide
FrequentistCalculator.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id: FrequentistCalculator.cxx 37084 2010-11-29 21:37:13Z moneta $
2// Author: Kyle Cranmer, Sven Kreiss 23/05/10
3/*************************************************************************
4 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
15#include "RooMinimizer.h"
16#include "RooMinuit.h"
17#include "RooProfileLL.h"
18
19/** \class RooStats::FrequentistCalculator
20 \ingroup Roostats
21
22Does a frequentist hypothesis test.
23
24Hypothesis Test Calculator using a full frequentist procedure for sampling the
25test statistic distribution.
26The nuisance parameters are fixed to their MLEs.
27The use of ToyMCSampler as the TestStatSampler is assumed.
28
29*/
30
32
33using namespace RooStats;
34using namespace std;
35
36////////////////////////////////////////////////////////////////////////////////
37
38void FrequentistCalculator::PreHook() const {
39 if (fFitInfo != NULL) {
40 delete fFitInfo;
41 fFitInfo = NULL;
42 }
43 if (fStoreFitInfo) {
44 fFitInfo = new RooArgSet();
45 }
46}
47
48////////////////////////////////////////////////////////////////////////////////
49
51}
52
53////////////////////////////////////////////////////////////////////////////////
54
55int FrequentistCalculator::PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const {
56
57 // ****** any TestStatSampler ********
58
59 // create profile keeping everything but nuisance parameters fixed
60 RooArgSet * allParams = fNullModel->GetPdf()->getParameters(*fData);
61 RemoveConstantParameters(allParams);
62
63 // note: making nll or profile class variables can only be done in the constructor
64 // as all other hooks are const (which has to be because GetHypoTest is const). However,
65 // when setting it only in constructor, they would have to be changed every time SetNullModel
66 // or SetAltModel is called. Simply put, converting them into class variables breaks
67 // encapsulation.
68
69 bool doProfile = true;
70 RooArgSet allButNuisance(*allParams);
72 allButNuisance.remove(*fNullModel->GetNuisanceParameters());
74 oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Null." << endl;
75 *allParams = *fConditionalMLEsNull;
76 // LM: fConditionalMLEsNull must be nuisance parameters otherwise an error message will be printed
77 allButNuisance.add( *fConditionalMLEsNull );
80 remain.remove(*fConditionalMLEsNull,true,true);
81 if( remain.getSize() == 0 ) doProfile = false;
82 }
83 }
84 }else{
85 doProfile = false;
86 }
87 if (doProfile) {
88 oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Null." << endl;
91
92 RooArgSet conditionalObs;
94 RooArgSet globalObs;
96
97 auto& config = GetGlobalRooStatsConfig();
100 RooFit::ConditionalObservables(conditionalObs),
101 RooFit::Offset(config.useLikelihoodOffset));
102 RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
103 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
104
105 // Hack to extract a RooFitResult
106 if (fStoreFitInfo) {
107 RooFitResult *result = profile->minimizer()->save();
108 RooArgSet * detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitNull_");
109 fFitInfo->addOwned(*detOutput);
110 delete detOutput;
111 delete result;
112 }
113
114 delete profile;
115 delete nll;
117
118 // set in test statistics conditional and global observables
119 // (needed to get correct model likelihood)
120 TestStatistic * testStatistic = nullptr;
121 auto testStatSampler = GetTestStatSampler();
122 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
123 if (testStatistic) {
124 testStatistic->SetConditionalObservables(conditionalObs);
125 testStatistic->SetGlobalObservables(globalObs);
126 }
127
128 }
129
130 // add nuisance parameters to parameter point
132 parameterPoint->add(*fNullModel->GetNuisanceParameters());
133
134 delete allParams;
135
136
137 // ***** ToyMCSampler specific *******
138
139 // check whether TestStatSampler is a ToyMCSampler
140 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
141 if(toymcs) {
142 oocoutI((TObject*)0,InputArguments) << "Using a ToyMCSampler. Now configuring for Null." << endl;
143
144 // variable number of toys
145 if(fNToysNull >= 0) toymcs->SetNToys(fNToysNull);
146
147 // set the global observables to be generated by the ToyMCSampler
149
150 // adaptive sampling
151 if(fNToysNullTail) {
152 oocoutI((TObject*)0,InputArguments) << "Adaptive Sampling" << endl;
153 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
154 toymcs->SetToysRightTail(fNToysNullTail, obsTestStat);
155 }else{
156 toymcs->SetToysLeftTail(fNToysNullTail, obsTestStat);
157 }
158 }else{
159 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
160 }
161
163 }
164
165 return 0;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169
170int FrequentistCalculator::PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const {
171
172 // ****** any TestStatSampler ********
173
174 // create profile keeping everything but nuisance parameters fixed
175 RooArgSet * allParams = fAltModel->GetPdf()->getParameters(*fData);
176 RemoveConstantParameters(allParams);
177
178 bool doProfile = true;
179 RooArgSet allButNuisance(*allParams);
181 allButNuisance.remove(*fAltModel->GetNuisanceParameters());
182 if( fConditionalMLEsAlt ) {
183 oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Alt." << endl;
184 *allParams = *fConditionalMLEsAlt;
185 // LM: fConditionalMLEsAlt must be nuisance parameters otherwise an error message will be printed
186 allButNuisance.add( *fConditionalMLEsAlt );
189 remain.remove(*fConditionalMLEsAlt,true,true);
190 if( remain.getSize() == 0 ) doProfile = false;
191 }
192 }
193 }else{
194 doProfile = false;
195 }
196 if (doProfile) {
197 oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
200
201 RooArgSet conditionalObs;
203 RooArgSet globalObs;
205
206 const auto& config = GetGlobalRooStatsConfig();
208 RooFit::GlobalObservables(globalObs),
209 RooFit::ConditionalObservables(conditionalObs),
210 RooFit::Offset(config.useLikelihoodOffset));
211
212 RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
213 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
214
215 // Hack to extract a RooFitResult
216 if (fStoreFitInfo) {
217 RooFitResult *result = profile->minimizer()->save();
218 RooArgSet * detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitAlt_");
219 fFitInfo->addOwned(*detOutput);
220 delete detOutput;
221 delete result;
222 }
223
224 delete profile;
225 delete nll;
227
228 // set in test statistics conditional and global observables
229 // (needed to get correct model likelihood)
230 TestStatistic * testStatistic = nullptr;
231 auto testStatSampler = GetTestStatSampler();
232 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
233 if (testStatistic) {
234 testStatistic->SetConditionalObservables(conditionalObs);
235 testStatistic->SetGlobalObservables(globalObs);
236 }
237
238 }
239
240 // add nuisance parameters to parameter point
242 parameterPoint->add(*fAltModel->GetNuisanceParameters());
243
244 delete allParams;
245
246 // ***** ToyMCSampler specific *******
247
248 // check whether TestStatSampler is a ToyMCSampler
249 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
250 if(toymcs) {
251 oocoutI((TObject*)0,InputArguments) << "Using a ToyMCSampler. Now configuring for Alt." << endl;
252
253 // variable number of toys
254 if(fNToysAlt >= 0) toymcs->SetNToys(fNToysAlt);
255
256 // set the global observables to be generated by the ToyMCSampler
258
259 // adaptive sampling
260 if(fNToysAltTail) {
261 oocoutI((TObject*)0,InputArguments) << "Adaptive Sampling" << endl;
262 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
263 toymcs->SetToysLeftTail(fNToysAltTail, obsTestStat);
264 }else{
265 toymcs->SetToysRightTail(fNToysAltTail, obsTestStat);
266 }
267 }else{
268 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
269 }
270
271 }
272
273 return 0;
274}
#define oocoutI(o, a)
Definition: RooMsgService.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:90
#define ClassImp(name)
Definition: Rtypes.h:361
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't match any of...
Definition: RooAbsArg.cxx:544
Int_t getSize() const
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:918
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:515
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
virtual Bool_t addOwned(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:92
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
static RooMsgService & instance()
Return reference to singleton instance.
void setGlobalKillBelow(RooFit::MsgLevel level)
RooFit::MsgLevel globalKillBelow() const
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:26
MINIMIZER * minimizer()
Definition: RooProfileLL.h:38
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.
Does a frequentist hypothesis test.
int PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const
configure TestStatSampler for the Null run
int PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const
configure TestStatSampler for the Alt run
const ModelConfig * GetNullModel(void) const
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
const RooArgSet * GetConditionalObservables() const
get RooArgSet for conditional observables (return NULL if not existing)
Definition: ModelConfig.h:252
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:255
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
void LoadSnapshot() const
load the snapshot from ws if it exists
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31
virtual void SetConditionalObservables(const RooArgSet &)
interface to set conditional observables. If a test statistics needs them it will re-implement this f...
Definition: TestStatistic.h:51
virtual void SetGlobalObservables(const RooArgSet &)
interface to set global observables. If a test statistics needs them it will re-implement this functi...
Definition: TestStatistic.h:54
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:72
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:147
void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold)
Definition: ToyMCSampler.h:226
void SetToysLeftTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:216
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:175
void SetToysRightTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:221
Mother of all ROOT objects.
Definition: TObject.h:37
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(const RooArgSet &globs)
RooCmdArg CloneData(Bool_t flag)
RooCmdArg ConditionalObservables(const RooArgSet &set)
RooCmdArg Offset(Bool_t flag=kTRUE)
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...