Logo ROOT   6.14/05
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 
23 #include "Math/PdfFuncMathCore.h"
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
30 public:
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 
55 private:
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.
71 template <typename Sorter>
73 public:
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 
168 private:
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 
184  bool operator()(const BinomialProbHelper& l, const BinomialProbHelper& r) const {
185  return l.LRatio() > r.LRatio();
186  }
187 };
188 
189 class FeldmanCousinsBinomialInterval : public BinomialNeymanInterval<FeldmanCousinsSorter> {
190  //const char* name() const { return "Feldman-Cousins"; }
191 
192 };
193 
194 
195 
196 
197 #endif
static long int sum(long int i)
Definition: Factory.cxx:2258
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double LRatio() const
void Set(double l, double u)
Double_t x[n]
Definition: legend1.C:17
double pow(double, double)
bool operator()(const BinomialProbHelper &l, const BinomialProbHelper &r) const
double Rho() const
ROOT::R::TRInterface & r
Definition: Object.C:4
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
bool Find_rho_set(const double rho, const int ntot, int &x_l, int &x_r) const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Namespace for new Math classes and functions.
auto * l
Definition: textangle.C:4
BinomialProbHelper(double rho, int x, int n)
bool Neyman(const int ntot, const int nrho, double *rho, double *x_l, double *x_r)
void Init(double alpha)
const Int_t n
Definition: legend1.C:16
void Calculate(const double X, const double n)
double Prob() const