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),
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/// Implement noncentral binomial confidence intervals using the Neyman
66/// construction. The Sorter class gives the ordering of points,
67/// i.e. it must be a functor implementing a greater-than relationship
68/// between two prob_helper instances. See feldman_cousins for an example.
69template <typename Sorter>
71public:
72
74 fLower(0),
75 fUpper(1),
76 fAlpha(0)
77 {}
78
79 void Init(double alpha) {
80 fAlpha = alpha;
81 }
82
83 // Given a true value of rho and ntot trials, calculate the
84 // acceptance set [x_l, x_r] for use in a Neyman construction.
85 bool Find_rho_set(const double rho, const int ntot, int& x_l, int& x_r) const {
86 // Get the binomial probabilities for every x = 0..n, and sort them
87 // in decreasing order, determined by the Sorter class.
88 std::vector<BinomialProbHelper> probs;
89 for (int i = 0; i <= ntot; ++i)
90 probs.emplace_back(BinomialProbHelper(rho, i, ntot));
91 std::sort(probs.begin(), probs.end(), fSorter);
92
93 // Add up the probabilities until the total is 1 - alpha or
94 // bigger, adding the biggest point first, then the next biggest,
95 // etc. "Biggest" is given by the Sorter class and is taken care
96 // of by the sort above. JMTBAD need to find equal probs and use
97 // the sorter to differentiate between them.
98 const double target = 1 - fAlpha;
99 // An invalid interval.
100 x_l = ntot;
101 x_r = 0;
102 double sum = 0;
103 for (int i = 0; i <= ntot && sum < target; ++i) {
104 sum += probs[i].Prob();
105 const int& x = probs[i].X();
106 if (x < x_l) x_l = x;
107 if (x > x_r) x_r = x;
108 }
109
110 return x_l <= x_r;
111 }
112
113 // Construct nrho acceptance sets in rho = [0,1] given ntot trials
114 // and put the results in already-allocated x_l and x_r.
115 bool Neyman(const int ntot, const int nrho, double* rho, double* x_l, double* x_r) {
116 int xL, xR;
117 for (int i = 0; i < nrho; ++i) {
118 rho[i] = double(i)/nrho;
119 Find_rho_set(rho[i], ntot, xL, xR);
120 x_l[i] = xL;
121 x_r[i] = xR;
122 }
123 return true;
124 }
125
126 // Given X successes and n trials, calculate the interval using the
127 // rho acceptance sets implemented above.
128 void Calculate(const double X, const double n) {
129 Set(0, 1);
130
131 const double tol = 1e-9;
132 double rho_min, rho_max, rho;
133 rho_min = rho_max = rho = 0.;
134 int x_l, x_r;
135
136 // Binary search for the smallest rho whose acceptance set has right
137 // endpoint X; this is the lower endpoint of the rho interval.
138 rho_min = 0; rho_max = 1;
139 while (std::abs(rho_max - rho_min) > tol) {
140 rho = (rho_min + rho_max)/2;
141 Find_rho_set(rho, int(n), x_l, x_r);
142 if (x_r < X)
143 rho_min = rho;
144 else
145 rho_max = rho;
146 }
147 fLower = rho;
148
149 // Binary search for the largest rho whose acceptance set has left
150 // endpoint X; this is the upper endpoint of the rho interval.
151 rho_min = 0; rho_max = 1;
152 while (std::abs(rho_max - rho_min) > tol) {
153 rho = (rho_min + rho_max)/2;
154 Find_rho_set(rho, int(n), x_l, x_r);
155 if (x_l > X)
156 rho_max = rho;
157 else
158 rho_min = rho;
159 }
160 fUpper = rho;
161 }
162
163 double Lower() const { return fLower; }
164 double Upper() const { return fUpper; }
165
166private:
167 Sorter fSorter;
168
169 double fLower;
170 double fUpper;
171
172 double fAlpha;
173
174 void Set(double l, double u) { fLower = l; fUpper = u; }
175
176};
177
178
179
180
183 return l.LRatio() > r.LRatio();
184 }
185};
186
187class FeldmanCousinsBinomialInterval : public BinomialNeymanInterval<FeldmanCousinsSorter> {
188 //const char* name() const { return "Feldman-Cousins"; }
189
190};
191
192
193
194
195#endif
#define e(i)
Definition: RSha256.hxx:103
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
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 r
Implement noncentral binomial confidence intervals using the Neyman construction.
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)
Helper class impelementing the binomial probability and the likelihood ratio used for ordering the in...
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.
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1778
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1792
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
bool operator()(const BinomialProbHelper &l, const BinomialProbHelper &r) const
TLine l
Definition: textangle.C:4
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345