Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
11 #include "RooStats/FrequentistCalculator.h"
12 #include "RooStats/ToyMCSampler.h"
13 #include "RooStats/RooStatsUtils.h"
14 #include "RooStats/DetailedOutputAggregator.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 
31 ClassImp(RooStats::FrequentistCalculator);
32 
33 using namespace RooStats;
34 using namespace std;
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 
38 void FrequentistCalculator::PreHook() const {
39  if (fFitInfo != NULL) {
40  delete fFitInfo;
41  fFitInfo = NULL;
42  }
43  if (fStoreFitInfo) {
44  fFitInfo = new RooArgSet();
45  }
46 }
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 
50 void FrequentistCalculator::PostHook() const {
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;
89  RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
90  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
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;
116  RooMsgService::instance().setGlobalKillBelow(msglevel);
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;
198  RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
199  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
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;
226  RooMsgService::instance().setGlobalKillBelow(msglevel);
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 }