Logo ROOT   6.08/07
Reference Guide
RooPoisson.cxx
Go to the documentation of this file.
1  /*****************************************************************************
2  * Project: RooFit *
3  * *
4  * Simple Poisson PDF
5  * author: Kyle Cranmer <cranmer@cern.ch>
6  * *
7  *****************************************************************************/
8 
9 /**
10 \file RooPoisson.cxx
11 \class RooPoisson
12 \ingroup Roofit
13 
14 Poisson pdf
15 **/
16 
17 #include <iostream>
18 
19 #include "RooPoisson.h"
20 #include "RooAbsReal.h"
21 #include "RooAbsCategory.h"
22 
23 #include "RooRandom.h"
24 #include "RooMath.h"
25 #include "TMath.h"
26 #include "Math/ProbFuncMathCore.h"
27 
28 #include "TError.h"
29 
30 using namespace std;
31 
33 
34 
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 /// Constructor
38 
39 RooPoisson::RooPoisson(const char *name, const char *title,
40  RooAbsReal& _x,
41  RooAbsReal& _mean,
42  Bool_t noRounding) :
43  RooAbsPdf(name,title),
44  x("x","x",this,_x),
45  mean("mean","mean",this,_mean),
46  _noRounding(noRounding),
47  _protectNegative(false)
48 {
49 }
50 
51 
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// Copy constructor
55 
56  RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :
57  RooAbsPdf(other,name),
58  x("x",this,other.x),
59  mean("mean",this,other.mean),
60  _noRounding(other._noRounding),
61  _protectNegative(other._protectNegative)
62 {
63 }
64 
65 
66 
67 
68 ////////////////////////////////////////////////////////////////////////////////
69 /// Implementation in terms of the TMath Poisson function
70 
72 {
73  Double_t k = _noRounding ? x : floor(x);
74  if(_protectNegative && mean<0)
75  return 1e-3;
76  return TMath::Poisson(k,mean) ;
77 }
78 
79 
80 
81 ////////////////////////////////////////////////////////////////////////////////
82 /// calculate and return the negative log-likelihood of the Poisson
83 
85 {
86  return RooAbsPdf::getLogVal(s) ;
87 // Double_t prob = getVal(s) ;
88 // return prob ;
89 
90  // Make inputs to naming conventions of RooAbsPdf::extendedTerm
91  Double_t expected=mean ;
92  Double_t observed=x ;
93 
94  // Explicitly handle case Nobs=Nexp=0
95  if (fabs(expected)<1e-10 && fabs(observed)<1e-10) {
96  return 0 ;
97  }
98 
99  // Explicitly handle case Nexp=0
100  if (fabs(observed)<1e-10) {
101  return -1*expected;
102  }
103 
104  // Michaels code for log(poisson) in RooAbsPdf::extendedTer with an approximated log(observed!) term
105  Double_t extra=0;
106  if(observed<1000000) {
107  extra = -observed*log(expected)+expected+TMath::LnGamma(observed+1.);
108  } else {
109  //if many observed events, use Gauss approximation
110  Double_t sigma_square=expected;
111  Double_t diff=observed-expected;
112  extra=-log(sigma_square)/2 + (diff*diff)/(2*sigma_square);
113  }
114 
115 // if (fabs(extra)>100 || log(prob)>100) {
116 // cout << "RooPoisson::getLogVal(" << GetName() << ") mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
117 // }
118 
119 // if (fabs(extra+log(prob))>1) {
120 // cout << "RooPoisson::getLogVal(" << GetName() << ") WARNING mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
121 // }
122 
123  //return log(prob);
124  return -extra-analyticalIntegral(1,0) ; //log(prob);
125 
126 }
127 
128 
129 ////////////////////////////////////////////////////////////////////////////////
130 
131 Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
132 {
133  if (matchArgs(allVars,analVars,x)) return 1 ;
134  if (matchArgs(allVars, analVars, mean)) return 2;
135  return 0 ;
136 }
137 
138 
139 
140 ////////////////////////////////////////////////////////////////////////////////
141 
142 Double_t RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const
143 {
144  R__ASSERT(code == 1 || code == 2) ;
145 
146  if(_protectNegative && mean<0)
147  return exp(-2*mean); // make it fall quickly
148 
149  if (code == 1) {
150  // Implement integral over x as summation. Add special handling in case
151  // range boundaries are not on integer values of x
152  Double_t xmin = x.min(rangeName) ;
153  Double_t xmax = x.max(rangeName) ;
154 
155  // Protect against negative lower boundaries
156  if (xmin<0) xmin=0 ;
157 
158  Int_t ixmin = Int_t (xmin) ;
159  Int_t ixmax = Int_t (xmax)+1 ;
160 
161  Double_t fracLoBin = 1-(xmin-ixmin) ;
162  Double_t fracHiBin = 1-(ixmax-xmax) ;
163 
164  if (!x.hasMax()) {
165  if (xmin<1e-6) {
166  return 1 ;
167  } else {
168 
169  // Return 1 minus integral from 0 to x.min()
170 
171  if(ixmin == 0){ // first bin
172  return TMath::Poisson(0, mean)*(xmin-0);
173  }
174  Double_t sum(0) ;
175  sum += TMath::Poisson(0,mean)*fracLoBin ;
177  sum += TMath::Poisson(ixmin-1,mean)*fracHiBin ;
178  return 1-sum ;
179  }
180  }
181 
182  if(ixmin == ixmax-1){ // first bin
183  return TMath::Poisson(ixmin, mean)*(xmax-xmin);
184  }
185 
186  Double_t sum(0) ;
187  sum += TMath::Poisson(ixmin,mean)*fracLoBin ;
188  if (RooNumber::isInfinite(xmax)){
189  sum+= 1.-ROOT::Math::poisson_cdf(ixmin,mean) ;
190  } else {
191  sum+= ROOT::Math::poisson_cdf(ixmax-2, mean) - ROOT::Math::poisson_cdf(ixmin,mean) ;
192  sum += TMath::Poisson(ixmax-1,mean)*fracHiBin ;
193  }
194 
195  return sum ;
196 
197  } else if(code == 2) {
198 
199  // the integral with respect to the mean is the integral of a gamma distribution
200  Double_t mean_min = mean.min(rangeName);
201  Double_t mean_max = mean.max(rangeName);
202 
203  Double_t ix;
204  if(_noRounding) ix = x + 1;
205  else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
206 
207  return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
208  }
209 
210  return 0;
211 
212 }
213 
214 
215 
216 
217 
218 
219 
220 ////////////////////////////////////////////////////////////////////////////////
221 /// Advertise internal generator in x
222 
223 Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
224 {
225  if (matchArgs(directVars,generateVars,x)) return 1 ;
226  return 0 ;
227 }
228 
229 
230 
231 ////////////////////////////////////////////////////////////////////////////////
232 /// Implement internal generator using TRandom::Poisson
233 
235 {
236  R__ASSERT(code==1) ;
237  Double_t xgen ;
238  while(1) {
240  if (xgen<=x.max() && xgen>=x.min()) {
241  x = xgen ;
242  break;
243  }
244  }
245  return;
246 }
247 
248 
Poisson pdf.
Definition: RooPoisson.h:19
static long int sum(long int i)
Definition: Factory.cxx:1786
float xmin
Definition: THbookFile.cxx:93
Double_t Floor(Double_t x)
Definition: TMath.h:473
RooRealProxy x
Definition: RooPoisson.h:40
Bool_t _noRounding
Definition: RooPoisson.h:42
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
void generateEvent(Int_t code)
Implement internal generator using TRandom::Poisson.
Definition: RooPoisson.cxx:234
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
Definition: RooPoisson.cxx:131
#define R__ASSERT(e)
Definition: TError.h:98
RooRealProxy mean
Definition: RooPoisson.h:41
double poisson_cdf(unsigned int n, double mu)
Cumulative distribution function of the Poisson distribution Lower tail of the integral of the poisso...
virtual Double_t getLogVal(const RooArgSet *set=0) const
Return the log of the current value with given normalization An error message is printed if the argum...
Definition: RooAbsPdf.cxx:606
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
STL namespace.
Double_t x[n]
Definition: legend1.C:17
Double_t getLogVal(const RooArgSet *set=0) const
calculate and return the negative log-likelihood of the Poisson
Definition: RooPoisson.cxx:84
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
Double_t evaluate() const
Implementation in terms of the TMath Poisson function.
Definition: RooPoisson.cxx:71
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Definition: RooPoisson.cxx:142
double floor(double)
float xmax
Definition: THbookFile.cxx:93
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par) The Poisson PDF is implemented by means of Eule...
Definition: TMath.cxx:571
#define ClassImp(name)
Definition: Rtypes.h:279
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
Bool_t _protectNegative
Definition: RooPoisson.h:43
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
Bool_t hasMax(const char *rname=0) const
Definition: RooRealProxy.h:59
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Advertise internal generator in x.
Definition: RooPoisson.cxx:223
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:489
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:362
double exp(double)
char name[80]
Definition: TGX11.cxx:109
double log(double)