Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
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 std::unique_ptr<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 std::unique_ptr<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 std::unique_ptr<RooAbsArg> profileOwner{nll->createProfile(allButNuisance)};
102 auto profile = dynamic_cast<RooProfileLL*>(profileOwner.get());
103 // set minimier options
105 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
106
107 // Hack to extract a RooFitResult
108 if (fStoreFitInfo) {
109 std::unique_ptr<RooFitResult> result {profile->minimizer()->save()};
110 std::unique_ptr<RooArgSet> detOutput {DetailedOutputAggregator::GetAsArgSet(result.get(), "fitNull_")};
111 fFitInfo->addOwned(*detOutput);
112 }
113
115
116 // set in test statistics conditional and global observables
117 // (needed to get correct model likelihood)
118 TestStatistic * testStatistic = nullptr;
119 auto testStatSampler = GetTestStatSampler();
120 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
121 if (testStatistic) {
122 testStatistic->SetConditionalObservables(conditionalObs);
123 testStatistic->SetGlobalObservables(globalObs);
124 }
125
126 }
127
128 // add nuisance parameters to parameter point
130 parameterPoint->add(*fNullModel->GetNuisanceParameters());
131
132
133 // ***** ToyMCSampler specific *******
134
135 // check whether TestStatSampler is a ToyMCSampler
136 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
137 if(toymcs) {
138 oocoutI(nullptr,InputArguments) << "Using a ToyMCSampler. Now configuring for Null." << endl;
139
140 // variable number of toys
141 if(fNToysNull >= 0) toymcs->SetNToys(fNToysNull);
142
143 // set the global observables to be generated by the ToyMCSampler
145
146 // adaptive sampling
147 if(fNToysNullTail) {
148 oocoutI(nullptr,InputArguments) << "Adaptive Sampling" << endl;
149 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
150 toymcs->SetToysRightTail(fNToysNullTail, obsTestStat);
151 }else{
152 toymcs->SetToysLeftTail(fNToysNullTail, obsTestStat);
153 }
154 }else{
155 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
156 }
157
159 }
160
161 return 0;
162}
163
164////////////////////////////////////////////////////////////////////////////////
165
166int FrequentistCalculator::PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const {
167
168 // ****** any TestStatSampler ********
169
170 // create profile keeping everything but nuisance parameters fixed
171 std::unique_ptr<RooArgSet> allParams{fAltModel->GetPdf()->getParameters(*fData)};
172 RemoveConstantParameters(&*allParams);
173
174 bool doProfile = true;
175 RooArgSet allButNuisance(*allParams);
177 allButNuisance.remove(*fAltModel->GetNuisanceParameters());
178 if( fConditionalMLEsAlt ) {
179 oocoutI(nullptr,InputArguments) << "Using given conditional MLEs for Alt." << endl;
180 allParams->assign(*fConditionalMLEsAlt);
181 // LM: fConditionalMLEsAlt must be nuisance parameters otherwise an error message will be printed
182 allButNuisance.add( *fConditionalMLEsAlt );
185 remain.remove(*fConditionalMLEsAlt,true,true);
186 if( remain.empty() ) doProfile = false;
187 }
188 }
189 }else{
190 doProfile = false;
191 }
192 if (doProfile) {
193 oocoutI(nullptr,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
196
197 RooArgSet conditionalObs;
199 RooArgSet globalObs;
201
202 const auto& config = GetGlobalRooStatsConfig();
203 std::unique_ptr<RooAbsReal> nll{fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(false), RooFit::Constrain(*allParams),
204 RooFit::GlobalObservables(globalObs),
205 RooFit::ConditionalObservables(conditionalObs),
206 RooFit::Offset(config.useLikelihoodOffset))};
207
208 std::unique_ptr<RooAbsReal> profileOwner{nll->createProfile(allButNuisance)};
209 auto profile = dynamic_cast<RooProfileLL*>(profileOwner.get());
210 // set minimizer options
211 profile->minimizer()->setPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()-1); // use -1 to make more silent
212 profile->getVal(); // this will do fit and set nuisance parameters to profiled values
213
214 // Hack to extract a RooFitResult
215 if (fStoreFitInfo) {
216 std::unique_ptr<RooFitResult> result {profile->minimizer()->save()};
217 std::unique_ptr<RooArgSet> detOutput {DetailedOutputAggregator::GetAsArgSet(result.get(), "fitAlt_")};
218 fFitInfo->addOwned(*detOutput);
219 }
220
222
223 // set in test statistics conditional and global observables
224 // (needed to get correct model likelihood)
225 TestStatistic * testStatistic = nullptr;
226 auto testStatSampler = GetTestStatSampler();
227 if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
228 if (testStatistic) {
229 testStatistic->SetConditionalObservables(conditionalObs);
230 testStatistic->SetGlobalObservables(globalObs);
231 }
232
233 }
234
235 // add nuisance parameters to parameter point
237 parameterPoint->add(*fAltModel->GetNuisanceParameters());
238
239 // ***** ToyMCSampler specific *******
240
241 // check whether TestStatSampler is a ToyMCSampler
242 ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
243 if(toymcs) {
244 oocoutI(nullptr,InputArguments) << "Using a ToyMCSampler. Now configuring for Alt." << endl;
245
246 // variable number of toys
247 if(fNToysAlt >= 0) toymcs->SetNToys(fNToysAlt);
248
249 // set the global observables to be generated by the ToyMCSampler
251
252 // adaptive sampling
253 if(fNToysAltTail) {
254 oocoutI(nullptr,InputArguments) << "Adaptive Sampling" << endl;
255 if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
256 toymcs->SetToysLeftTail(fNToysAltTail, obsTestStat);
257 }else{
258 toymcs->SetToysRightTail(fNToysAltTail, obsTestStat);
259 }
260 }else{
261 toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
262 }
263
264 }
265
266 return 0;
267}
#define oocoutI(o, a)
#define ClassImp(name)
Definition Rtypes.h:377
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
RooFit::OwningPtr< 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...
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
RooFit::OwningPtr< RooAbsReal > createNLL(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Construct representation of -log(L) of PDF with given dataset.
Definition RooAbsPdf.h:162
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
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...
RooMinimizer * minimizer()
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)
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
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)
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
virtual void SetConditionalObservables(const RooArgSet &)
interface to set conditional observables. If a test statistics needs them it will re-implement this f...
virtual void SetGlobalObservables(const RooArgSet &)
interface to set global observables. If a test statistics needs them it will re-implement this functi...
ToyMCSampler is an implementation of the TestStatSampler interface.
virtual void SetNToys(const Int_t ntoy)
void SetToysBothTails(double toys, double low_threshold, double high_threshold)
void SetToysRightTail(double toys, double threshold)
void SetToysLeftTail(double toys, double threshold)
void SetGlobalObservables(const RooArgSet &o) override
specify the conditional observables
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.
Namespace for the RooStats classes.
Definition Asimov.h:19
void RemoveConstantParameters(RooArgSet *set)
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...