Logo ROOT   6.10/09
Reference Guide
HypoTestResult.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Sven Kreiss
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  * Project: RooStats
13  * Package: RooFit/RooStats
14  * @(#)root/roofit/roostats:$Id$
15  * Authors:
16  * Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Sven Kreiss
17  *
18  *****************************************************************************/
19 
20 
21 /** \class RooStats::HypoTestResult
22  \ingroup Roostats
23 
24 HypoTestResult is a base class for results from hypothesis tests.
25 Any tool inheriting from HypoTestCalculator can return a HypoTestResult.
26 As such, it stores a p-value for the null-hypothesis (eg. background-only)
27 and an alternate hypothesis (eg. signal+background).
28 The p-values can also be transformed into confidence levels
29 (\f$CL_{b}\f$, \f$CL_{s+b}\f$) in a trivial way.
30 The ratio of the \f$CL_{s+b}\f$ to \f$CL_{b}\f$ is often called
31 \f$CL_{s}\f$, and is considered useful, though it is not a probability.
32 Finally, the p-value of the null can be transformed into a number of
33 equivalent Gaussian sigma using the Significance method.
34 
35 The p-value of the null for a given test statistic is rigorously defined and
36 this is the starting point for the following conventions.
37 
38 ### Conventions used in this class
39 
40 The p-value for the null and alternate are on the **same side** of the
41 observed value of the test statistic. This is the more standard
42 convention and avoids confusion when doing inverted tests.
43 
44 For exclusion, we also want the formula \f$CL_{s} = CL_{s+b} / CL_{b}\f$
45 to hold which therefore defines our conventions for \f$CL_{s+b}\f$ and
46 \f$CL_{b}\f$. \f$CL_{s}\f$ was specifically invented for exclusion
47 and therefore all quantities need be related through the assignments
48 as they are for exclusion: \f$CL_{s+b} = p_{s+b}\f$; \f$CL_{b} = p_{b}\f$. This
49 is derived by considering the scenarios of a powerful and not powerful
50 inverted test, where for the not so powerful test, \f$CL_{s}\f$ must be
51 close to one.
52 
53 For results of Hypothesis tests,
54 \f$CL_{s}\f$ has no similar direct interpretation as for exclusion and can
55 be larger than one.
56 
57 */
58 
61 #include "RooAbsReal.h"
62 
63 #include "RooStats/RooStatsUtils.h"
64 
65 #include <limits>
66 #define NaN numeric_limits<float>::quiet_NaN()
67 #define IsNaN(a) TMath::IsNaN(a)
68 
70 
71 using namespace RooStats;
72 using namespace std;
73 
74 ////////////////////////////////////////////////////////////////////////////////
75 /// Default constructor
76 
78  TNamed(name,name),
79  fNullPValue(NaN), fAlternatePValue(NaN),
80  fNullPValueError(0), fAlternatePValueError(0),
81  fTestStatisticData(NaN),
82  fAllTestStatisticsData(NULL),
83  fNullDistr(NULL), fAltDistr(NULL),
84  fNullDetailedOutput(NULL), fAltDetailedOutput(NULL), fFitInfo(NULL),
85  fPValueIsRightTail(kTRUE),
86  fBackgroundIsAlt(kFALSE)
87 {
88 }
89 
90 ////////////////////////////////////////////////////////////////////////////////
91 /// Alternate constructor
92 
94  TNamed(name,name),
95  fNullPValue(nullp), fAlternatePValue(altp),
103 {
104 }
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 /// copy constructor
108 
110  TNamed(other),
119 {
120  this->Append( &other );
121 }
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// Destructor
125 
127 {
128  if( fNullDistr ) delete fNullDistr;
129  if( fAltDistr ) delete fAltDistr;
130 
133 
135 }
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// assignment operator
139 
141  if (this == &other) return *this;
142  SetName(other.GetName());
143  SetTitle(other.GetTitle());
144  fNullPValue = other.fNullPValue;
149 
152  if( fNullDistr ) { delete fNullDistr; fNullDistr = NULL; }
153  if( fAltDistr ) { delete fAltDistr; fAltDistr = NULL; }
156  if (fFitInfo) { delete fFitInfo; fFitInfo = NULL; }
157 
160 
161  this->Append( &other );
162 
163  return *this;
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 /// Add additional toy-MC experiments to the current results.
168 /// Use the data test statistics of the added object if it is not already
169 /// set (otherwise, ignore the new one).
170 
172  if(fNullDistr)
174  else
176 
177  if(fAltDistr)
178  fAltDistr->Add(other->GetAltDistribution());
179  else
181 
182 
183  if( fNullDetailedOutput ) {
185  }else{
187  }
188 
189  if( fAltDetailedOutput ) {
191  }else{
193  }
194 
195  if( fFitInfo ) {
196  if( other->GetFitInfo() ) fFitInfo->append( *other->GetFitInfo() );
197  }else{
198  if( other->GetFitInfo() ) fFitInfo = new RooDataSet( *other->GetFitInfo() );
199  }
200 
201  // if no data is present use the other HypoTestResult's data
203 
206 }
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 
211  fAltDistr = alt;
213 }
214 
215 ////////////////////////////////////////////////////////////////////////////////
216 
218  fNullDistr = null;
220 }
221 
222 ////////////////////////////////////////////////////////////////////////////////
223 
225  fTestStatisticData = tsd;
226 
229 }
230 
231 ////////////////////////////////////////////////////////////////////////////////
232 
235  delete fAllTestStatisticsData;
237  }
238  if (tsd) fAllTestStatisticsData = (const RooArgList*)tsd->snapshot();
239 
242  if( firstTS ) SetTestStatisticData( firstTS->getVal() );
243  }
244 }
245 
246 ////////////////////////////////////////////////////////////////////////////////
247 
249  fPValueIsRightTail = pr;
250 
253 }
254 
255 ////////////////////////////////////////////////////////////////////////////////
256 
258  return !IsNaN(fTestStatisticData);
259 }
260 
261 ////////////////////////////////////////////////////////////////////////////////
262 
264  // compute error on Null pvalue
265  return fNullPValueError;
266 }
267 
268 ////////////////////////////////////////////////////////////////////////////////
269 /// compute \f$CL_{b}\f$ error
270 /// \f$CL_{b}\f$ = 1 - NullPValue()
271 /// must use opposite condition that routine above
272 
275 }
276 
277 ////////////////////////////////////////////////////////////////////////////////
278 
281 }
282 
283 ////////////////////////////////////////////////////////////////////////////////
284 /// Taylor expansion series approximation for standard deviation (error propagation)
285 
288 }
289 
290 ////////////////////////////////////////////////////////////////////////////////
291 /// Returns an estimate of the error on \f$CL_{s}\f$ through combination of the
292 /// errors on \f$CL_{b}\f$ and \f$CL_{s+b}\f$:
293 /// \f[
294 /// \sigma_{CL_s} = CL_s
295 /// \sqrt{\left( \frac{\sigma_{CL_{s+b}}}{CL_{s+b}} \right)^2 + \left( \frac{\sigma_{CL_{b}}}{CL_{b}} \right)^2}
296 /// \f]
297 
299  if(!fAltDistr || !fNullDistr) return 0.0;
300 
301  // unsigned const int n_b = fNullDistr->GetSamplingDistribution().size();
302  // unsigned const int n_sb = fAltDistr->GetSamplingDistribution().size();
303 
304  // if CLb() == 0 CLs = -1 so return a -1 error
305  if (CLb() == 0 ) return -1;
306 
307  double cl_b_err2 = pow(CLbError(),2);
308  double cl_sb_err2 = pow(CLsplusbError(),2);
309 
310  return TMath::Sqrt(cl_sb_err2 + cl_b_err2 * pow(CLs(),2))/CLb();
311 }
312 
313 ////////////////////////////////////////////////////////////////////////////////
314 /// updates the pvalue if sufficient data is available
315 
316 void HypoTestResult::UpdatePValue(const SamplingDistribution* distr, Double_t &pvalue, Double_t &perror, Bool_t /*isNull*/) {
317  if(IsNaN(fTestStatisticData)) return;
318  if(!distr) return;
319 
320  /* Got to be careful for discrete distributions:
321  * To get the right behaviour for limits, the p-value must
322  * include the value of fTestStatistic both for Alt and Null cases
323  */
324  if(fPValueIsRightTail) {
325  pvalue = distr->IntegralAndError(perror, fTestStatisticData, RooNumber::infinity(), kTRUE,
326  kTRUE , kTRUE ); // always closed interval [ fTestStatistic, inf ]
327 
328  }else{
329  pvalue = distr->IntegralAndError(perror, -RooNumber::infinity(), fTestStatisticData, kTRUE,
330  kTRUE, kTRUE ); // // always closed [ -inf, fTestStatistic ]
331  }
332 }
333 
334 ////////////////////////////////////////////////////////////////////////////////
335 /// Print out some information about the results
336 /// Note: use Alt/Null labels for the hypotheses here as the Null
337 /// might be the s+b hypothesis.
338 
340 {
341  bool fromToys = (fAltDistr || fNullDistr);
342 
343  std::cout << std::endl << "Results " << GetName() << ": " << endl;
344  std::cout << " - Null p-value = " << NullPValue();
345  if (fromToys) std::cout << " +/- " << NullPValueError();
346  std::cout << std::endl;
347  std::cout << " - Significance = " << Significance();
348  if (fromToys) std::cout << " +/- " << SignificanceError() << " sigma";
349  std::cout << std::endl;
350  if(fAltDistr)
351  std::cout << " - Number of Alt toys: " << fAltDistr->GetSize() << std::endl;
352  if(fNullDistr)
353  std::cout << " - Number of Null toys: " << fNullDistr->GetSize() << std::endl;
354 
355  if (HasTestStatisticData() ) std::cout << " - Test statistic evaluated on data: " << fTestStatisticData << std::endl;
356  std::cout << " - CL_b: " << CLb();
357  if (fromToys) std::cout << " +/- " << CLbError();
358  std::cout << std::endl;
359  std::cout << " - CL_s+b: " << CLsplusb();
360  if (fromToys) std::cout << " +/- " << CLsplusbError();
361  std::cout << std::endl;
362  std::cout << " - CL_s: " << CLs();
363  if (fromToys) std::cout << " +/- " << CLsError();
364  std::cout << std::endl;
365 
366  return;
367 }
Bool_t HasTestStatisticData(void) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
void SetAltDistribution(SamplingDistribution *alt)
virtual Double_t CLb() const
Convert NullPValue into a "confidence level".
const char Option_t
Definition: RtypesCore.h:62
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
RooDataSet * GetAltDetailedOutput(void) const
const RooArgList * fAllTestStatisticsData
RooDataSet * GetFitInfo(void) const
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
HypoTestResult is a base class for results from hypothesis tests.
Double_t NullPValueError() const
The error on the Null p-value.
HypoTestResult & operator=(const HypoTestResult &other)
assignment operator
RooDataSet * fNullDetailedOutput
void UpdatePValue(const SamplingDistribution *distr, Double_t &pvalue, Double_t &perror, Bool_t pIsRightTail)
updates the pvalue if sufficient data is available
bool Bool_t
Definition: RtypesCore.h:59
SamplingDistribution * GetAltDistribution(void) const
STL namespace.
#define NULL
Definition: RtypesCore.h:88
Double_t CLbError() const
The error on the "confidence level" of the null hypothesis.
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
virtual Double_t Significance() const
familiar name for the Null p-value in terms of 1-sided Gaussian significance
null_t< F > null()
Double_t CLsError() const
The error on the ratio .
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Bool_t GetPValueIsRightTail(void) const
double pow(double, double)
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
void SetAllTestStatisticsData(const RooArgList *tsd)
Int_t getSize() const
virtual ~HypoTestResult()
destructor
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 Double_t NullPValue() const
Return p-value for null hypothesis.
Double_t GetTestStatisticData(void) const
SamplingDistribution * fAltDistr
#define IsNaN(a)
static Double_t infinity()
Return internal infinity representation.
Definition: RooNumber.cxx:49
Int_t GetSize() const
size of samples
SamplingDistribution * GetNullDistribution(void) const
void Add(const SamplingDistribution *other)
merge two sampling distributions
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.
HypoTestResult(const char *name=0)
default constructor
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
virtual Double_t CLs() const
is simply (not a method, but a quantity)
const Bool_t kFALSE
Definition: RtypesCore.h:92
void SetPValueIsRightTail(Bool_t pr)
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
RooDataSet * GetNullDetailedOutput(void) const
Namespace for the RooStats classes.
Definition: Asimov.h:20
SamplingDistribution * fNullDistr
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
void append(RooDataSet &data)
Add all data points of given data set to this data set.
Double_t CLsplusbError() const
The error on the "confidence level" of the alternative hypothesis.
Double_t IntegralAndError(Double_t &error, Double_t low, Double_t high, Bool_t normalize=kTRUE, Bool_t lowClosed=kTRUE, Bool_t highClosed=kFALSE) const
numerical integral in these limits including error estimation
Double_t SignificanceError() const
The error on the significance, computed from NullPValueError via error propagation.
RooDataSet * fAltDetailedOutput
Bool_t GetBackGroundIsAlt(void) const
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
#define NaN
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
const Bool_t kTRUE
Definition: RtypesCore.h:91
void SetNullDistribution(SamplingDistribution *null)
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48