ROOT   Reference Guide
Searching...
No Matches
TEfficiencyHelper.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: L. Moneta Nov 2010
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2010 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10// helper class for binomial Neyman intervals
11// author Jordan Tucker
12// integration in CMSSW: Luca Lista
13// modified and integrated in ROOT: Lorenzo Moneta
14
15
16#ifndef TEFFiciencyHelper_h
17#define TEFFiciencyHelper_h
18
19#include <algorithm>
20#include <cmath>
21#include <vector>
22
24
25
26// Helper class impelementing the
27// binomial probability and the likelihood ratio
28// used for ordering the interval in the FeldmanCousins interval class
30public:
31 BinomialProbHelper(double rho, int x, int n)
32 : fRho(rho), fX(x), fN(n),
34 fProb(ROOT::Math::binomial_pdf(x, rho, n)) {
35 // Cache the likelihood ratio L(\rho)/L(\hat{\rho}), too.
36 if (x == 0)
37 fLRatio = pow(1 - rho, n);
38 else if (x == n)
39 fLRatio = pow(rho, n);
40 else {
41 // Impose this criterion: if any if the two terms is zero, the product is
42 // zero and not NaN.
43 const double term1 = pow(rho/fRho_hat, x);
44 const double term2 = pow((1 - rho)/(1 - fRho_hat), n - x);
45 fLRatio = (term1 == 0. || term2 == 0.) ? 0. : term1 * term2;
46 }
47 }
48
49 double Rho () const { return fRho; };
50 int X () const { return fX; };
51 int N () const { return fN; };
52 double Prob () const { return fProb; };
53 double LRatio() const { return fLRatio; };
54
55private:
56 double fRho;
57 int fX;
58 int fN;
59 double fRho_hat;
60 double fProb;
61 double fLRatio;
62};
63
64
65
66// Implement noncentral binomial confidence intervals using the Neyman
67// construction. The Sorter class gives the ordering of points,
68// i.e. it must be a functor implementing a greater-than relationship
69// between two prob_helper instances. See feldman_cousins for an
70// example.
71template <typename Sorter>
73public:
74
76 fLower(0),
77 fUpper(1),
78 fAlpha(0)
79 {}
80
81 void Init(double alpha) {
82 fAlpha = alpha;
83 }
84
85 // Given a true value of rho and ntot trials, calculate the
86 // acceptance set [x_l, x_r] for use in a Neyman construction.
87 bool Find_rho_set(const double rho, const int ntot, int& x_l, int& x_r) const {
88 // Get the binomial probabilities for every x = 0..n, and sort them
89 // in decreasing order, determined by the Sorter class.
90 std::vector<BinomialProbHelper> probs;
91 for (int i = 0; i <= ntot; ++i)
92 probs.emplace_back(BinomialProbHelper(rho, i, ntot));
93 std::sort(probs.begin(), probs.end(), fSorter);
94
95 // Add up the probabilities until the total is 1 - alpha or
96 // bigger, adding the biggest point first, then the next biggest,
97 // etc. "Biggest" is given by the Sorter class and is taken care
98 // of by the sort above. JMTBAD need to find equal probs and use
99 // the sorter to differentiate between them.
100 const double target = 1 - fAlpha;
101 // An invalid interval.
102 x_l = ntot;
103 x_r = 0;
104 double sum = 0;
105 for (int i = 0; i <= ntot && sum < target; ++i) {
106 sum += probs[i].Prob();
107 const int& x = probs[i].X();
108 if (x < x_l) x_l = x;
109 if (x > x_r) x_r = x;
110 }
111
112 return x_l <= x_r;
113 }
114
115 // Construct nrho acceptance sets in rho = [0,1] given ntot trials
116 // and put the results in already-allocated x_l and x_r.
117 bool Neyman(const int ntot, const int nrho, double* rho, double* x_l, double* x_r) {
118 int xL, xR;
119 for (int i = 0; i < nrho; ++i) {
120 rho[i] = double(i)/nrho;
121 Find_rho_set(rho[i], ntot, xL, xR);
122 x_l[i] = xL;
123 x_r[i] = xR;
124 }
125 return true;
126 }
127
128 // Given X successes and n trials, calculate the interval using the
129 // rho acceptance sets implemented above.
130 void Calculate(const double X, const double n) {
131 Set(0, 1);
132
133 const double tol = 1e-9;
134 double rho_min, rho_max, rho;
135 rho_min = rho_max = rho = 0.;
136 int x_l, x_r;
137
138 // Binary search for the smallest rho whose acceptance set has right
139 // endpoint X; this is the lower endpoint of the rho interval.
140 rho_min = 0; rho_max = 1;
141 while (std::abs(rho_max - rho_min) > tol) {
142 rho = (rho_min + rho_max)/2;
143 Find_rho_set(rho, int(n), x_l, x_r);
144 if (x_r < X)
145 rho_min = rho;
146 else
147 rho_max = rho;
148 }
149 fLower = rho;
150
151 // Binary search for the largest rho whose acceptance set has left
152 // endpoint X; this is the upper endpoint of the rho interval.
153 rho_min = 0; rho_max = 1;
154 while (std::abs(rho_max - rho_min) > tol) {
155 rho = (rho_min + rho_max)/2;
156 Find_rho_set(rho, int(n), x_l, x_r);
157 if (x_l > X)
158 rho_max = rho;
159 else
160 rho_min = rho;
161 }
162 fUpper = rho;
163 }
164
165 double Lower() const { return fLower; }
166 double Upper() const { return fUpper; }
167
168private:
169 Sorter fSorter;
170
171 double fLower;
172 double fUpper;
173
174 double fAlpha;
175
176 void Set(double l, double u) { fLower = l; fUpper = u; }
177
178};
179
180
181
182
185 return l.LRatio() > r.LRatio();
186 }
187};
188
189class FeldmanCousinsBinomialInterval : public BinomialNeymanInterval<FeldmanCousinsSorter> {
190 //const char* name() const { return "Feldman-Cousins"; }
191
192};
193
194
195
196
197#endif
double
ROOT::R::TRInterface & r
Definition Object.C:4
#define e(i)
Definition RSha256.hxx:103
double pow(double, double)
bool Neyman(const int ntot, const int nrho, double *rho, double *x_l, double *x_r)
bool Find_rho_set(const double rho, const int ntot, int &x_l, int &x_r) const
void Init(double alpha)
void Set(double l, double u)
void Calculate(const double X, const double n)
double LRatio() const
BinomialProbHelper(double rho, int x, int n)
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
bool operator()(const BinomialProbHelper &l, const BinomialProbHelper &r) const
auto * l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345