Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Asymptotics.cxx
Go to the documentation of this file.
1#include "xRooFit/xRooFit.h"
2
3#include <cmath>
4#include "Math/ProbFunc.h"
7
9
11
12double xRooFit::Asymptotics::k(const IncompatFunc &compatRegions, double pValue, double poiVal, double poiPrimeVal,
13 double sigma, double low, double high)
14{
15
16 // determine the pll value corresponding to nSigma expected - i.e. where the altPValue equals e.g. 50% for nSigma=0,
17 // find the solution (wrt x) of: FitManager::altPValue(x, var(poi), alt_val, _sigma_mu, _compatibilityFunction) -
18 // targetPValue = 0
19 double targetTailIntegral = pValue; // ROOT::Math::normal_cdf(*nSigma);
20
21 // check how much of the alt distribution density is in the delta function @ 0
22 // if more than 1 - target is in there, if so then pll must be 0
23 double prob_in_delta = Phi_m(poiVal, poiPrimeVal, std::numeric_limits<double>::infinity(), sigma, compatRegions);
24 // also get a contribution to the delta function for mu_hat < mu_L IF mu==mu_L
25 if (poiVal == low) {
26 // since mu_hat is gaussian distributed about mu_prime with std-dev = sigma
27 // the integral is Phi( mu_L - mu_prime / (sigma) )
28 double mu_L = low;
29 prob_in_delta += ROOT::Math::normal_cdf((mu_L - poiPrimeVal) / sigma);
30 }
31
32 if (prob_in_delta > 1 - targetTailIntegral) {
33 return 0;
34 }
35
36 struct TailIntegralFunction {
37 TailIntegralFunction(double _poiVal, double _alt_val, double _sigma_mu, double _low, double _high,
38 IncompatFunc _compatibilityFunction, double _target)
39 : poiVal(_poiVal),
40 alt_val(_alt_val),
41 sigma_mu(_sigma_mu),
42 low(_low),
43 high(_high),
44 target(_target),
45 cFunc(_compatibilityFunction)
46 {
47 }
48
49 double operator()(double x) const
50 {
51 double val = PValue(cFunc, x, poiVal, alt_val, sigma_mu, low, high);
52 if (val < 0)
53 kInvalid = true;
54 return val - target;
55 }
56
57 double poiVal, alt_val, sigma_mu, low, high, target;
58 IncompatFunc cFunc;
59 mutable bool kInvalid = false;
60 };
61
62 TailIntegralFunction f(poiVal, poiPrimeVal, sigma, low, high, compatRegions, targetTailIntegral);
65
66 auto tmpLvl = gErrorIgnoreLevel;
68 double _pll = 500.;
69 double currVal(1.);
70 int tryCount(0);
71 double _prev_pll = _pll;
72 do {
73 currVal = wf(_pll);
74 if (currVal > 1e-4) {
75 _pll = 2. * (_pll + 1.); // goto bigger pll scale
76 } else if (currVal < -1e-4) {
77 _pll /= 2.; // goto smaller pll scale
78 }
79 // std::cout << "pll = " << _pll << " currVal = " << currVal << std::endl;
80 brf.SetFunction(wf, 0, _pll);
81 if (brf.Solve()) {
82 _prev_pll = _pll;
83 _pll = brf.Root();
84 }
85 if (f.kInvalid) { // happens if problem evaluating PValue (e.g. sigma was nan)
86 _pll = std::numeric_limits<double>::quiet_NaN();
87 break;
88 }
89 // std::cout << " -- " << brf.Root() << " " << FitManager::altPValue(_pll, mu, alt_val, sigma, pllModifier()) << "
90 // >> " << wf(_pll) << std::endl;
91 tryCount++;
92 if (tryCount > 20) {
93 gErrorIgnoreLevel = tmpLvl;
94 Warning("Asymptotics::k", "Reached limit pValue=%g pll=%g", pValue, _pll);
95 break;
96 }
97 } while (std::abs(wf(_pll)) > 1e-4 && std::abs(wf(_pll)) < std::abs(wf(_prev_pll)) * 0.99);
98 gErrorIgnoreLevel = tmpLvl;
99 return _pll;
100}
101
102double xRooFit::Asymptotics::PValue(const IncompatFunc &compatRegions, double k, double poiVal, double poi_primeVal,
103 double sigma, double lowBound, double upBound)
104{
105 // uncapped test statistic is equal to onesidednegative when k is positive, and equal to 1.0 - difference between
106 // twosided and onesidednegative when k is negative ...
107 if (compatRegions == IncompatibilityFunction(Uncapped, poiVal)) {
108 // if(k==0) return 0.5;
109 if (k > 0)
110 return PValue(OneSidedNegative, k, poiVal, poi_primeVal, sigma, lowBound, upBound);
111 return 1. - (PValue(TwoSided, -k, poiVal, poi_primeVal, sigma, lowBound, upBound) -
112 PValue(OneSidedNegative, -k, poiVal, poi_primeVal, sigma, lowBound, upBound));
113 }
114
115 // if(k<0) return 1.;
116 if (k <= 0) {
117 if (compatRegions == IncompatibilityFunction(OneSidedNegative, poiVal) && std::abs(poiVal - poi_primeVal) < 1e-9)
118 return 0.5; // when doing discovery (one-sided negative) use a 0.5 pValue
119 return 1.; // case to catch the delta function that ends up at exactly 0 for the one-sided tests
120 }
121
122 if (sigma == 0) {
123 if (lowBound != -std::numeric_limits<double>::infinity() || upBound != std::numeric_limits<double>::infinity()) {
124 return -1;
125 } else if (std::abs(poiVal - poi_primeVal) > 1e-12) {
126 return -1;
127 }
128 }
129
130 // get the poi value that defines the test statistic, and the poi_prime hypothesis we are testing
131 // when setting limits, these are often the same value
132
133 double Lambda_y = 0;
134 if (std::abs(poiVal - poi_primeVal) > 1e-12)
135 Lambda_y = (poiVal - poi_primeVal) / sigma;
136
137 if (std::isnan(Lambda_y))
138 return -1;
139
140 double k_low = (lowBound == -std::numeric_limits<double>::infinity()) ? std::numeric_limits<double>::infinity()
141 : pow((poiVal - lowBound) / sigma, 2);
142 double k_high = (upBound == std::numeric_limits<double>::infinity()) ? std::numeric_limits<double>::infinity()
143 : pow((upBound - poiVal) / sigma, 2);
144
145 double out = Phi_m(poiVal, poi_primeVal, std::numeric_limits<double>::infinity(), sigma, compatRegions) - 1;
146
147 if (out < -1) {
148 // compatibility function is unsupported, return negative
149 return -2;
150 }
151
152 // go through the 4 'regions' ... only two of which will apply
153 if (k <= k_high) {
154 out += ROOT::Math::gaussian_cdf(sqrt(k) + Lambda_y) -
155 Phi_m(poiVal, poi_primeVal, Lambda_y + sqrt(k), sigma, compatRegions);
156 } else {
157 double Lambda_high = (poiVal - upBound) * (poiVal + upBound - 2. * poi_primeVal) / (sigma * sigma);
158 double sigma_high = 2. * (upBound - poiVal) / sigma;
159 out += ROOT::Math::gaussian_cdf((k - Lambda_high) / sigma_high) -
160 Phi_m(poiVal, poi_primeVal, (k - Lambda_high) / sigma_high, sigma, compatRegions);
161 }
162
163 if (k <= k_low) {
164 out += ROOT::Math::gaussian_cdf(sqrt(k) - Lambda_y) +
165 Phi_m(poiVal, poi_primeVal, Lambda_y - sqrt(k), sigma, compatRegions);
166 } else {
167 double Lambda_low = (poiVal - lowBound) * (poiVal + lowBound - 2. * poi_primeVal) / (sigma * sigma);
168 double sigma_low = 2. * (poiVal - lowBound) / sigma;
169 out += ROOT::Math::gaussian_cdf((k - Lambda_low) / sigma_low) +
170 Phi_m(poiVal, poi_primeVal, (Lambda_low - k) / sigma_low, sigma, compatRegions);
171 /*out += ROOT::Math::gaussian_cdf((k-Lambda_low)/sigma_low) +
172 2*Phi_m(poiVal,poi_primeVal,(Lambda_low - k_low)==0 ? 0 : ((Lambda_low -
173 k_low)/sigma_low),sigma,compatRegions)
174 - Phi_m(poiVal,poi_primeVal,(Lambda_low - k)/sigma_low,sigma,compatFunc);
175*/
176
177 // handle case where poiVal = lowBound (e.g. testing mu=0 when lower bound is mu=0).
178 // sigma_low will be 0 and gaussian_cdf will end up being 1, but we need it to converge instead
179 // to 0.5 so that pValue(k=0) converges to 1 rather than 0.5.
180 // handle this by 'adding' back in the lower bound
181 // TODO: Think more about this?
182 /*if (sigma_low == 0) {
183 out -= 0.5;
184 }*/
185 }
186
187 return 1. - out;
188}
189
190double
191xRooFit::Asymptotics::Phi_m(double /*mu*/, double mu_prime, double a, double sigma, const IncompatFunc &compatRegions)
192{
193
194 if (sigma == 0)
195 sigma = 1e-100; // avoid nans if sigma is 0
196
197 // want to evaluate gaussian integral in regions where IncompatFunc = 0
198
199 double out = 0;
200 double lowEdge = std::numeric_limits<double>::quiet_NaN();
201 for (auto &transition : compatRegions) {
202 if (transition.first >= (a * sigma + mu_prime))
203 break;
204 if (transition.second == 0 && std::isnan(lowEdge)) {
205 lowEdge = transition.first;
206 } else if (!std::isnan(lowEdge)) {
207 out += ROOT::Math::gaussian_cdf((transition.first - mu_prime) / sigma) -
208 ROOT::Math::gaussian_cdf((lowEdge - mu_prime) / sigma);
209 lowEdge = std::numeric_limits<double>::quiet_NaN();
210 }
211 }
212 if (!std::isnan(lowEdge)) { // also catches case where last transition is before a
213 out += ROOT::Math::gaussian_cdf(a) - ROOT::Math::gaussian_cdf((lowEdge - mu_prime) / sigma);
214 }
215
216 return out;
217}
218
219int xRooFit::Asymptotics::CompatFactor(const IncompatFunc &func, double mu_hat)
220{
221 if (std::isnan(mu_hat))
222 return 1; // nan is never compatible
223 int out = 1;
224 for (auto &transition : func) {
225 if (transition.first > mu_hat)
226 break;
227 out = transition.second;
228 }
229 return out;
230}
231
232// RooRealVar xRooFit::Asymptotics::FindLimit(TGraph* pVals, double target_pVal) {
233// auto target_sig = RooStats::PValueToSignificance(target_pVal);
234// TGraph sig;
235// for(int i=0;i<pVals->GetN();i++) {
236// sig.SetPoint(i,pVals->GetPointX(i),RooStats::PValueToSignificance(pVals->GetPointY(i))-target_sig);
237// }
238// sig.Sort(); // ensure points are in x order
239// // now loop over and find where function crosses zero
240// for(int i=0;i<sig.GetN();i++) {
241// if (sig.GetPointY(i)>=0) {
242// if (i==0) return RooRealVar("limit","limit",std::numeric_limits<double>::quiet_NaN());
243// double prev_x = sig.GetPointX(i-1);
244// double next_x = sig.GetPointX(i);
245// double prev_y = sig.GetPointY(i-1);
246// double next_y = sig.GetPointY(i);
247// return RooRealVar("limit","limit", prev_x + (next_x-prev_x)*(-prev_y)/(next_y-prev_y));
248// }
249// }
250// return RooRealVar("limit","limit",std::numeric_limits<double>::quiet_NaN());
251// }
252
#define f(i)
Definition RSha256.hxx:104
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
constexpr Int_t kFatal
Definition TError.h:49
Int_t gErrorIgnoreLevel
Error handling routines.
Definition TError.cxx:31
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:229
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t target
TRObject operator()(const T1 &t1) const
@ kInvalid
Definition TSystem.h:77
std::vector< std::pair< double, int > > IncompatFunc
Definition xRooFit.h:142
Class for finding the root of a one dimensional function using the Brent algorithm.
bool SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup) override
Sets the function for the rest of the algorithms.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10) override
Returns the X value corresponding to the function value fy for (xmin<x<xmax).
double Root() const override
Returns root value.
Template class to wrap any C++ callable object which takes one argument i.e.
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition RVec.hxx:1846
const Double_t sigma
Double_t x[n]
Definition legend1.C:17
double gaussian_cdf(double x, double sigma=1, double x0=0)
Alternative name for same function.
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
#define BEGIN_XROOFIT_NAMESPACE
Definition Config.h:24
#define END_XROOFIT_NAMESPACE
Definition Config.h:25