Logo ROOT   6.12/07
Reference Guide
HypoTestCalculatorGeneric.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
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 /** \class RooStats::HypoTestCalculatorGeneric
12  \ingroup Roostats
13 
14 Common base class for the Hypothesis Test Calculators.
15 It is not designed to use directly but via its derived classes
16 
17 Same purpose as HybridCalculatorOriginal, but different implementation.
18 
19 This is the "generic" version that works with any TestStatSampler. The
20 HybridCalculator derives from this class but explicitly uses the
21 ToyMCSampler as its TestStatSampler.
22 
23 */
24 
26 #include "RooStats/ToyMCSampler.h"
30 
31 #include "RooAddPdf.h"
32 
33 #include "RooRandom.h"
34 
35 
37 
38 using namespace RooStats;
39 using namespace std;
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 /// Constructor. When test stat sampler is not provided
43 /// uses ToyMCSampler and RatioOfProfiledLikelihoodsTestStat
44 /// and nToys = 1000.
45 /// User can : GetTestStatSampler()->SetNToys( # )
46 
48  const RooAbsData &data,
49  const ModelConfig &altModel,
50  const ModelConfig &nullModel,
51  TestStatSampler *sampler
52  ) :
53  fAltModel(&altModel),
54  fNullModel(&nullModel),
55  fData(&data),
56  fTestStatSampler(sampler),
57  fDefaultSampler(0),
58  fDefaultTestStat(0),
59  fAltToysSeed(0)
60 {
61  if(!sampler){
63  = new RatioOfProfiledLikelihoodsTestStat(*nullModel.GetPdf(),
64  *altModel.GetPdf(),
65  altModel.GetSnapshot());
66 
69  }
70 
71 
72 }
73 
74 ////////////////////////////////////////////////////////////////////////////////
75 /// common setup for both models
76 
81 
82  // for this model
83  model.LoadSnapshot();
85  fTestStatSampler->SetPdf(*model.GetPdf());
87  // global observables or nuisance pdf will be set by the derived classes
88  // (e.g. Frequentist or HybridCalculator)
89 }
90 
91 ////////////////////////////////////////////////////////////////////////////////
92 
96 }
97 
98 ////////////////////////////////////////////////////////////////////////////////
99 /// several possibilities:
100 /// no prior nuisance given and no nuisance parameters: ok
101 /// no prior nuisance given but nuisance parameters: error
102 /// prior nuisance given for some nuisance parameters:
103 /// - nuisance parameters are constant, so they don't float in test statistic
104 /// - nuisance parameters are floating, so they do float in test statistic
105 
107 
108  // initial setup
109  PreHook();
110  const_cast<ModelConfig*>(fNullModel)->GuessObsAndNuisance(*fData);
111  const_cast<ModelConfig*>(fAltModel)->GuessObsAndNuisance(*fData);
112 
113  const RooArgSet * nullSnapshot = fNullModel->GetSnapshot();
114  if(nullSnapshot == NULL) {
115  oocoutE((TObject*)0,Generation) << "Null model needs a snapshot. Set using modelconfig->SetSnapshot(poi)." << endl;
116  return 0;
117  }
118 
119  // CheckHook
120  if(CheckHook() != 0) {
121  oocoutE((TObject*)0,Generation) << "There was an error in CheckHook(). Stop." << endl;
122  return 0;
123  }
124 
126  oocoutE((TObject*)0,InputArguments) << "Test Statistic Sampler or Test Statistics not defined. Stop." << endl;
127  return 0;
128  }
129 
130  // get a big list of all variables for convenient switching
131  RooArgSet *nullParams = fNullModel->GetPdf()->getParameters(*fData);
132  RooArgSet *altParams = fAltModel->GetPdf()->getParameters(*fData);
133  // save all parameters so we can set them back to what they were
134  RooArgSet *bothParams = fNullModel->GetPdf()->getParameters(*fData);
135  bothParams->add(*altParams,false);
136  RooArgSet *saveAll = (RooArgSet*) bothParams->snapshot();
137 
138  // check whether we have a ToyMCSampler and if so, keep a pointer to it
139  ToyMCSampler* toymcs = dynamic_cast<ToyMCSampler*>( fTestStatSampler );
140 
141 
142  // evaluate test statistic on data
143  RooArgSet nullP(*nullSnapshot);
144  double obsTestStat;
145 
146  RooArgList* allTS = NULL;
147  if( toymcs ) {
148  allTS = toymcs->EvaluateAllTestStatistics(*const_cast<RooAbsData*>(fData), nullP);
149  if (!allTS) return 0;
150  //oocoutP((TObject*)0,Generation) << "All Test Statistics on data: " << endl;
151  //allTS->Print("v");
152  RooRealVar* firstTS = (RooRealVar*)allTS->at(0);
153  obsTestStat = firstTS->getVal();
154  if (allTS->getSize()<=1) {
155  delete allTS;
156  allTS= 0; // don't save
157  }
158  }else{
159  obsTestStat = fTestStatSampler->EvaluateTestStatistic(*const_cast<RooAbsData*>(fData), nullP);
160  }
161  oocoutP((TObject*)0,Generation) << "Test Statistic on data: " << obsTestStat << endl;
162 
163  // set parameters back ... in case the evaluation of the test statistic
164  // modified something (e.g. a nuisance parameter that is not randomized
165  // must be set here)
166  *bothParams = *saveAll;
167 
168 
169 
170  // Generate sampling distribution for null
172  RooArgSet paramPointNull(*fNullModel->GetParametersOfInterest());
173  if(PreNullHook(&paramPointNull, obsTestStat) != 0) {
174  oocoutE((TObject*)0,Generation) << "PreNullHook did not return 0." << endl;
175  }
176  SamplingDistribution* samp_null = NULL;
177  RooDataSet* detOut_null = NULL;
178  if(toymcs) {
179  detOut_null = toymcs->GetSamplingDistributions(paramPointNull);
180  if( detOut_null ) {
181  samp_null = new SamplingDistribution( detOut_null->GetName(), detOut_null->GetTitle(), *detOut_null );
182  if (detOut_null->get()->getSize()<=1) {
183  delete detOut_null;
184  detOut_null= 0;
185  }
186  }
187  }else samp_null = fTestStatSampler->GetSamplingDistribution(paramPointNull);
188 
189  // set parameters back
190  *bothParams = *saveAll;
191 
192  // Generate sampling distribution for alternate
194  RooArgSet paramPointAlt(*fAltModel->GetParametersOfInterest());
195  if(PreAltHook(&paramPointAlt, obsTestStat) != 0) {
196  oocoutE((TObject*)0,Generation) << "PreAltHook did not return 0." << endl;
197  }
198  SamplingDistribution* samp_alt = NULL;
199  RooDataSet* detOut_alt = NULL;
200  if(toymcs) {
201 
202  // case of re-using same toys for every points
203  // set a given seed
204  unsigned int prevSeed = 0;
205  if (fAltToysSeed > 0) {
206  prevSeed = RooRandom::integer(std::numeric_limits<unsigned int>::max()-1)+1; // want to avoid zero value
208  }
209 
210  detOut_alt = toymcs->GetSamplingDistributions(paramPointAlt);
211  if( detOut_alt ) {
212  samp_alt = new SamplingDistribution( detOut_alt->GetName(), detOut_alt->GetTitle(), *detOut_alt );
213  if (detOut_alt->get()->getSize()<=1) {
214  delete detOut_alt;
215  detOut_alt= 0;
216  }
217  }
218 
219  // restore the seed
220  if (prevSeed > 0) {
221  RooRandom::randomGenerator()->SetSeed(prevSeed);
222  }
223 
224  }else samp_alt = fTestStatSampler->GetSamplingDistribution(paramPointAlt);
225 
226 
227  // create result
228  string resultname = "HypoTestCalculator_result";
229  HypoTestResult* res = new HypoTestResult(resultname.c_str());
231  res->SetTestStatisticData(obsTestStat);
232  res->SetAltDistribution(samp_alt);
233  res->SetNullDistribution(samp_null);
234  res->SetAltDetailedOutput( detOut_alt );
235  res->SetNullDetailedOutput( detOut_null );
236  res->SetAllTestStatisticsData( allTS );
237 
238  const RooArgSet *aset = GetFitInfo();
239  if (aset != NULL) {
240  RooDataSet *dset = new RooDataSet(NULL, NULL, *aset);
241  dset->add(*aset);
242  res->SetFitInfo( dset );
243  }
244 
245  *bothParams = *saveAll;
246  delete allTS;
247  delete bothParams;
248  delete saveAll;
249  delete altParams;
250  delete nullParams;
251  delete nullSnapshot;
252  PostHook();
253  return res;
254 }
255 
256 ////////////////////////////////////////////////////////////////////////////////
257 /// to re-use same toys for alternate hypothesis
258 
260  fAltToysSeed = RooRandom::integer(std::numeric_limits<unsigned int>::max()-1)+1;
261 }
void SetupSampler(const ModelConfig &model) const
common setup for both models
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
void UseSameAltToys()
to re-use same toys for alternate hypothesis
void SetAltDistribution(SamplingDistribution *alt)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:237
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 SetParametersForTestStat(const RooArgSet &)=0
void SetFitInfo(RooDataSet *d)
void LoadSnapshot() const
load the snapshot from ws if it exists
virtual TestStatistic * GetTestStatistic() const =0
virtual RooDataSet * GetSamplingDistributions(RooArgSet &paramPoint)
Use for serial and parallel runs.
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
TestStatSampler is an interface class for a tools which produce RooStats SamplingDistributions.
HypoTestResult is a base class for results from hypothesis tests.
virtual void SetPdf(RooAbsPdf &)=0
static UInt_t integer(UInt_t max, TRandom *generator=randomGenerator())
Return an integer uniformly distributed from [0,n-1].
Definition: RooRandom.cxx:102
virtual void SetSamplingDistName(const char *name)=0
STL namespace.
virtual void SetNuisanceParameters(const RooArgSet &)=0
Common base class for the Hypothesis Test Calculators.
#define oocoutP(o, a)
Definition: RooMsgService.h:45
virtual void SetObservables(const RooArgSet &)=0
#define oocoutE(o, a)
Definition: RooMsgService.h:47
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:589
virtual RooArgList * EvaluateAllTestStatistics(RooAbsData &data, const RooArgSet &poi)
Evaluate all test statistics, returning result and any detailed output.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
virtual SamplingDistribution * GetSamplingDistribution(RooArgSet &paramsOfInterest)=0
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
void SetAllTestStatisticsData(const RooArgList *tsd)
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
Int_t getSize() const
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
virtual bool PValueIsRightTail(void) const
Defines the sign convention of the test statistic. Overwrite function if necessary.
Definition: TestStatistic.h:45
void SetAltDetailedOutput(RooDataSet *d)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
void SetTestStatisticData(const Double_t tsd)
This class simply holds a sampling distribution of some test statistic.
TestStatistic that returns the ratio of profiled likelihoods.
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event &#39;index&#39;.
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0)
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
void SetPValueIsRightTail(Bool_t pr)
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:222
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:225
#define ClassImp(name)
Definition: Rtypes.h:359
virtual int PreNullHook(RooArgSet *, double) const
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:537
virtual const RooArgSet * GetFitInfo() const
HypoTestCalculatorGeneric(const RooAbsData &data, const ModelConfig &altModel, const ModelConfig &nullModel, TestStatSampler *sampler=0)
Constructor.
Mother of all ROOT objects.
Definition: TObject.h:37
void SetNullDetailedOutput(RooDataSet *d)
virtual Double_t EvaluateTestStatistic(RooAbsData &data, RooArgSet &paramsOfInterest)=0
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:228
virtual int PreAltHook(RooArgSet *, double) const
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing) ...
void SetNullDistribution(SamplingDistribution *null)
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48