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 "RooProfileLL.h"
17
18/** \class RooStats::FrequentistCalculator
19 \ingroup Roostats
20
21Does a frequentist hypothesis test.
22
23Hypothesis Test Calculator using a full frequentist procedure for sampling the
24test statistic distribution.
25The nuisance parameters are fixed to their MLEs.
26The use of ToyMCSampler as the TestStatSampler is assumed.
27
28*/
29
31
32using namespace RooStats;
33using namespace std;
34
35////////////////////////////////////////////////////////////////////////////////
36
37void FrequentistCalculator::PreHook() const {
38 if (fFitInfo != nullptr) {
39 delete fFitInfo;
40 fFitInfo = nullptr;
41 }
42 if (fStoreFitInfo) {
43 fFitInfo = new RooArgSet();
44 }
45}
46
47////////////////////////////////////////////////////////////////////////////////
48
50}
51
52////////////////////////////////////////////////////////////////////////////////
53
54int FrequentistCalculator::PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const {
55
56 // ****** any TestStatSampler ********
57
58 // create profile keeping everything but nuisance parameters fixed
59 RooArgSet * allParams = fNullModel->GetPdf()->getParameters(*fData);
60 RemoveConstantParameters(allParams);
61
62 // note: making nll or profile class variables can only be done in the constructor
63 // as all other hooks are const (which has to be because GetHypoTest is const). However,
64 // when setting it only in constructor, they would have to be changed every time SetNullModel
65 // or SetAltModel is called. Simply put, converting them into class variables breaks
66 // encapsulation.
67
68 bool doProfile = true;
69 RooArgSet allButNuisance(*allParams);
71 allButNuisance.remove(*fNullModel->GetNuisanceParameters());
73 oocoutI(nullptr,InputArguments) << "Using given conditional MLEs for Null." << endl;
74 allParams->assign(*fConditionalMLEsNull);
75 // LM: fConditionalMLEsNull must be nuisance parameters otherwise an error message will be printed
76 allButNuisance.add( *fConditionalMLEsNull );
79 remain.remove(*fConditionalMLEsNull,true,true);
80 if( remain.empty() ) doProfile = false;
81 }
82 }
83 }else{
84 doProfile = false;
85 }
86 if (doProfile) {
87 oocoutI(nullptr,InputArguments) << "Profiling conditional MLEs for Null." << endl;
90
91 RooArgSet conditionalObs;
93 RooArgSet globalObs;
95
96 auto& config = GetGlobalRooStatsConfig();
97 RooAbsReal* nll = fNullModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(false), RooFit::Constrain(*allParams),
99 RooFit::ConditionalObservables(conditionalObs),
100 RooFit::Offset(config.useLikelihoodOffset));
101 RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
102 // set minimier options
104 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
105
106 // Hack to extract a RooFitResult
107 if (fStoreFitInfo) {
108 RooFitResult *result = profile->minimizer()->save();
109 RooArgSet * detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitNull_");
110 fFitInfo->addOwned(*detOutput);
111 delete detOutput;
112 delete result;
113 }
114
115 delete profile;
116 delete nll;
118
119 // set in test statistics conditional and global observables
120 // (needed to get correct model likelihood)
121 TestStatistic * testStatistic = nullptr;
122 auto testStatSampler = GetTestStatSampler();
123 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
124 if (testStatistic) {
125 testStatistic->SetConditionalObservables(conditionalObs);
126 testStatistic->SetGlobalObservables(globalObs);
127 }
128
129 }
130
131 // add nuisance parameters to parameter point
133 parameterPoint->add(*fNullModel->GetNuisanceParameters());
134
135 delete allParams;
136
137
138 // ***** ToyMCSampler specific *******
139
140 // check whether TestStatSampler is a ToyMCSampler
141 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
142 if(toymcs) {
143 oocoutI(nullptr,InputArguments) << "Using a ToyMCSampler. Now configuring for Null." << endl;
144
145 // variable number of toys
146 if(fNToysNull >= 0) toymcs->SetNToys(fNToysNull);
147
148 // set the global observables to be generated by the ToyMCSampler
150
151 // adaptive sampling
152 if(fNToysNullTail) {
153 oocoutI(nullptr,InputArguments) << "Adaptive Sampling" << endl;
154 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
155 toymcs->SetToysRightTail(fNToysNullTail, obsTestStat);
156 }else{
157 toymcs->SetToysLeftTail(fNToysNullTail, obsTestStat);
158 }
159 }else{
160 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
161 }
162
164 }
165
166 return 0;
167}
168
169////////////////////////////////////////////////////////////////////////////////
170
171int FrequentistCalculator::PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const {
172
173 // ****** any TestStatSampler ********
174
175 // create profile keeping everything but nuisance parameters fixed
176 RooArgSet * allParams = fAltModel->GetPdf()->getParameters(*fData);
177 RemoveConstantParameters(allParams);
178
179 bool doProfile = true;
180 RooArgSet allButNuisance(*allParams);
182 allButNuisance.remove(*fAltModel->GetNuisanceParameters());
183 if( fConditionalMLEsAlt ) {
184 oocoutI(nullptr,InputArguments) << "Using given conditional MLEs for Alt." << endl;
185 allParams->assign(*fConditionalMLEsAlt);
186 // LM: fConditionalMLEsAlt must be nuisance parameters otherwise an error message will be printed
187 allButNuisance.add( *fConditionalMLEsAlt );
190 remain.remove(*fConditionalMLEsAlt,true,true);
191 if( remain.empty() ) doProfile = false;
192 }
193 }
194 }else{
195 doProfile = false;
196 }
197 if (doProfile) {
198 oocoutI(nullptr,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
201
202 RooArgSet conditionalObs;
204 RooArgSet globalObs;
206
207 const auto& config = GetGlobalRooStatsConfig();
208 RooAbsReal* nll = fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(false), RooFit::Constrain(*allParams),
209 RooFit::GlobalObservables(globalObs),
210 RooFit::ConditionalObservables(conditionalObs),
211 RooFit::Offset(config.useLikelihoodOffset));
212
213 RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
214 // set minimizer options
215 profile->minimizer()->setPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()-1); // use -1 to make more silent
216 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
217
218 // Hack to extract a RooFitResult
219 if (fStoreFitInfo) {
220 RooFitResult *result = profile->minimizer()->save();
222 fFitInfo->addOwned(*detOutput);
223 delete detOutput;
224 delete result;
225 }
226
227 delete profile;
228 delete nll;
230
231 // set in test statistics conditional and global observables
232 // (needed to get correct model likelihood)
233 TestStatistic * testStatistic = nullptr;
234 auto testStatSampler = GetTestStatSampler();
235 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
236 if (testStatistic) {
237 testStatistic->SetConditionalObservables(conditionalObs);
238 testStatistic->SetGlobalObservables(globalObs);
239 }
240
241 }
242
243 // add nuisance parameters to parameter point
245 parameterPoint->add(*fAltModel->GetNuisanceParameters());
246
247 delete allParams;
248
249 // ***** ToyMCSampler specific *******
250
251 // check whether TestStatSampler is a ToyMCSampler
252 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
253 if(toymcs) {
254 oocoutI(nullptr,InputArguments) << "Using a ToyMCSampler. Now configuring for Alt." << endl;
255
256 // variable number of toys
257 if(fNToysAlt >= 0) toymcs->SetNToys(fNToysAlt);
258
259 // set the global observables to be generated by the ToyMCSampler
261
262 // adaptive sampling
263 if(fNToysAltTail) {
264 oocoutI(nullptr,InputArguments) << "Adaptive Sampling" << endl;
265 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
266 toymcs->SetToysLeftTail(fNToysAltTail, obsTestStat);
267 }else{
268 toymcs->SetToysRightTail(fNToysAltTail, obsTestStat);
269 }
270 }else{
271 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
272 }
273
274 }
275
276 return 0;
277}
#define oocoutI(o, a)
Definition: RooMsgService.h:49
#define ClassImp(name)
Definition: Rtypes.h:375
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
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) 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:541
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
bool empty() const
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
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...
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:994
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:480
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooFitResult * save(const char *name=nullptr, const char *title=nullptr)
Save and return a RooFitResult snapshot of current minimizer status.
void setPrintLevel(int newLevel)
Change the MINUIT internal printing level.
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:22
RooMinimizer * minimizer()
Definition: RooProfileLL.h:33
static RooArgSet * GetAsArgSet(RooFitResult *result, TString prefix="", bool withErrorsAndPulls=false)
Translate the given fit result to a RooArgSet in a generic way.
Does a frequentist hypothesis test.
int PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const override
configure TestStatSampler for the Null run
int PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const override
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 nullptr if not existing)
Definition: ModelConfig.h:249
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
Definition: ModelConfig.h:252
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
Definition: ModelConfig.h:237
void LoadSnapshot() const
load the snapshot from ws if it exists
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
Definition: ModelConfig.h:231
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:67
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:141
void SetToysBothTails(double toys, double low_threshold, double high_threshold)
Definition: ToyMCSampler.h:220
void SetToysRightTail(double toys, double threshold)
Definition: ToyMCSampler.h:215
void SetToysLeftTail(double toys, double threshold)
Definition: ToyMCSampler.h:210
void SetGlobalObservables(const RooArgSet &o) override
specify the conditional observables
Definition: ToyMCSampler.h:169
RooCmdArg Constrain(const RooArgSet &params)
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
RooCmdArg Offset(bool flag=true)
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:59
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
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...