Logo ROOT   6.12/07
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  RooAbsReal* nll = fNullModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams),
98  RooFit::GlobalObservables(globalObs),
100  RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
101  profile->getVal(); // this will do fit and set nuisance parameters to profiled values
102 
103  // Hack to extract a RooFitResult
104  if (fStoreFitInfo) {
105  RooFitResult *result = profile->minimizer()->save();
106  RooArgSet * detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitNull_");
107  fFitInfo->addOwned(*detOutput);
108  delete detOutput;
109  delete result;
110  }
111 
112  delete profile;
113  delete nll;
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
129  if(fNullModel->GetNuisanceParameters())
130  parameterPoint->add(*fNullModel->GetNuisanceParameters());
131 
132  delete allParams;
133 
134 
135  // ***** ToyMCSampler specific *******
136 
137  // check whether TestStatSampler is a ToyMCSampler
138  ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
139  if(toymcs) {
140  oocoutI((TObject*)0,InputArguments) << "Using a ToyMCSampler. Now configuring for Null." << endl;
141 
142  // variable number of toys
143  if(fNToysNull >= 0) toymcs->SetNToys(fNToysNull);
144 
145  // set the global observables to be generated by the ToyMCSampler
146  toymcs->SetGlobalObservables(*fNullModel->GetGlobalObservables());
147 
148  // adaptive sampling
149  if(fNToysNullTail) {
150  oocoutI((TObject*)0,InputArguments) << "Adaptive Sampling" << endl;
151  if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
152  toymcs->SetToysRightTail(fNToysNullTail, obsTestStat);
153  }else{
154  toymcs->SetToysLeftTail(fNToysNullTail, obsTestStat);
155  }
156  }else{
157  toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
158  }
159 
160  GetNullModel()->LoadSnapshot();
161  }
162 
163  return 0;
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 
168 int FrequentistCalculator::PreAltHook(RooArgSet *parameterPoint, double obsTestStat) const {
169 
170  // ****** any TestStatSampler ********
171 
172  // create profile keeping everything but nuisance parameters fixed
173  RooArgSet * allParams = fAltModel->GetPdf()->getParameters(*fData);
174  RemoveConstantParameters(allParams);
175 
176  bool doProfile = true;
177  RooArgSet allButNuisance(*allParams);
178  if( fAltModel->GetNuisanceParameters() ) {
179  allButNuisance.remove(*fAltModel->GetNuisanceParameters());
180  if( fConditionalMLEsAlt ) {
181  oocoutI((TObject*)0,InputArguments) << "Using given conditional MLEs for Alt." << endl;
182  *allParams = *fConditionalMLEsAlt;
183  // LM: fConditionalMLEsAlt must be nuisance parameters otherwise an error message will be printed
184  allButNuisance.add( *fConditionalMLEsAlt );
185  if (fAltModel->GetNuisanceParameters()) {
186  RooArgSet remain(*fAltModel->GetNuisanceParameters());
187  remain.remove(*fConditionalMLEsAlt,true,true);
188  if( remain.getSize() == 0 ) doProfile = false;
189  }
190  }
191  }else{
192  doProfile = false;
193  }
194  if (doProfile) {
195  oocoutI((TObject*)0,InputArguments) << "Profiling conditional MLEs for Alt." << endl;
198 
199  RooArgSet conditionalObs;
200  if (fAltModel->GetConditionalObservables()) conditionalObs.add(*fAltModel->GetConditionalObservables());
201  RooArgSet globalObs;
202  if (fAltModel->GetGlobalObservables()) globalObs.add(*fAltModel->GetGlobalObservables());
203 
204 
205  RooAbsReal* nll = fAltModel->GetPdf()->createNLL(*const_cast<RooAbsData*>(fData), RooFit::CloneData(kFALSE), RooFit::Constrain(*allParams),
206  RooFit::GlobalObservables(globalObs),
208 
209  RooProfileLL* profile = dynamic_cast<RooProfileLL*>(nll->createProfile(allButNuisance));
210  profile->getVal(); // this will do fit and set nuisance parameters to profiled values
211 
212  // Hack to extract a RooFitResult
213  if (fStoreFitInfo) {
214  RooFitResult *result = profile->minimizer()->save();
215  RooArgSet * detOutput = DetailedOutputAggregator::GetAsArgSet(result, "fitAlt_");
216  fFitInfo->addOwned(*detOutput);
217  delete detOutput;
218  delete result;
219  }
220 
221  delete profile;
222  delete nll;
224 
225  // set in test statistics conditional and global observables
226  // (needed to get correct model likelihood)
227  TestStatistic * testStatistic = nullptr;
228  auto testStatSampler = GetTestStatSampler();
229  if (testStatSampler) testStatistic = testStatSampler->GetTestStatistic();
230  if (testStatistic) {
231  testStatistic->SetConditionalObservables(conditionalObs);
232  testStatistic->SetGlobalObservables(globalObs);
233  }
234 
235  }
236 
237  // add nuisance parameters to parameter point
238  if(fAltModel->GetNuisanceParameters())
239  parameterPoint->add(*fAltModel->GetNuisanceParameters());
240 
241  delete allParams;
242 
243  // ***** ToyMCSampler specific *******
244 
245  // check whether TestStatSampler is a ToyMCSampler
246  ToyMCSampler *toymcs = dynamic_cast<ToyMCSampler*>(GetTestStatSampler());
247  if(toymcs) {
248  oocoutI((TObject*)0,InputArguments) << "Using a ToyMCSampler. Now configuring for Alt." << endl;
249 
250  // variable number of toys
251  if(fNToysAlt >= 0) toymcs->SetNToys(fNToysAlt);
252 
253  // set the global observables to be generated by the ToyMCSampler
254  toymcs->SetGlobalObservables(*fAltModel->GetGlobalObservables());
255 
256  // adaptive sampling
257  if(fNToysAltTail) {
258  oocoutI((TObject*)0,InputArguments) << "Adaptive Sampling" << endl;
259  if(GetTestStatSampler()->GetTestStatistic()->PValueIsRightTail()) {
260  toymcs->SetToysLeftTail(fNToysAltTail, obsTestStat);
261  }else{
262  toymcs->SetToysRightTail(fNToysAltTail, obsTestStat);
263  }
264  }else{
265  toymcs->SetToysBothTails(0, 0, obsTestStat); // disable adaptive sampling
266  }
267 
268  }
269 
270  return 0;
271 }
RooCmdArg Offset(Bool_t flag=kTRUE)
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:86
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
Class RooProfileLL implements the profile likelihood estimator for a given likelihood and set of para...
Definition: RooProfileLL.h:26
void SetToysLeftTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:216
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:175
RooCmdArg CloneData(Bool_t flag)
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
#define oocoutI(o, a)
Definition: RooMsgService.h:44
RooFit::MsgLevel globalKillBelow() const
MINIMIZER * minimizer()
Definition: RooProfileLL.h:38
static RooMsgService & instance()
Return reference to singleton instance.
STL namespace.
void SetToysRightTail(Double_t toys, Double_t threshold)
Definition: ToyMCSampler.h:221
void SetToysBothTails(Double_t toys, Double_t low_threshold, Double_t high_threshold)
Definition: ToyMCSampler.h:226
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
RooCmdArg GlobalObservables(const RooArgSet &globs)
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:146
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
void setGlobalKillBelow(RooFit::MsgLevel level)
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
const Bool_t kFALSE
Definition: RtypesCore.h:88
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:359
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
bool IsNLLOffset()
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
void RemoveConstantParameters(RooArgSet *set)
Definition: RooStatsUtils.h:62
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:461
Does a frequentist hypothesis test.
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.
RooCmdArg ConditionalObservables(const RooArgSet &set)
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31
RooCmdArg Constrain(const RooArgSet &params)