Logo ROOT   6.08/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 /**
12 Same purpose as HybridCalculatorOriginal, but different implementation.
13 
14 This is the "generic" version that works with any TestStatSampler. The
15 HybridCalculator derives from this class but explicitly uses the
16 ToyMCSampler as its TestStatSampler.
17 */
18 
20 #include "RooStats/ToyMCSampler.h"
24 
25 #include "RooAddPdf.h"
26 
27 #include "RooRandom.h"
28 
29 
31 
32 using namespace RooStats;
33 using namespace std;
34 
35 ////////////////////////////////////////////////////////////////////////////////
36 /// Constructor. When test stat sampler is not provided
37 /// uses ToyMCSampler and RatioOfProfiledLikelihoodsTestStat
38 /// and nToys = 1000.
39 /// User can : GetTestStatSampler()->SetNToys( # )
40 
42  const RooAbsData &data,
43  const ModelConfig &altModel,
44  const ModelConfig &nullModel,
45  TestStatSampler *sampler
46  ) :
47  fAltModel(&altModel),
48  fNullModel(&nullModel),
49  fData(&data),
50  fTestStatSampler(sampler),
51  fDefaultSampler(0),
52  fDefaultTestStat(0),
53  fAltToysSeed(0)
54 {
55  if(!sampler){
57  = new RatioOfProfiledLikelihoodsTestStat(*nullModel.GetPdf(),
58  *altModel.GetPdf(),
59  altModel.GetSnapshot());
60 
63  }
64 
65 
66 }
67 
68 ////////////////////////////////////////////////////////////////////////////////
69 /// common setup for both models
70 
75 
76  // for this model
77  model.LoadSnapshot();
79  fTestStatSampler->SetPdf(*model.GetPdf());
81  // global observables or nuisanance pdf will be set by the derived classes
82  // (e.g. Frequentist or HybridCalculator)
83 }
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 
95  // several possibilities:
96  // no prior nuisance given and no nuisance parameters: ok
97  // no prior nuisance given but nuisance parameters: error
98  // prior nuisance given for some nuisance parameters:
99  // - nuisance parameters are constant, so they don't float in test statistic
100  // - nuisance parameters are floating, so they do float in test statistic
101 
102  // initial setup
103  PreHook();
104  const_cast<ModelConfig*>(fNullModel)->GuessObsAndNuisance(*fData);
105  const_cast<ModelConfig*>(fAltModel)->GuessObsAndNuisance(*fData);
106 
107  const RooArgSet * nullSnapshot = fNullModel->GetSnapshot();
108  if(nullSnapshot == NULL) {
109  oocoutE((TObject*)0,Generation) << "Null model needs a snapshot. Set using modelconfig->SetSnapshot(poi)." << endl;
110  return 0;
111  }
112 
113  // CheckHook
114  if(CheckHook() != 0) {
115  oocoutE((TObject*)0,Generation) << "There was an error in CheckHook(). Stop." << endl;
116  return 0;
117  }
118 
120  oocoutE((TObject*)0,InputArguments) << "Test Statistic Sampler or Test Statistics not defined. Stop." << endl;
121  return 0;
122  }
123 
124  // get a big list of all variables for convenient switching
125  RooArgSet *nullParams = fNullModel->GetPdf()->getParameters(*fData);
126  RooArgSet *altParams = fAltModel->GetPdf()->getParameters(*fData);
127  // save all parameters so we can set them back to what they were
128  RooArgSet *bothParams = fNullModel->GetPdf()->getParameters(*fData);
129  bothParams->add(*altParams,false);
130  RooArgSet *saveAll = (RooArgSet*) bothParams->snapshot();
131 
132  // check whether we have a ToyMCSampler and if so, keep a pointer to it
133  ToyMCSampler* toymcs = dynamic_cast<ToyMCSampler*>( fTestStatSampler );
134 
135 
136  // evaluate test statistic on data
137  RooArgSet nullP(*nullSnapshot);
138  double obsTestStat;
139 
140  RooArgList* allTS = NULL;
141  if( toymcs ) {
142  allTS = toymcs->EvaluateAllTestStatistics(*const_cast<RooAbsData*>(fData), nullP);
143  if (!allTS) return 0;
144  //oocoutP((TObject*)0,Generation) << "All Test Statistics on data: " << endl;
145  //allTS->Print("v");
146  RooRealVar* firstTS = (RooRealVar*)allTS->at(0);
147  obsTestStat = firstTS->getVal();
148  if (allTS->getSize()<=1) {
149  delete allTS;
150  allTS= 0; // don't save
151  }
152  }else{
153  obsTestStat = fTestStatSampler->EvaluateTestStatistic(*const_cast<RooAbsData*>(fData), nullP);
154  }
155  oocoutP((TObject*)0,Generation) << "Test Statistic on data: " << obsTestStat << endl;
156 
157  // set parameters back ... in case the evaluation of the test statistic
158  // modified something (e.g. a nuisance parameter that is not randomized
159  // must be set here)
160  *bothParams = *saveAll;
161 
162 
163 
164  // Generate sampling distribution for null
166  RooArgSet paramPointNull(*fNullModel->GetParametersOfInterest());
167  if(PreNullHook(&paramPointNull, obsTestStat) != 0) {
168  oocoutE((TObject*)0,Generation) << "PreNullHook did not return 0." << endl;
169  }
170  SamplingDistribution* samp_null = NULL;
171  RooDataSet* detOut_null = NULL;
172  if(toymcs) {
173  detOut_null = toymcs->GetSamplingDistributions(paramPointNull);
174  if( detOut_null ) {
175  samp_null = new SamplingDistribution( detOut_null->GetName(), detOut_null->GetTitle(), *detOut_null );
176  if (detOut_null->get()->getSize()<=1) {
177  delete detOut_null;
178  detOut_null= 0;
179  }
180  }
181  }else samp_null = fTestStatSampler->GetSamplingDistribution(paramPointNull);
182 
183  // set parameters back
184  *bothParams = *saveAll;
185 
186  // Generate sampling distribution for alternate
188  RooArgSet paramPointAlt(*fAltModel->GetParametersOfInterest());
189  if(PreAltHook(&paramPointAlt, obsTestStat) != 0) {
190  oocoutE((TObject*)0,Generation) << "PreAltHook did not return 0." << endl;
191  }
192  SamplingDistribution* samp_alt = NULL;
193  RooDataSet* detOut_alt = NULL;
194  if(toymcs) {
195 
196  // case of re-using same toys for every points
197  // set a given seed
198  unsigned int prevSeed = 0;
199  if (fAltToysSeed > 0) {
200  prevSeed = RooRandom::integer(std::numeric_limits<unsigned int>::max()-1)+1; // want to avoid zero value
202  }
203 
204  detOut_alt = toymcs->GetSamplingDistributions(paramPointAlt);
205  if( detOut_alt ) {
206  samp_alt = new SamplingDistribution( detOut_alt->GetName(), detOut_alt->GetTitle(), *detOut_alt );
207  if (detOut_alt->get()->getSize()<=1) {
208  delete detOut_alt;
209  detOut_alt= 0;
210  }
211  }
212 
213  // restore the seed
214  if (prevSeed > 0) {
215  RooRandom::randomGenerator()->SetSeed(prevSeed);
216  }
217 
218  }else samp_alt = fTestStatSampler->GetSamplingDistribution(paramPointAlt);
219 
220 
221  // create result
222  string resultname = "HypoTestCalculator_result";
223  HypoTestResult* res = new HypoTestResult(resultname.c_str());
225  res->SetTestStatisticData(obsTestStat);
226  res->SetAltDistribution(samp_alt);
227  res->SetNullDistribution(samp_null);
228  res->SetAltDetailedOutput( detOut_alt );
229  res->SetNullDetailedOutput( detOut_null );
230  res->SetAllTestStatisticsData( allTS );
231 
232  const RooArgSet *aset = GetFitInfo();
233  if (aset != NULL) {
234  RooDataSet *dset = new RooDataSet(NULL, NULL, *aset);
235  dset->add(*aset);
236  res->SetFitInfo( dset );
237  }
238 
239  *bothParams = *saveAll;
240  delete allTS;
241  delete bothParams;
242  delete saveAll;
243  delete altParams;
244  delete nullParams;
245  delete nullSnapshot;
246  PostHook();
247  return res;
248 }
249 
250 ////////////////////////////////////////////////////////////////////////////////
251 /// to re-use same toys for alternate hypothesis
252 
254  fAltToysSeed = RooRandom::integer(std::numeric_limits<unsigned int>::max()-1)+1;
255 }
256 
257 
258 
void SetupSampler(const ModelConfig &model) const
common setup for both models
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
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:52
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:259
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
virtual TestStatistic * GetTestStatistic() const =0
virtual RooDataSet * GetSamplingDistributions(RooArgSet &paramPoint)
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:46
virtual void SetObservables(const RooArgSet &)=0
#define oocoutE(o, a)
Definition: RooMsgService.h:48
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
virtual RooArgList * EvaluateAllTestStatistics(RooAbsData &data, const RooArgSet &poi)
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:37
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:47
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:99
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:244
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:247
#define ClassImp(name)
Definition: Rtypes.h:279
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:560
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:250
#define NULL
Definition: Rtypes.h:82
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:52