ROOT   6.10/09 Reference Guide
PdfFuncMathCore.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: Andras Zsenei & Lorenzo Moneta 06/2005
3
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10
11
12
13 #include "Math/Math.h"
14 #include "Math/SpecFuncMathCore.h"
15 #include <limits>
16
17
18 namespace ROOT {
19 namespace Math {
20
21
22  double beta_pdf(double x, double a, double b) {
23  if (x < 0 || x > 1.0) return 0;
24  if (x == 0 ) {
25  // need this wor Windows
26  if (a < 1) return std::numeric_limits<double>::infinity();
27  else if (a > 1) return 0;
28  else if ( a == 1) return b; // to avoid a nan from log(0)*0
29  }
30  if (x == 1 ) {
31  // need this wor Windows
32  if (b < 1) return std::numeric_limits<double>::infinity();
33  else if (b > 1) return 0;
34  else if ( b == 1) return a; // to avoid a nan from log(0)*0
35  }
37  std::log(x) * (a -1.) + ROOT::Math::log1p(-x ) * (b - 1.) );
38  }
39
40  double binomial_pdf(unsigned int k, double p, unsigned int n) {
41
42  if (k > n) {
43  return 0.0;
44  } else {
45
46  double coeff = ROOT::Math::lgamma(n+1) - ROOT::Math::lgamma(k+1) - ROOT::Math::lgamma(n-k+1);
47  return std::exp(coeff + k * std::log(p) + (n - k) * ROOT::Math::log1p(-p));
48  }
49  }
50
51  double negative_binomial_pdf(unsigned int k, double p, double n) {
52  // impelment in term of gamma function
53
54  if (n < 0) return 0.0;
55  if (p < 0 || p > 1.0) return 0.0;
56
57  double coeff = ROOT::Math::lgamma(k+n) - ROOT::Math::lgamma(k+1.0) - ROOT::Math::lgamma(n);
58  return std::exp(coeff + n * std::log(p) + double(k) * ROOT::Math::log1p(-p));
59
60  }
61
62
63  double breitwigner_pdf(double x, double gamma, double x0) {
64
65  double gammahalf = gamma/2.0;
66  return gammahalf/(M_PI * ((x-x0)*(x-x0) + gammahalf*gammahalf));
67
68
69  }
70
71
72
73  double cauchy_pdf(double x, double b, double x0) {
74
75  return b/(M_PI * ((x-x0)*(x-x0) + b*b));
76
77  }
78
79
80
81  double chisquared_pdf(double x, double r, double x0) {
82
83  if ((x-x0) < 0) {
84  return 0.0;
85  }
86  double a = r/2 -1.;
87  // let return inf for case x = x0 and treat special case of r = 2 otherwise will return nan
88  if (x == x0 && a == 0) return 0.5;
89
90  return std::exp ((r/2 - 1) * std::log((x-x0)/2) - (x-x0)/2 - ROOT::Math::lgamma(r/2))/2;
91
92  }
93
94  double crystalball_function(double x, double alpha, double n, double sigma, double mean) {
95  // evaluate the crystal ball function
96  if (sigma < 0.) return 0.;
97  double z = (x - mean)/sigma;
98  if (alpha < 0) z = -z;
99  double abs_alpha = std::abs(alpha);
100  // double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
101  // double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
102  // double N = 1./(sigma*(C+D));
103  if (z > - abs_alpha)
104  return std::exp(- 0.5 * z * z);
105  else {
106  //double A = std::pow(n/abs_alpha,n) * std::exp(-0.5*abs_alpha*abs_alpha);
107  double nDivAlpha = n/abs_alpha;
108  double AA = std::exp(-0.5*abs_alpha*abs_alpha);
109  double B = nDivAlpha -abs_alpha;
110  double arg = nDivAlpha/(B-z);
111  return AA * std::pow(arg,n);
112  }
113  }
114  double crystalball_pdf(double x, double alpha, double n, double sigma, double mean) {
115  // evaluation of the PDF ( is defined only for n >1)
116  if (sigma < 0.) return 0.;
117  if ( n <= 1) return std::numeric_limits<double>::quiet_NaN(); // pdf is not normalized for n <=1
118  double abs_alpha = std::abs(alpha);
119  double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
120  double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
121  double N = 1./(sigma*(C+D));
122  return N * crystalball_function(x,alpha,n,sigma,mean);
123  }
124
125
126  double exponential_pdf(double x, double lambda, double x0) {
127
128  if ((x-x0) < 0) {
129  return 0.0;
130  } else {
131  return lambda * std::exp (-lambda * (x-x0));
132  }
133
134  }
135
136
137
138  double fdistribution_pdf(double x, double n, double m, double x0) {
139
140  // function is defined only for both n and m > 0
141  if (n < 0 || m < 0)
142  return std::numeric_limits<double>::quiet_NaN();
143  if ((x-x0) < 0)
144  return 0.0;
145
146  return std::exp((n/2) * std::log(n) + (m/2) * std::log(m) + ROOT::Math::lgamma((n+m)/2) - ROOT::Math::lgamma(n/2) - ROOT::Math::lgamma(m/2)
147  + (n/2 -1) * std::log(x-x0) - ((n+m)/2) * std::log(m + n*(x-x0)) );
148
149  }
150
151
152
153  double gamma_pdf(double x, double alpha, double theta, double x0) {
154
155  if ((x-x0) < 0) {
156  return 0.0;
157  } else if ((x-x0) == 0) {
158
159  if (alpha == 1) {
160  return 1.0/theta;
161  } else {
162  return 0.0;
163  }
164
165  } else if (alpha == 1) {
166  return std::exp(-(x-x0)/theta)/theta;
167  } else {
168  return std::exp((alpha - 1) * std::log((x-x0)/theta) - (x-x0)/theta - ROOT::Math::lgamma(alpha))/theta;
169  }
170
171  }
172
173
174
175  double gaussian_pdf(double x, double sigma, double x0) {
176
177  double tmp = (x-x0)/sigma;
178  return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
179  }
180
181  double bigaussian_pdf(double x, double y, double sigmax , double sigmay , double rho , double x0 , double y0 ) {
182
183  double u = (x-x0)/sigmax;
184  double v = (y-y0)/sigmay;
185  double c = 1. - rho*rho;
186  double z = u*u - 2.*rho*u*v + v*v;
187  return 1./(2 * M_PI * sigmax * sigmay * std::sqrt(c) ) * std::exp(- z / (2. * c) );
188  }
189
190
191  double landau_pdf(double x, double xi, double x0) {
192  // LANDAU pdf : algorithm from CERNLIB G110 denlan
193  // same algorithm is used in GSL
194
195  static double p1[5] = {0.4259894875,-0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253};
196  static double q1[5] = {1.0 ,-0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063};
197
198  static double p2[5] = {0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211};
199  static double q2[5] = {1.0 , 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714};
200
201  static double p3[5] = {0.1788544503, 0.09359161662,0.006325387654, 0.00006611667319,-0.000002031049101};
202  static double q3[5] = {1.0 , 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675};
203
204  static double p4[5] = {0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186};
205  static double q4[5] = {1.0 , 106.8615961, 337.6496214, 2016.712389, 1597.063511};
206
207  static double p5[5] = {1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910};
208  static double q5[5] = {1.0 , 156.9424537, 3745.310488, 9834.698876, 66924.28357};
209
210  static double p6[5] = {1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109};
211  static double q6[5] = {1.0 , 651.4101098, 56974.73333, 165917.4725, -2815759.939};
212
213  static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
214
215  static double a2[2] = {-1.845568670,-4.284640743};
216
217  if (xi <= 0) return 0;
218  double v = (x - x0)/xi;
219  double u, ue, us, denlan;
220  if (v < -5.5) {
221  u = std::exp(v+1.0);
222  if (u < 1e-10) return 0.0;
223  ue = std::exp(-1/u);
224  us = std::sqrt(u);
225  denlan = 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
226  } else if(v < -1) {
227  u = std::exp(-v-1);
228  denlan = std::exp(-u)*std::sqrt(u)*
229  (p1[0]+(p1[1]+(p1[2]+(p1[3]+p1[4]*v)*v)*v)*v)/
230  (q1[0]+(q1[1]+(q1[2]+(q1[3]+q1[4]*v)*v)*v)*v);
231  } else if(v < 1) {
232  denlan = (p2[0]+(p2[1]+(p2[2]+(p2[3]+p2[4]*v)*v)*v)*v)/
233  (q2[0]+(q2[1]+(q2[2]+(q2[3]+q2[4]*v)*v)*v)*v);
234  } else if(v < 5) {
235  denlan = (p3[0]+(p3[1]+(p3[2]+(p3[3]+p3[4]*v)*v)*v)*v)/
236  (q3[0]+(q3[1]+(q3[2]+(q3[3]+q3[4]*v)*v)*v)*v);
237  } else if(v < 12) {
238  u = 1/v;
239  denlan = u*u*(p4[0]+(p4[1]+(p4[2]+(p4[3]+p4[4]*u)*u)*u)*u)/
240  (q4[0]+(q4[1]+(q4[2]+(q4[3]+q4[4]*u)*u)*u)*u);
241  } else if(v < 50) {
242  u = 1/v;
243  denlan = u*u*(p5[0]+(p5[1]+(p5[2]+(p5[3]+p5[4]*u)*u)*u)*u)/
244  (q5[0]+(q5[1]+(q5[2]+(q5[3]+q5[4]*u)*u)*u)*u);
245  } else if(v < 300) {
246  u = 1/v;
247  denlan = u*u*(p6[0]+(p6[1]+(p6[2]+(p6[3]+p6[4]*u)*u)*u)*u)/
248  (q6[0]+(q6[1]+(q6[2]+(q6[3]+q6[4]*u)*u)*u)*u);
249  } else {
250  u = 1/(v-v*std::log(v)/(v+1));
251  denlan = u*u*(1+(a2[0]+a2[1]*u)*u);
252  }
253  return denlan/xi;
254
255  }
256
257
258
259
260
261  double lognormal_pdf(double x, double m, double s, double x0) {
262
263  if ((x-x0) <= 0) {
264  return 0.0;
265  } else {
266  double tmp = (std::log((x-x0)) - m)/s;
267  return 1.0 / ((x-x0) * std::fabs(s) * std::sqrt(2 * M_PI)) * std::exp(-(tmp * tmp) /2);
268  }
269
270  }
271
272
273
274  double normal_pdf(double x, double sigma, double x0) {
275
276  double tmp = (x-x0)/sigma;
277  return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
278
279  }
280
281
282
283  double poisson_pdf(unsigned int n, double mu) {
284
285  if (n > 0)
286  return std::exp (n*std::log(mu) - ROOT::Math::lgamma(n+1) - mu);
287  else {
288  // when n = 0 and mu = 0, 1 is returned
289  if (mu >= 0) return std::exp(-mu);
290  // return a nan for mu < 0 since it does not make sense
291  return std::log(mu);
292  }
293
294  }
295
296
297
298  double tdistribution_pdf(double x, double r, double x0) {
299
300  return (std::exp (ROOT::Math::lgamma((r + 1.0)/2.0) - ROOT::Math::lgamma(r/2.0)) / std::sqrt (M_PI * r))
301  * std::pow ((1.0 + (x-x0)*(x-x0)/r), -(r + 1.0)/2.0);
302
303  }
304
305
306
307  double uniform_pdf(double x, double a, double b, double x0) {
308
309  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! when a=b
310
311  if ((x-x0) < b && (x-x0) >= a) {
312  return 1.0/(b - a);
313  } else {
314  return 0.0;
315  }
316
317  }
318
319
320 } // namespace Math
321 } // namespace ROOT
322
323
324
325
326
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
static double B[]
double erf(double x)
Error function encountered in integrating the normal distribution.
double tdistribution_pdf(double x, double r, double x0=0)
Probability density function of Student&#39;s t-distribution.
static double p3(double t, double a, double b, double c, double d)
double chisquared_pdf(double x, double r, double x0=0)
Probability density function of the distribution with degrees of freedom.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double beta_pdf(double x, double a, double b)
Probability density function of the beta distribution.
double bigaussian_pdf(double x, double y, double sigmax=1, double sigmay=1, double rho=0, double x0=0, double y0=0)
Probability density function of the bi-dimensional (Gaussian) distribution.
#define N
TArc * a
Definition: textangle.C:12
double exponential_pdf(double x, double lambda, double x0=0)
Probability density function of the exponential distribution.
double crystalball_function(double x, double alpha, double n, double sigma, double x0=0)
Crystal ball function.
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution: with where .
static double p2(double t, double a, double b, double c)
double pow(double, double)
double gaussian_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
const Double_t sigma
double gamma_pdf(double x, double alpha, double theta, double x0=0)
Probability density function of the gamma distribution.
double breitwigner_pdf(double x, double gamma, double x0=0)
Probability density function of Breit-Wigner distribution, which is similar, just a different definit...
double negative_binomial_pdf(unsigned int k, double p, double n)
Probability density function of the negative binomial distribution.
double gamma(double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static double C[]
TRandom2 r(17)
#define M_PI
Definition: Rotated.cxx:105
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
SVector< double, 2 > v
Definition: Dict.h:5
double fdistribution_pdf(double x, double n, double m, double x0=0)
Probability density function of the F-distribution.
TMarker * m
Definition: textangle.C:8
static double p1(double t, double a, double b)
double log1p(double x)
declarations for functions which are not implemented by some compilers
Definition: Math.h:100
double cauchy_pdf(double x, double b=1, double x0=0)
Probability density function of the Cauchy distribution which is also called Lorentzian distribution...
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double uniform_pdf(double x, double a, double b, double x0=0)
Probability density function of the uniform (flat) distribution.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double crystalball_pdf(double x, double alpha, double n, double sigma, double x0=0)
pdf definition of the crystal_ball which is defined only for n > 1 otherwise integral is diverging ...
double lgamma(double x)
Calculates the logarithm of the gamma function.
double exp(double)
const Int_t n
Definition: legend1.C:16
double log(double)