Logo ROOT   6.12/07
Reference Guide
HypoTestInverterOriginal.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id$
2 // Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
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::HypoTestInverterOriginal
12  \ingroup Roostats
13 
14 This class is now deprecated and to be replaced by the HypoTestInverter.
15 HypoTestInverterOriginal class for performing an hypothesis test inversion by scanning the hypothesis test results of the
16 HybridCalculator for various values of the parameter of interest. By looking at the confidence level curve of
17 the result an upper limit, where it intersects the desired confidence level, can be derived.
18 The class implements the RooStats::IntervalCalculator interface and returns an RooStats::HypoTestInverterResult class.
19 The result is a SimpleInterval, which via the method UpperLimit returns to the user the upper limit value.
20 
21 The HypoTestInverterOriginal implements various option for performing the scan. HypoTestInverterOriginal::RunFixedScan will scan using a fixed grid the parameter of interest. HypoTestInverterOriginal::RunAutoScan will perform an automatic scan to find optimally the curve and it will stop until the desired precision is obtained.
22 The confidence level value at a given point can be done via HypoTestInverterOriginal::RunOnePoint.
23 The class can scan the CLs+b values or alternatively CLs (if the method HypoTestInverterOriginal::UseCLs has been called).
24 
25 
26 New contributions to this class have been written by Matthias Wolf (advanced AutoRun algorithm)
27 **/
28 
29 // include header file of this class
31 
33 #include "RooStats/HybridResult.h"
36 
37 // include other header files
38 #include "RooAbsPdf.h"
39 #include "RooAbsData.h"
40 #include "RooRealVar.h"
41 #include "TMath.h"
42 
44 
45 using namespace RooStats;
46 using namespace std;
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 
51  fCalculator0(0),
52  fScannedVariable(0),
53  fResults(0),
54  fUseCLs(false),
55  fSize(0)
56 {
57  // default constructor (doesn't do anything)
58 }
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 
63  RooRealVar& scannedVariable, double size ) :
64  TNamed( ),
65  fCalculator0(&myhc0),
66  fScannedVariable(&scannedVariable),
67  fResults(0),
68  fUseCLs(false),
69  fSize(size)
70 {
71  // constructor from a reference to an HypoTestCalculator
72  // (it must be an HybridCalculator type) and a RooRealVar for the variable
73  SetName("HypoTestInverterOriginal");
74 
75 
77  if (hc == 0) {
78  Fatal("HypoTestInverterOriginal","Using non HybridCalculatorOriginal class IS NOT SUPPORTED");
79  }
80 
81 }
82 
83 ////////////////////////////////////////////////////////////////////////////////
84 
86 {
87  // destructor
88 
89  // delete the HypoTestInverterResult
90  if (fResults) delete fResults;
91 }
92 
93 ////////////////////////////////////////////////////////////////////////////////
94 
96  // create a new HypoTestInverterResult to hold all computed results
97  if (fResults == 0) {
98  TString results_name = this->GetName();
99  results_name += "_results";
101  fResults->SetTitle("HypoTestInverterOriginal Result");
102  }
104 }
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 
108 bool HypoTestInverterOriginal::RunAutoScan( double xMin, double xMax, double target, double epsilon, unsigned int numAlgorithm )
109 {
110  /// Search for the value of the parameter of interest (vary the
111  /// hypothesis being tested) in the specified range [xMin,xMax]
112  /// until the confidence level is compatible with the target value
113  /// within one time the estimated error (and the estimated error
114  /// should also become smaller than the specified parameter epsilon)
115 
116  // various sanity checks on the input parameters
117  if ( xMin>=xMax || xMin< fScannedVariable->getMin() || xMax>fScannedVariable->getMax() ) {
118  std::cout << "Error: problem with the specified range\n";
119  return false;
120  }
121  if ( target<=0 || target>=1 ) {
122  std::cout << "Error: problem with target value\n";
123  return false;
124  }
125  if ( epsilon>0.5-fabs(0.5-target) ) {
126  std::cout << "Error: problem with error value\n";
127  return false;
128  }
129  if ( numAlgorithm!=0 && numAlgorithm!=1 ) {
130  std::cout << "Error: invalid interpolation algorithm\n";
131  return false;
132  }
133 
134  CreateResults();
135 
136  // if ( TMath::AreEqualRel(target,1-Size()/2,DBL_EPSILON) ) { // to uncomment for ROOT 5.26
137  if ( fabs(1-target/(1-Size()/2))<DBL_EPSILON ) {
139  std::cout << "Target matches lower limit: de-activate interpolation in HypoTestInverterResult\n";
140  }
141  // if ( TMath::AreEqualRel(target,Size()/2,DBL_EPSILON) ) { // to uncomment for ROOT 5.26
142  if ( fabs(1-target/((Size()/2)))<DBL_EPSILON ) {
144  std::cout << "Target matches upper limit: de-activate interpolation in HypoTestInverterResult\n";
145  }
146 
147  // parameters of the algorithm that are hard-coded
148  const double nSigma = 1; // number of times the estimated error the final p-value should be from the target
149 
150  // backup some values to be restored at the end
151  const unsigned int nToys_backup = ((HybridCalculatorOriginal*)fCalculator0)->GetNumberOfToys();
152 
153  // check the 2 hypothesis tests specified as extrema in the constructor
154  double leftX = xMin;
155  if (!RunOnePoint(leftX)) return false;
156  double leftCL = fResults->GetYValue(fResults->ArraySize()-1);
157  double leftCLError = fResults->GetYError(fResults->ArraySize()-1);
158 
159  double rightX = xMax;
160  if (!RunOnePoint(rightX)) return false;
161  double rightCL = fResults->GetYValue(fResults->ArraySize()-1);
162  double rightCLError = fResults->GetYError(fResults->ArraySize()-1);
163 
164  if ( rightCL>target && leftCL>target ) {
165  std::cout << "The confidence level at both boundaries are both too large ( " << leftCL << " and " << rightCL << std::endl << "Run again with other boundaries or larger toy-MC statistics\n";
166  return false;
167  }
168  if ( rightCL<target && leftCL<target ) {
169  std::cout << "The confidence level at both boundaries are both too small ( " << leftCL << " and " << rightCL << std::endl << "Run again with other boundaries or larger toy-MC statistics\n";
170  return false;
171  }
172 
173  unsigned int nIteration = 2; // number of iteration performed by the algorithm
174  bool quitThisLoop = false; // flag to interrupt the search and quit cleanly
175 
176  double centerCL = 0;
177  double centerCLError = 0;
178 
179  // search for the value of the searched variable where the CL is
180  // within 1 sigma of the desired level and sigma smaller than
181  // epsilon.
182  do {
183  double x = 0;
184 
185  // safety checks
186  if (leftCL==rightCL) {
187  std::cout << "This cannot (and should not) happen... quit\n";
188  quitThisLoop = true;
189  } else if (leftX==rightX) {
190  std::cout << "This cannot (and should not) happen... quit\n";
191  quitThisLoop = true;
192  } else {
193 
194  // apply chosen type of interpolation algorithm
195  if (numAlgorithm==0) {
196  // exponential interpolation
197 
198  // add safety checks
199  if (!leftCL) leftCL = DBL_EPSILON;
200  if (!rightCL) rightCL = DBL_EPSILON;
201 
202  double a = (log(leftCL) - log(rightCL)) / (leftX - rightX);
203  double b = leftCL / exp(a * leftX);
204  x = (log(target) - log(b)) / a;
205 
206  // to do: do not allow next iteration outside the xMin,xMax interval
207  if (x<xMin || x>xMax || TMath::IsNaN(x)) {
208  std::cout << "Extrapolated value out of range or nan: exits\n";
209  quitThisLoop = true;
210  }
211  } else if (numAlgorithm==1) {
212  // linear interpolation
213 
214  double a = (leftCL-rightCL)/(leftX-rightX);
215  double b = leftCL-a*leftX;
216  x = (target-b)/a;
217 
218  if (x<xMin || x>xMax || TMath::IsNaN(x)) {
219  std::cout << "Extrapolated value out of range or nan: exits\n";
220  quitThisLoop = true;
221  }
222  } // end of interpolation algorithms
223  }
224 
225  if ( x==leftX || x==rightX ) {
226  std::cout << "Error: exit because interpolated value equals to a previous iteration\n";
227  quitThisLoop = true;
228  }
229 
230  // perform another hypothesis-test for value x
231  bool success = false;
232  if (!quitThisLoop) success = RunOnePoint(x);
233 
234  if (success) {
235 
236  nIteration++; // succeeded, increase the iteration counter
237  centerCL = fResults->GetYValue(fResults->ArraySize()-1);
238  centerCLError = fResults->GetYError(fResults->ArraySize()-1);
239 
240  // replace either the left or right point by this new point
241 
242  // test if the interval points are on different sides, then
243  // replace the one on the correct side with the center
244  if ( (leftCL > target) == (rightCL < target) ) {
245  if ( (centerCL > target) == (leftCL > target) ) {
246  leftX = x;
247  leftCL = centerCL;
248  leftCLError = centerCLError;
249  } else {
250  rightX = x;
251  rightCL = centerCL;
252  rightCLError = centerCLError;
253  }
254  // Otherwise replace the point farest away from target (measured in
255  // sigmas)
256  } else if ( (fabs(leftCL - target) / leftCLError) >
257  (fabs(rightCL - target) / rightCLError) ) {
258  leftX = x;
259  leftCL = centerCL;
260  leftCLError = centerCLError;
261  } else {
262  rightX = x;
263  rightCL = centerCL;
264  rightCLError = centerCLError;
265  }
266 
267  // if a point is found compatible with the target CL but with too
268  // large error, increase the number of toyMC
269  if ( fabs(centerCL-target) < nSigma*centerCLError && centerCLError > epsilon ) {
270  do {
271 
272  int nToys = ((HybridCalculatorOriginal*)fCalculator0)->GetNumberOfToys(); // current number of toys
273  int nToysTarget = (int) TMath::Max(nToys*1.5, 1.2*nToys*pow(centerCLError/epsilon,2)); // estimated number of toys until the target precision is reached
274 
275  std::cout << "Increasing the number of toys to: " << nToysTarget << std::endl;
276 
277  // run again the same point with more toyMC (run the complement number of toys)
278  ((HybridCalculatorOriginal*)fCalculator0)->SetNumberOfToys(nToysTarget-nToys);
279 
280  if (!RunOnePoint(x)) quitThisLoop=true;
281  nIteration++; // succeeded, increase the iteration counter
282  centerCL = fResults->GetYValue(fResults->ArraySize()-1);
283  centerCLError = fResults->GetYError(fResults->ArraySize()-1);
284 
285  // set the number of toys to reach the target
286  ((HybridCalculatorOriginal*)fCalculator0)->SetNumberOfToys(nToysTarget);
287  } while ( fabs(centerCL-target) < nSigma*centerCLError && centerCLError > epsilon && quitThisLoop==false ); // run this block again if it's still compatible with the target and the error still too large
288  }
289 
290  if (leftCL==rightCL) {
291  std::cout << "Algorithm failed: left and right CL are equal (no interpolation possible or more toy-MC statistics needed)\n";
292  quitThisLoop = true;
293  }
294  } // end running one more iteration
295 
296  } while ( ( fabs(centerCL-target) > nSigma*centerCLError || centerCLError > epsilon ) && quitThisLoop==false ); // end of the main 'do' loop
297 
298  // restore some parameters that might have been changed by the algorithm
299  ((HybridCalculatorOriginal*)fCalculator0)->SetNumberOfToys(nToys_backup);
300 
301  if ( quitThisLoop==true ) {
302  // abort and return 'false' to indicate fail status
303  std::cout << "Aborted the search because something happened\n";
304  return false;
305  }
306 
307  std::cout << "Converged in " << fResults->ArraySize() << " iterations\n";
308 
309  // finished: return 'true' for success status
310  return true;
311 }
312 
313 ////////////////////////////////////////////////////////////////////////////////
314 
315 bool HypoTestInverterOriginal::RunFixedScan( int nBins, double xMin, double xMax )
316 {
317  // Run a Fixed scan in npoints between min and max
318 
319  CreateResults();
320  // safety checks
321  if ( nBins<=0 ) {
322  std::cout << "Please provide nBins>0\n";
323  return false;
324  }
325  if ( nBins==1 && xMin!=xMax ) {
326  std::cout << "nBins==1 -> I will run for xMin (" << xMin << ")\n";
327  }
328  if ( xMin==xMax && nBins>1 ) {
329  std::cout << "xMin==xMax -> I will enforce nBins==1\n";
330  nBins = 1;
331  }
332  if ( xMin>xMax ) {
333  std::cout << "Please provide xMin (" << xMin << ") smaller that xMax (" << xMax << ")\n";
334  return false;
335  }
336 
337  for (int i=0; i<nBins; i++) {
338  double thisX = xMin+i*(xMax-xMin)/(nBins-1);
339  bool status = RunOnePoint(thisX);
340 
341  // check if failed status
342  if ( status==false ) {
343  std::cout << "Loop interrupted because of failed status\n";
344  return false;
345  }
346  }
347 
348  return true;
349 }
350 
351 ////////////////////////////////////////////////////////////////////////////////
352 /// run only one point
353 
355 {
356 
357  CreateResults();
358 
359  // check if thisX is in the range specified for fScannedVariable
360  if ( thisX<fScannedVariable->getMin() ) {
361  std::cout << "Out of range: using the lower bound on the scanned variable rather than " << thisX<< "\n";
362  thisX = fScannedVariable->getMin();
363  }
364  if ( thisX>fScannedVariable->getMax() ) {
365  std::cout << "Out of range: using the upper bound on the scanned variable rather than " << thisX<< "\n";
366  thisX = fScannedVariable->getMax();
367  }
368 
369  double oldValue = fScannedVariable->getVal();
370 
371  fScannedVariable->setVal(thisX);
372  std::cout << "Running for " << fScannedVariable->GetName() << " = " << thisX << endl;
373 
374  // compute the results
375  HypoTestResult* myHybridResult = fCalculator0->GetHypoTest();
376 
377  double lastXtested;
378  if ( fResults->ArraySize()!=0 ) lastXtested = fResults->GetXValue(fResults->ArraySize()-1);
379  else lastXtested = -999;
380 
381  if ( lastXtested==thisX ) {
382 
383  std::cout << "Merge with previous result\n";
384  HybridResult* latestResult = (HybridResult*) fResults->GetResult(fResults->ArraySize()-1);
385  latestResult->Add((HybridResult*)myHybridResult);
386  delete myHybridResult;
387 
388  } else {
389 
390  // fill the results in the HypoTestInverterResult array
391  fResults->fXValues.push_back(thisX);
392  fResults->fYObjects.Add(myHybridResult);
393  }
394 
395 
396  std::cout << "computed: " << fResults->GetYValue(fResults->ArraySize()-1) << endl;
397 
398  fScannedVariable->setVal(oldValue);
399 
400  return true;
401 }
virtual Double_t getMin(const char *name=0) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Double_t getMax(const char *name=0) const
bool RunAutoScan(double xMin, double xMax, double target, double epsilon=0.005, unsigned int numAlgorithm=0)
void Add(HybridResult *other)
add additional toy-MC experiments to the current results use the data test statistics of the added ob...
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:140
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
HypoTestResult is a base class for results from hypothesis tests.
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
HybridCalculatorOriginal class.
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
double GetYError(int index) const
function to return the estimated error on the value of the confidence level for the i^th entry in the...
Bool_t IsNaN(Double_t x)
Definition: TMath.h:777
STL namespace.
Double_t x[n]
Definition: legend1.C:17
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
double pow(double, double)
void UseCLs(bool on=true)
flag to switch between using CLsb (default) or CLs as confidence level
bool fInterpolateLowerLimit
two sided scan (look for lower/upper limit)
HypoTestCalculator is an interface class for a tools which produce RooStats HypoTestResults.
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
virtual Double_t Size() const
Get the size of the test (eg. rate of Type I error)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
auto * a
Definition: textangle.C:12
bool RunOnePoint(double thisX)
run only one point
REAL epsilon
Definition: triangle.c:617
virtual HypoTestResult * GetHypoTest() const =0
Namespace for the RooStats classes.
Definition: Asimov.h:20
This class is now deprecated and to be replaced by the HypoTestInverter.
#define ClassImp(name)
Definition: Rtypes.h:359
HypoTestInverterResult class holds the array of hypothesis test results and compute a confidence inte...
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results ...
virtual void Add(TObject *obj)
Definition: TList.h:87
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
bool RunFixedScan(int nBins, double xMin, double xMax)
double exp(double)
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
int ArraySize() const
number of entries in the results array
Class encapsulating the result of the HybridCalculatorOriginal.
Definition: HybridResult.h:25
double log(double)
std::vector< double > fXValues
number of points used to build expected p-values