36 using namespace RooStats;
37 using namespace RooFit;
42 double HypoTestInverterResult::fgAsymptoticMaxSigma = 5;
45 HypoTestInverterResult::HypoTestInverterResult(const
char *
name ) :
49 fInterpolateLowerLimit(true),
50 fInterpolateUpperLimit(true),
51 fFittedLowerLimit(false),
52 fFittedUpperLimit(false),
53 fInterpolOption(kLinear),
56 fCLsCleanupThreshold(0.005)
66 fUseCLs(other.fUseCLs),
67 fIsTwoSided(other.fIsTwoSided),
68 fInterpolateLowerLimit(other.fInterpolateLowerLimit),
69 fInterpolateUpperLimit(other.fInterpolateUpperLimit),
70 fFittedLowerLimit(other.fFittedLowerLimit),
71 fFittedUpperLimit(other.fFittedUpperLimit),
72 fInterpolOption(other.fInterpolOption),
73 fLowerLimitError(other.fLowerLimitError),
74 fUpperLimitError(other.fUpperLimitError),
75 fCLsCleanupThreshold(other.fCLsCleanupThreshold)
83 for (
int i = 0; i < nOther; ++i)
115 for (
int i=0; i < nOther; ++i) {
133 fInterpolateLowerLimit(true),
134 fInterpolateUpperLimit(true),
135 fFittedLowerLimit(
false),
136 fFittedUpperLimit(
false),
137 fInterpolOption(kLinear),
138 fLowerLimitError(-1),
139 fUpperLimitError(-1),
140 fCLsCleanupThreshold(0.005)
147 fParameters.removeAll();
148 fParameters.takeOwnership();
163 std::vector<double> qv;
172 bool resultIsAsymptotic(
false);
177 resultIsAsymptotic =
true;
181 int nPointsRemoved(0);
183 double CLsobsprev(1.0);
184 std::vector<double>::iterator itr =
fXValues.begin();
201 if (resultIsAsymptotic) {
203 double dsig = 2.*maxSigma / (values.size() -1) ;
204 int i0 = (int)
TMath::Floor ( ( -nsig2 + maxSigma )/dsig + 0.5 );
205 int i1 = (int)
TMath::Floor ( ( -nsig1 + maxSigma )/dsig + 0.5 );
206 int i2 = (int)
TMath::Floor ( ( maxSigma )/dsig + 0.5 );
207 int i3 = (int)
TMath::Floor ( ( nsig1 + maxSigma )/dsig + 0.5 );
208 int i4 = (int)
TMath::Floor ( ( nsig2 + maxSigma )/dsig + 0.5 );
216 double *
z =
const_cast<double *
>( &values[0] );
224 for (
int j=0; j<5; ++j) { qv[j]=q[j]; }
231 double CLsobs = qv[5];
235 bool removeThisPoint(
false);
238 if (!removeThisPoint && resultIsAsymptotic && i>=1 && CLsobs>CLsobsprev) {
240 removeThisPoint =
true;
241 }
else { CLsobsprev = CLsobs; }
244 if (!removeThisPoint && i>=1 && CLsobs>=0.9999) {
246 removeThisPoint =
true;
252 if (removeThisPoint) {
269 return nPointsRemoved;
295 if (nOther == 0)
return true;
303 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging result from " << otherResult.
GetName()
304 <<
" in " <<
GetName() << std::endl;
309 if (addExpPValues || mergeExpPValues)
310 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - merging also the expected p-values from pseudo-data" << std::endl;
317 for (
int i = 0; i < nOther; ++i)
324 for (
int i = 0; i < nOther; ++i) {
325 double otherVal = otherResult.
fXValues[i];
327 if (otherHTR == 0)
continue;
328 bool sameXFound =
false;
329 for (
int j = 0; j < nThis; ++j) {
336 thisHTR->
Append(otherHTR);
338 if (mergeExpPValues) {
343 if (thisNToys != otherNToys )
344 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::Add expexted p values have been generated with different toys " << thisNToys <<
" , " << otherNToys << std::endl;
363 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new number of points is " <<
fXValues.size()
366 oocoutI(
this,
Eval) <<
"HypoTestInverterResult::Add - new toys/point is "
386 if (!r)
return false;
513 const double tol = 1.E-12;
524 struct InterpolatedGraph {
525 InterpolatedGraph(
const TGraph &
g,
double target,
const char * interpOpt) :
526 fGraph(g), fTarget(target), fInterpOpt(interpOpt) {}
530 return fGraph.Eval(x, (
TSpline*) 0, fInterpOpt) - fTarget;
532 const TGraph & fGraph;
547 InterpolatedGraph
f(graph,y0,opt);
552 const double * y = graph.GetY();
553 int n = graph.GetN();
555 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult::GetGraphX - need at least 2 points for interpolation (n=" << n <<
")\n";
556 return (n>0) ? y[0] : 0;
570 if (axmin >= axmax ) {
571 xmin = graph.GetX()[0];
572 xmax = graph.GetX()[n-1];
580 if ( (ymax < y0 && !lowSearch) || ( ymin > y0 && lowSearch) ) {
584 if ( (ymax < y0 && lowSearch) || ( ymin > y0 && !lowSearch) ) {
590 #ifdef ISREALLYNEEDED //??
592 bool isCross =
false;
595 for (
int i = 0; i<
n; ++i) {
597 graph.GetPoint(i,xp,yp);
598 if (xp >= xmin && xp <= xmax) {
600 if (prod < 0 && !first) {
608 return (lowSearch) ? varmin : varmax;
617 bool ret = brf.
Solve(100, 1.
E-16, 1.
E-6);
619 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - interpolation failed - return inf" << std::endl;
622 double limit = brf.
Root();
626 if (lowSearch) std::cout <<
"lower limit search : ";
627 else std::cout <<
"Upper limit search : ";
628 std::cout <<
"do interpolation between " << xmin <<
" and " << xmax <<
" limit is " << limit << std::endl;
633 if (axmin >= axmax) {
636 std::cout <<
"do new interpolation dividing from " << index <<
" and " << y[index] << std::endl;
639 if (lowSearch && index >= 1 && (y[0] - y0) * ( y[index]- y0) < 0) {
641 limit =
GetGraphX(graph, y0, lowSearch, graph.GetX()[0], graph.GetX()[index] );
643 else if (!lowSearch && index < n-2 && (y[n-1] - y0) * ( y[index+1]- y0) < 0) {
645 limit =
GetGraphX(graph, y0, lowSearch, graph.GetX()[index+1], graph.GetX()[n-1] );
661 <<
"Interpolate the upper limit between the 2 results closest to the target confidence level"
674 double val = (lowSearch) ? xmin : xmax;
675 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::FindInterpolatedLimit"
676 <<
" - not enough points to get the inverted interval - return "
685 std::vector<unsigned int> index(n );
689 for (
int i = 0; i <
n; ++i)
701 double * itrmax = std::max_element(graph.GetY() , graph.GetY() +
n);
702 double ymax = *itrmax;
703 int iymax = itrmax - graph.GetY();
704 double xwithymax = graph.GetX()[iymax];
713 xmin = ( graph.GetY()[0] <= target ) ? graph.GetX()[0] : varmin;
726 xmax = ( graph.GetY()[n-1] <= target ) ? graph.GetX()[n-1] : varmax;
740 if (iymax <= (n-1)/2 ) {
751 std::cout <<
" found xmin, xmax = " << xmin <<
" " << xmax <<
" for search " << lowSearch << std::endl;
762 if (upI < 1)
return xmin;
768 if (lowI >= n-1)
return xmax;
776 std::cout <<
"finding " << lowSearch <<
" limit betweem " << xmin <<
" " << xmax << endl;
779 double limit =
GetGraphX(graph, target, lowSearch, xmin, xmax);
785 std::cout <<
"limit is " << limit << std::endl;
799 double limit2 =
GetGraphX(graph, target, !lowSearch, xmin, xmax);
806 std::cout <<
"other limit is " << limit2 << std::endl;
825 int closestIndex = -1;
827 double smallestError = 2;
828 double bestValue = 2;
837 if ( dist < bestValue) {
842 if (bestIndex >=0)
return bestIndex;
850 std::vector<unsigned int> indx(n);
852 std::vector<double> xsorted( n);
853 for (
int i = 0; i <
n; ++i) xsorted[i] =
fXValues[indx[i] ];
857 std::cout <<
"finding closest point to " << xtarget <<
" is " << index1 <<
" " << indx[index1] << std::endl;
861 if (index1 < 0)
return indx[0];
862 if (index1 >= n-1)
return indx[n-1];
863 int index2 = index1 +1;
865 if (mode == 2)
return (
GetXValue(indx[index1]) <
GetXValue(indx[index2])) ? indx[index1] : indx[index2];
866 if (mode == 3)
return (
GetXValue(indx[index1]) >
GetXValue(indx[index2])) ? indx[index1] : indx[index2];
909 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError"
910 <<
"Empty result \n";
915 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimateError"
916 <<
" only points - return its error\n";
926 std::cout <<
"calculate estimate error " << type <<
" between " << xmin <<
" and " << xmax << std::endl;
931 std::vector<unsigned int> indx(
fXValues.size());
942 graph.SetPointError(ip, 0.,
GetYError(indx[i]) );
947 if (graph.GetN() < 2) {
948 if (
np >= 2)
oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimatedError - no valid points - cannot estimate the " << type <<
" limit error " << std::endl;
961 TF1 fct(
"fct",
"exp([0] * (x - [2] ) + [1] * (x-[2])**2)", minX, maxX);
962 double scale = maxX-minX;
964 fct.SetParameters( 2./scale, 0.1/scale, graph.GetX()[0] );
965 fct.SetParLimits(0,0,100./scale);
966 fct.SetParLimits(1,0, 10./scale); }
968 fct.SetParameters( -2./scale, -0.1/scale );
969 fct.SetParLimits(0,-100./scale, 0);
970 fct.SetParLimits(1,-100./scale, 0); }
972 if (graph.GetN() < 3) fct.FixParameter(1,0.);
980 TCanvas *
c1 =
new TCanvas();
981 std::cout <<
"fitting for limit " << type <<
"between " << minX <<
" , " << maxX <<
" points considered " << graph.GetN() << std::endl;
982 int fitstat = graph.Fit(&fct,
" EX0");
983 graph.SetMarkerStyle(20);
989 int fitstat = graph.Fit(&fct,
"Q EX0");
997 double m = fct.Derivative(
GetXValue(index) );
1002 oocoutW(
this,
Eval) <<
"HypoTestInverterResult::CalculateEstimatedError - cannot estimate the " << type <<
" limit error " << std::endl;
1011 std::cout <<
"closes point to the limit is " << index <<
" " <<
GetXValue(index) <<
" and has error " <<
GetYError(index) << std::endl;
1040 if (!firstResult)
return 0;
1047 if (!result)
return 0;
1054 if (index < 0 || index >=
ArraySize() )
return 0;
1060 static bool useFirstB =
false;
1062 int bIndex = (useFirstB) ? 0 : index;
1069 if (bDistribution && sbDistribution) {
1080 for (
int i = 0; i < bDistribution->
GetSize(); ++i) {
1089 const int npoints = 11;
1091 const double dsig = 2* smax/ (npoints-1) ;
1092 std::vector<double>
values(npoints);
1093 for (
int i = 0; i < npoints; ++i) {
1094 double nsig = -smax + dsig*i;
1096 if (pval < 0) {
return 0;}
1099 return new SamplingDistribution(
"Asymptotic expected values",
"Asymptotic expected values",values);
1111 oocoutE(
this,
Eval) <<
"HypoTestInverterResult::GetLimitDistribution"
1112 <<
" not enought points - return 0 " << std::endl;
1116 ooccoutD(
this,
Eval) <<
"HypoTestInverterResult - computing limit distribution...." << std::endl;
1121 std::vector<SamplingDistribution*> distVec( npoints );
1123 for (
unsigned int i = 0; i < distVec.size(); ++i) {
1128 distVec[i]->InverseCDF(0);
1129 sum += distVec[i]->GetSize();
1132 int size = int( sum/ npoints);
1135 ooccoutW(
this,
InputArguments) <<
"HypoTestInverterResult - set a minimum size of 10 for limit distribution" << std::endl;
1143 std::vector< std::vector<double> > quantVec(npoints );
1144 for (
int i = 0; i < npoints; ++i) {
1146 if (!distVec[i])
continue;
1149 std::vector<double> pvalues = distVec[i]->GetSamplingDistribution();
1150 delete distVec[i]; distVec[i] = 0;
1151 std::sort(pvalues.begin(), pvalues.end());
1156 quantVec[i] = std::vector<double>(size);
1157 for (
int ibin = 0; ibin < size; ++ibin) {
1159 p[0] =
std::min( (ibin+1) * 1./
double(size), 1.0);
1162 (quantVec[i])[ibin] = q[0];
1168 std::vector<unsigned int> index( npoints );
1175 std::vector<double> limits(size);
1177 for (
int j = 0; j < size; ++j ) {
1180 for (
int k = 0; k < npoints ; ++k) {
1181 if (quantVec[index[k]].size() > 0 )
1182 g.SetPoint(g.GetN(),
GetXValue(index[k]), (quantVec[index[k]])[j] );
1185 limits[j] =
GetGraphX(g, target, lower);
1241 if (nEntries <= 0)
return (lower) ? 1 : 0;
1249 if (!limitDist)
return 0;
1251 if (values.size() <= 1)
return 0;
1273 std::vector<unsigned int> index(nEntries);
1276 for (
int j=0; j<nEntries; ++j) {
1280 ooccoutI(
this,
Eval) <<
"HypoTestInverterResult - cannot compute expected p value distribution for point, x = "
1281 <<
GetXValue(i) <<
" skip it " << std::endl;
1285 double * x =
const_cast<double *
>(&values[0]);
1287 g.SetPoint(g.GetN(),
fXValues[i], q[0] );
1291 ooccoutE(
this,
Eval) <<
"HypoTestInverterResult - cannot compute limits , not enough points, n = " << g.GetN() << std::endl;
1302 if (!limitDist)
return 0;
1304 double * x =
const_cast<double *
>(&values[0]);
SamplingDistribution * GetLimitDistribution(bool lower) const
void SetBackgroundAsAlt(Bool_t l=kTRUE)
double dist(Rotation3D const &r1, Rotation3D const &r2)
void SetAltDistribution(SamplingDistribution *alt)
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Double_t Floor(Double_t x)
virtual Double_t ConfidenceLevel() const
return confidence level
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
static double GetExpectedPValues(double pnull, double palt, double nsigma, bool usecls, bool oneSided=true)
function given the null and the alt p value - return the expected one given the N - sigma value ...
double CLbError(int index) const
return the observed CLb value for the i-th entry
HypoTestInverterResult(const char *name=0)
default constructor
int ArraySize() const
number of entries in the results array
virtual Double_t CLs() const
CLs is simply CLs+b/CLb (not a method, but a quantity)
int FindClosestPointIndex(double target, int mode=0, double xtarget=0)
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
Base class for spline implementation containing the Draw/Paint methods //.
std::vector< double > values
HypoTestResult is a base class for results from hypothesis tests.
void ToUpper()
Change string to upper case.
SamplingDistribution * GetSignalAndBackgroundTestStatDist(int index) const
virtual Double_t getMin(const char *name=0) const
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...
SamplingDistribution * GetAltDistribution(void) const
Int_t GetSize() const
size of samples
virtual Double_t CLsplusb() const
Convert AlternatePValue into a "confidence level".
ClassImp(TIterator) Bool_t TIterator return false
Compare two iterator objects.
SamplingDistribution * GetBackgroundTestStatDist(int index) const
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
virtual TObject * RemoveAt(Int_t idx)
Template class to wrap any C++ callable object which takes one argument i.e.
const char * Data() const
SimpleInterval & operator=(const SimpleInterval &other)
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
double FindInterpolatedLimit(double target, bool lowSearch=false, double xmin=1, double xmax=0)
double CLsError(int index) const
return the observed CLb value for the i-th entry
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities Parameters: x -the data sample n ...
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Double_t fConfidenceLevel
virtual ~HypoTestInverterResult()
destructor
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Double_t LowerLimit()
lower and upper bound of the confidence interval (to get upper/lower limits, multiply the size( = 1-c...
bool fInterpolateLowerLimit
two sided scan (look for lower/upper limit)
int ExclusionCleanup()
remove points that appear to have failed.
double GetGraphX(const TGraph &g, double y0, bool lowSearch, double &xmin, double &xmax) const
HypoTestResult * GetResult(int index) const
return a pointer to the i^th result object
RooRealVar represents a fundamental (non-derived) real valued object.
double CLsplusb(int index) const
return the observed CLsplusb value for the i-th entry
double Root() const
Returns root value.
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
SamplingDistribution * GetExpectedPValueDist(int index) const
return expected distribution of p-values (Cls or Clsplusb)
double GetExpectedLowerLimit(double nsig=0, const char *opt="") const
get Limit value correspnding at the desired nsigma level (0) is median -1 sigma is 1 sigma ...
double fCLsCleanupThreshold
TList fExpPValues
list of HypoTestResult for each point
virtual const char * GetName() const
Returns name of object.
Bool_t GetPValueIsRightTail(void) const
void SetTestStatisticData(const Double_t tsd)
This class simply holds a sampling distribution of some test statistic.
ClassImp(RooStats::HypoTestInverterResult) using namespace RooStats
virtual void Append(const HypoTestResult *other)
add values from another HypoTestResult
static double fgAsymptoticMaxSigma
Class for finding the root of a one dimensional function using the Brent algorithm.
void SetPValueIsRightTail(Bool_t pr)
double GetXValue(int index) const
function to return the value of the parameter of interest for the i^th entry in the results ...
TRObject operator()(const T1 &t1) const
virtual Int_t GetSize() const
double CLs(int index) const
return the observed CLb value for the i-th entry
HypoTestInverterResult class: holds the array of hypothesis test results and compute a confidence int...
Double_t LowerLimitEstimatedError()
rough estimation of the error on the computed bound of the confidence interval Estimate of lower limi...
double GetYValue(int index) const
function to return the value of the confidence level for the i^th entry in the results ...
T MaxElement(Long64_t n, const T *a)
InterpolOption_t fInterpolOption
SamplingDistribution * GetNullDistribution(void) const
double CLsplusbError(int index) const
return the observed CLsplusb value for the i-th entry
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
virtual Double_t getMax(const char *name=0) const
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Double_t UpperLimitEstimatedError()
Estimate of lower limit error function evaluates only a rought error on the lower limit...
virtual TObject * clone(const char *newname) const
double GetExpectedUpperLimit(double nsig=0, const char *opt="") const
get Limit value correspnding at the desired nsigma level (0) is median -1 sigma is 1 sigma ...
HypoTestInverterResult & operator=(const HypoTestInverterResult &other)
operator =
virtual void Add(TObject *obj)
double fLowerLimitError
interpolatation option (linear or spline)
virtual void RemoveAll(TCollection *col)
Remove all objects in collection col from this collection.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
double GetExpectedLimit(double nsig, bool lower, const char *opt="") const
bool Add(const HypoTestInverterResult &otherResult)
merge with the content of another HypoTestInverterResult object
int FindIndex(double xvalue) const
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
void SetNpx(int npx)
Set the number of point used to bracket root using a grid.
double CalculateEstimatedError(double target, bool lower=true, double xmin=1, double xmax=0)
double CLb(int index) const
return the observed CLb value for the i-th entry
bool fInterpolateUpperLimit
Long64_t BinarySearch(Long64_t n, const T *array, T value)
void SetNullDistribution(SamplingDistribution *null)
void SortItr(Iterator first, Iterator last, IndexIterator index, Bool_t down=kTRUE)
SamplingDistribution * GetNullTestStatDist(int index) const
same in terms of alt and null
std::vector< double > fXValues
max sigma value used to scan asymptotic expected p values
T MinElement(Long64_t n, const T *a)
Bool_t GetBackGroundIsAlt(void) const