Logo ROOT   master
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 }
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:88
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
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
#define oocoutI(o, a)
Definition: RooMsgService.h:45
RooFit::MsgLevel globalKillBelow() const
MINIMIZER * minimizer()
Definition: RooProfileLL.h:38
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
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.
STL namespace.
RooCmdArg GlobalObservables(const RooArgSet &globs)
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 Constrain(const RooArgSet &params)
virtual void SetNToys(const Int_t ntoy)
Definition: ToyMCSampler.h:146
int PreNullHook(RooArgSet *parameterPoint, double obsTestStat) const
configure TestStatSampler for the Null run
RooCmdArg ConditionalObservables(const RooArgSet &set)
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
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
RooCmdArg CloneData(Bool_t flag)
Namespace for the RooStats classes.
Definition: Asimov.h:19
#define ClassImp(name)
Definition: Rtypes.h:365
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
RooStatsConfig & GetGlobalRooStatsConfig()
Retrieve the config object which can be used to set flags for things like offsetting the likelihood o...
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:69
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:517
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.
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31