Logo ROOT   6.10/09
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 /** \class RooPoisson
10  \ingroup Roofit
11 
12 Poisson pdf
13 **/
14 
15 #include <iostream>
16 
17 #include "RooPoisson.h"
18 #include "RooAbsReal.h"
19 #include "RooAbsCategory.h"
20 
21 #include "RooRandom.h"
22 #include "RooMath.h"
23 #include "TMath.h"
24 #include "Math/ProbFuncMathCore.h"
25 
26 #include "TError.h"
27 
28 using namespace std;
29 
31 
32 ////////////////////////////////////////////////////////////////////////////////
33 /// Constructor
34 
35 RooPoisson::RooPoisson(const char *name, const char *title,
36  RooAbsReal& _x,
37  RooAbsReal& _mean,
38  Bool_t noRounding) :
39  RooAbsPdf(name,title),
40  x("x","x",this,_x),
41  mean("mean","mean",this,_mean),
42  _noRounding(noRounding),
43  _protectNegative(false)
44 {
45 }
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 /// Copy constructor
49 
50  RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :
51  RooAbsPdf(other,name),
52  x("x",this,other.x),
53  mean("mean",this,other.mean),
54  _noRounding(other._noRounding),
55  _protectNegative(other._protectNegative)
56 {
57 }
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// Implementation in terms of the TMath Poisson function
61 
63 {
64  Double_t k = _noRounding ? x : floor(x);
65  if(_protectNegative && mean<0)
66  return 1e-3;
67  return TMath::Poisson(k,mean) ;
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// calculate and return the negative log-likelihood of the Poisson
72 
74 {
75  return RooAbsPdf::getLogVal(s) ;
76 // Double_t prob = getVal(s) ;
77 // return prob ;
78 
79  // Make inputs to naming conventions of RooAbsPdf::extendedTerm
80  Double_t expected=mean ;
81  Double_t observed=x ;
82 
83  // Explicitly handle case Nobs=Nexp=0
84  if (fabs(expected)<1e-10 && fabs(observed)<1e-10) {
85  return 0 ;
86  }
87 
88  // Explicitly handle case Nexp=0
89  if (fabs(observed)<1e-10) {
90  return -1*expected;
91  }
92 
93  // Michaels code for log(poisson) in RooAbsPdf::extendedTer with an approximated log(observed!) term
94  Double_t extra=0;
95  if(observed<1000000) {
96  extra = -observed*log(expected)+expected+TMath::LnGamma(observed+1.);
97  } else {
98  //if many observed events, use Gauss approximation
99  Double_t sigma_square=expected;
100  Double_t diff=observed-expected;
101  extra=-log(sigma_square)/2 + (diff*diff)/(2*sigma_square);
102  }
103 
104 // if (fabs(extra)>100 || log(prob)>100) {
105 // cout << "RooPoisson::getLogVal(" << GetName() << ") mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
106 // }
107 
108 // if (fabs(extra+log(prob))>1) {
109 // cout << "RooPoisson::getLogVal(" << GetName() << ") WARNING mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
110 // }
111 
112  //return log(prob);
113  return -extra-analyticalIntegral(1,0) ; //log(prob);
114 
115 }
116 
117 ////////////////////////////////////////////////////////////////////////////////
118 
119 Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
120 {
121  if (matchArgs(allVars,analVars,x)) return 1 ;
122  if (matchArgs(allVars, analVars, mean)) return 2;
123  return 0 ;
124 }
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 
128 Double_t RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const
129 {
130  R__ASSERT(code == 1 || code == 2) ;
131 
132  if(_protectNegative && mean<0)
133  return exp(-2*mean); // make it fall quickly
134 
135  if (code == 1) {
136  // Implement integral over x as summation. Add special handling in case
137  // range boundaries are not on integer values of x
138  Double_t xmin = x.min(rangeName) ;
139  Double_t xmax = x.max(rangeName) ;
140 
141  // Protect against negative lower boundaries
142  if (xmin<0) xmin=0 ;
143 
144  Int_t ixmin = Int_t (xmin) ;
145  Int_t ixmax = Int_t (xmax)+1 ;
146 
147  Double_t fracLoBin = 1-(xmin-ixmin) ;
148  Double_t fracHiBin = 1-(ixmax-xmax) ;
149 
150  if (!x.hasMax()) {
151  if (xmin<1e-6) {
152  return 1 ;
153  } else {
154 
155  // Return 1 minus integral from 0 to x.min()
156 
157  if(ixmin == 0){ // first bin
158  return TMath::Poisson(0, mean)*(xmin-0);
159  }
160  Double_t sum(0) ;
161  sum += TMath::Poisson(0,mean)*fracLoBin ;
163  sum += TMath::Poisson(ixmin-1,mean)*fracHiBin ;
164  return 1-sum ;
165  }
166  }
167 
168  if(ixmin == ixmax-1){ // first bin
169  return TMath::Poisson(ixmin, mean)*(xmax-xmin);
170  }
171 
172  Double_t sum(0) ;
173  sum += TMath::Poisson(ixmin,mean)*fracLoBin ;
174  if (RooNumber::isInfinite(xmax)){
175  sum+= 1.-ROOT::Math::poisson_cdf(ixmin,mean) ;
176  } else {
177  sum+= ROOT::Math::poisson_cdf(ixmax-2, mean) - ROOT::Math::poisson_cdf(ixmin,mean) ;
178  sum += TMath::Poisson(ixmax-1,mean)*fracHiBin ;
179  }
180 
181  return sum ;
182 
183  } else if(code == 2) {
184 
185  // the integral with respect to the mean is the integral of a gamma distribution
186  Double_t mean_min = mean.min(rangeName);
187  Double_t mean_max = mean.max(rangeName);
188 
189  Double_t ix;
190  if(_noRounding) ix = x + 1;
191  else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
192 
193  return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
194  }
195 
196  return 0;
197 
198 }
199 
200 ////////////////////////////////////////////////////////////////////////////////
201 /// Advertise internal generator in x
202 
203 Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
204 {
205  if (matchArgs(directVars,generateVars,x)) return 1 ;
206  return 0 ;
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////////
210 /// Implement internal generator using TRandom::Poisson
211 
213 {
214  R__ASSERT(code==1) ;
215  Double_t xgen ;
216  while(1) {
218  if (xgen<=x.max() && xgen>=x.min()) {
219  x = xgen ;
220  break;
221  }
222  }
223  return;
224 }
Poisson pdf.
Definition: RooPoisson.h:19
static long int sum(long int i)
Definition: Factory.cxx:2162
float xmin
Definition: THbookFile.cxx:93
Double_t Floor(Double_t x)
Definition: TMath.h:600
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:212
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:119
#define R__ASSERT(e)
Definition: TError.h:96
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:73
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:62
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:128
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:336
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:203
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)
double log(double)