Logo ROOT  
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
14This class is now deprecated and to be replaced by the HypoTestInverter.
15HypoTestInverterOriginal class for performing an hypothesis test inversion by scanning the hypothesis test results of the
16HybridCalculator for various values of the parameter of interest. By looking at the confidence level curve of
17the result an upper limit, where it intersects the desired confidence level, can be derived.
18The class implements the RooStats::IntervalCalculator interface and returns an RooStats::HypoTestInverterResult class.
19The result is a SimpleInterval, which via the method UpperLimit returns to the user the upper limit value.
20
21The 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.
22The confidence level value at a given point can be done via HypoTestInverterOriginal::RunOnePoint.
23The class can scan the CLs+b values or alternatively CLs (if the method HypoTestInverterOriginal::UseCLs has been called).
24
25
26New contributions to this class have been written by Matthias Wolf (advanced AutoRun algorithm)
27**/
28
29// include header file of this class
31
36
37// include other header files
38#include "RooAbsPdf.h"
39#include "RooAbsData.h"
40#include "RooRealVar.h"
41#include "TMath.h"
42
44
45using namespace RooStats;
46using namespace std;
47
48////////////////////////////////////////////////////////////////////////////////
49
50HypoTestInverterOriginal::HypoTestInverterOriginal( ) :
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
108bool 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
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
315bool HypoTestInverterOriginal::RunFixedScan( int nBins, double xMin, double xMax )
316{
317 // Run a Fixed scan in npoints between min and max
318
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
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";
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}
#define b(i)
Definition: RSha256.hxx:100
#define ClassImp(name)
Definition: Rtypes.h:365
double pow(double, double)
double exp(double)
double log(double)
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:278
HybridCalculatorOriginal class.
Class encapsulating the result of the HybridCalculatorOriginal.
Definition: HybridResult.h:25
void Add(HybridResult *other)
add additional toy-MC experiments to the current results use the data test statistics of the added ob...
HypoTestCalculator is an interface class for a tools which produce RooStats HypoTestResults.
virtual HypoTestResult * GetHypoTest() const =0
This class is now deprecated and to be replaced by the HypoTestInverter.
bool RunFixedScan(int nBins, double xMin, double xMax)
bool RunAutoScan(double xMin, double xMax, double target, double epsilon=0.005, unsigned int numAlgorithm=0)
virtual Double_t ConfidenceLevel() const
Get the Confidence level for the test.
bool RunOnePoint(double thisX)
run only one point
virtual Double_t Size() const
Get the size of the test (eg. rate of Type I error)
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
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
int ArraySize() const
number of entries in the results array
bool fInterpolateLowerLimit
two sided scan (look for lower/upper limit)
std::vector< double > fXValues
number of points used to build expected p-values
void UseCLs(bool on=true)
flag to switch between using CLsb (default) or CLs as confidence level
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results
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...
HypoTestResult is a base class for results from hypothesis tests.
virtual void Add(TObject *obj)
Definition: TList.h:87
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
Basic string class.
Definition: TString.h:131
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Namespace for the RooStats classes.
Definition: Asimov.h:20
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617