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 
12 #include "RooStats/ToyMCSampler.h"
13 #include "RooStats/RooStatsUtils.h"
15 #include "RooMinimizer.h"
16 #include "RooMinuit.h"
17 #include "RooProfileLL.h"
18 
19 /** \class RooStats::FrequentistCalculator
20  \ingroup Roostats
21 
22 Does a frequentist hypothesis test.
23 
24 Hypothesis Test Calculator using a full frequentist procedure for sampling the
25 test statistic distribution.
26 The nuisance parameters are fixed to their MLEs.
27 The use of ToyMCSampler as the TestStatSampler is assumed.
28 
29 */
30 
32 
33 using namespace RooStats;
34 using namespace std;
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 
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 
55 int 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);
71  if( fNullModel->GetNuisanceParameters() ) {
72  allButNuisance.remove(*fNullModel->GetNuisanceParameters());
73  if( fConditionalMLEsNull ) {
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 );
78  if (fNullModel->GetNuisanceParameters()) {
79  RooArgSet remain(*fNullModel->GetNuisanceParameters());
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;
93  if (fNullModel->GetConditionalObservables()) conditionalObs.add(*fNullModel->GetConditionalObservables());
94  RooArgSet globalObs;
95  if (fNullModel->GetGlobalObservables()) globalObs.add(*fNullModel->GetGlobalObservables());
96 
97  auto& config = GetGlobalRooStatsConfig();
98  RooAbsReal* nll = fNullModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams),
99  RooFit::GlobalObservables(globalObs),
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
131  if(fNullModel->GetNuisanceParameters())
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
148  toymcs->SetGlobalObservables(*fNullModel->GetGlobalObservables());
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 
162  GetNullModel()->LoadSnapshot();
163  }
164 
165  return 0;
166 }
167 
168 ////////////////////////////////////////////////////////////////////////////////
169 
170 int 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);
180  if( fAltModel->GetNuisanceParameters() ) {
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 );
187  if (fAltModel->GetNuisanceParameters()) {
188  RooArgSet remain(*fAltModel->GetNuisanceParameters());
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;
202  if (fAltModel->GetConditionalObservables()) conditionalObs.add(*fAltModel->GetConditionalObservables());
203  RooArgSet globalObs;
204  if (fAltModel->GetGlobalObservables()) globalObs.add(*fAltModel->GetGlobalObservables());
205 
206  const auto& config = GetGlobalRooStatsConfig();
207  RooAbsReal* nll = fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams),
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
241  if(fAltModel->GetNuisanceParameters())
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
257  toymcs->SetGlobalObservables(*fAltModel->GetGlobalObservables());
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 }
RooStats::FrequentistCalculator::PreAltHook
int PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const
configure TestStatSampler for the Alt run
Definition: FrequentistCalculator.cxx:170
RooStats::ToyMCSampler::SetToysLeftTail
void SetToysLeftTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:224
RooMinuit.h
RooStats::DetailedOutputAggregator::GetAsArgSet
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.
Definition: DetailedOutputAggregator.cxx:50
RooStats::FrequentistCalculator::PostHook
void PostHook() const
Definition: FrequentistCalculator.cxx:50
RooAbsData
Definition: RooAbsData.h:46
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
oocoutI
#define oocoutI(o, a)
Definition: RooMsgService.h:45
RooStats::FrequentistCalculator::PreHook
void PreHook() const
Definition: FrequentistCalculator.cxx:38
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
DetailedOutputAggregator.h
RooProfileLL
Definition: RooProfileLL.h:26
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
ToyMCSampler.h
RooFit::Offset
RooCmdArg Offset(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:212
RooAbsCollection::remove
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
Definition: RooAbsCollection.cxx:584
RooStats::TestStatistic
Definition: TestStatistic.h:31
RooStats::ToyMCSampler::SetToysRightTail
void SetToysRightTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:229
RooArgSet::add
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
RooAbsReal
Definition: RooAbsReal.h:61
RooFit::MsgLevel
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
RooFitResult
Definition: RooFitResult.h:40
RooProfileLL::minimizer
MINIMIZER * minimizer()
Definition: RooProfileLL.h:38
RooFit::GlobalObservables
RooCmdArg GlobalObservables(const RooArgSet &globs)
Definition: RooGlobalFunc.cxx:201
RooStats::FrequentistCalculator
Definition: FrequentistCalculator.h:31
FrequentistCalculator.h
RooStats::ToyMCSampler::SetNToys
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:154
RooFit::FATAL
@ FATAL
Definition: RooGlobalFunc.h:65
RooFit::Constrain
RooCmdArg Constrain(const RooArgSet &params)
Definition: RooGlobalFunc.cxx:200
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooMsgService::setGlobalKillBelow
void setGlobalKillBelow(RooFit::MsgLevel level)
Definition: RooMsgService.h:160
RooStats::ToyMCSampler::SetToysBothTails
void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold)
Definition: ToyMCSampler.h:234
RooStats::TestStatistic::SetGlobalObservables
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
RooFit::ConditionalObservables
RooCmdArg ConditionalObservables(const RooArgSet &set)
Definition: RooGlobalFunc.cxx:196
RooAbsReal::createProfile
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function.
Definition: RooAbsReal.cxx:516
RooStatsUtils.h
RooMinimizer.h
RooStats::FrequentistCalculator::PreNullHook
int PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const
configure TestStatSampler for the Null run
Definition: FrequentistCalculator.cxx:55
RooStats
Definition: Asimov.h:19
RooStats::ToyMCSampler
Definition: ToyMCSampler.h:79
TObject
Definition: TObject.h:37
RooStats::RemoveConstantParameters
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:69
RooFit::CloneData
RooCmdArg CloneData(Bool_t flag)
Definition: RooGlobalFunc.cxx:209
RooStats::GetGlobalRooStatsConfig
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
Definition: RooStatsUtils.cxx:33
RooStats::TestStatistic::SetConditionalObservables
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
RooMsgService::instance
static RooMsgService & instance()
Return reference to singleton instance.
Definition: RooMsgService.cxx:363
RooMsgService::globalKillBelow
RooFit::MsgLevel globalKillBelow() const
Definition: RooMsgService.h:161
RooStats::ToyMCSampler::SetGlobalObservables
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:183
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooProfileLL.h
RooArgSet
Definition: RooArgSet.h:28