ROOT   Reference Guide
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),
33 fRho_hat(double(x)/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
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 Prob() const
double LRatio() const
double Rho() const
BinomialProbHelper(double rho, int x, int n)
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
VSD Structures.
Definition: StringConv.hxx:21
bool operator()(const BinomialProbHelper &l, const BinomialProbHelper &r) const
auto * l
Definition: textangle.C:4
static long int sum(long int i)
Definition: Factory.cxx:2276