CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BinomialNoncentralInterval.h
Go to the documentation of this file.
1 #ifndef PhysicsTools_RooStatsCms_BinomialNoncentralinterval_h
2 #define PhysicsTools_RooStatsCms_BinomialNoncentralinterval_h
3 /* \class BinomialNoncentralinterval
4  *
5  * \author Jordan Tucker
6  *
7  * integration in CMSSW: Luca Listta
8  *
9  */
10 
11 #include <algorithm>
12 #include <cmath>
13 #include <vector>
14 
15 #include "Math/PdfFuncMathCore.h"
16 
17 #if (defined (STANDALONE) or defined (__CINT__) )
18 #include "BinomialInterval.h"
19 #include "BinomialProbHelper.h"
20 #else
23 #endif
24 
25 // Implement noncentral binomial confidence intervals using the Neyman
26 // construction. The Sorter class gives the ordering of points,
27 // i.e. it must be a functor implementing a greater-than relationship
28 // between two prob_helper instances. See feldman_cousins for an
29 // example.
30 template <typename Sorter>
32  public:
33  // Given a true value of rho and ntot trials, calculate the
34  // acceptance set [x_l, x_r] for use in a Neyman construction.
35  bool find_rho_set(const double rho, const int ntot, int& x_l, int& x_r) const {
36  // Get the binomial probabilities for every x = 0..n, and sort them
37  // in decreasing order, determined by the Sorter class.
38  std::vector<BinomialProbHelper> probs;
39  for (int i = 0; i <= ntot; ++i)
40  probs.push_back(BinomialProbHelper(rho, i, ntot));
41  std::sort(probs.begin(), probs.end(), sorter_);
42 
43  // Add up the probabilities until the total is 1 - alpha or
44  // bigger, adding the biggest point first, then the next biggest,
45  // etc. "Biggest" is given by the Sorter class and is taken care
46  // of by the sort above. JMTBAD need to find equal probs and use
47  // the sorter to differentiate between them.
48  const double target = 1 - alpha_;
49  // An invalid interval.
50  x_l = ntot;
51  x_r = 0;
52  double sum = 0;
53  for (int i = 0; i <= ntot && sum < target; ++i) {
54  sum += probs[i].prob();
55  const int& x = probs[i].x();
56  if (x < x_l) x_l = x;
57  if (x > x_r) x_r = x;
58  }
59 
60  return x_l <= x_r;
61  }
62 
63  // Construct nrho acceptance sets in rho = [0,1] given ntot trials
64  // and put the results in already-allocated x_l and x_r.
65  bool neyman(const int ntot, const int nrho, double* rho, double* x_l, double* x_r) {
66  int xL, xR;
67  for (int i = 0; i < nrho; ++i) {
68  rho[i] = double(i)/nrho;
69  find_rho_set(rho[i], ntot, xL, xR);
70  x_l[i] = xL;
71  x_r[i] = xR;
72  }
73  return true;
74  }
75 
76  // Given X successes and n trials, calculate the interval using the
77  // rho acceptance sets implemented above.
78  void calculate(const double X, const double n) {
79  set(0, 1);
80 
81  const double tol = 1e-9;
82  double rho_min, rho_max, rho;
83  int x_l, x_r;
84 
85  // Binary search for the smallest rho whose acceptance set has right
86  // endpoint X; this is the lower endpoint of the rho interval.
87  rho_min = 0; rho_max = 1;
88  while (std::fabs(rho_max - rho_min) > tol) {
89  rho = (rho_min + rho_max)/2;
90  find_rho_set(rho, int(n), x_l, x_r);
91  if (x_r < X)
92  rho_min = rho;
93  else
94  rho_max = rho;
95  }
96  lower_ = rho;
97 
98  // Binary search for the largest rho whose acceptance set has left
99  // endpoint X; this is the upper endpoint of the rho interval.
100  rho_min = 0; rho_max = 1;
101  while (std::fabs(rho_max - rho_min) > tol) {
102  rho = (rho_min + rho_max)/2;
103  find_rho_set(rho, int(n), x_l, x_r);
104  if (x_l > X)
105  rho_max = rho;
106  else
107  rho_min = rho;
108  }
109  upper_ = rho;
110  }
111 
112  private:
113  Sorter sorter_;
114 
115 #if (defined (STANDALONE) or defined (__CINT__) )
116 ClassDefT(BinomialNoncentralInterval,1)
117 #endif
118 };
119 
120 #endif
int i
Definition: DBlmapReader.cc:9
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
Definition: DDAxes.h:10
#define X(str)
Definition: MuonsGrabber.cc:49
void set(double l, double u)
void calculate(const double X, const double n)
x
Definition: VDTMath.h:216