Logo ROOT  
Reference Guide
PdfFuncMathCore.h
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/**
14
15Probability density functions, cumulative distribution functions
16and their inverses (quantiles) for various statistical distributions (continuous and discrete).
17Whenever possible the conventions followed are those of the
18CRC Concise Encyclopedia of Mathematics, Second Edition
19(or <A HREF="http://mathworld.wolfram.com/">Mathworld</A>).
20By convention the distributions are centered around 0, so for
21example in the case of a Gaussian there is no parameter mu. The
22user must calculate the shift themselves if they wish.
23
24MathCore provides the majority of the probability density functions, of the
25cumulative distributions and of the quantiles (inverses of the cumulatives).
26Additional distributions are also provided by the
27<A HREF="../../MathMore/html/group__StatFunc.html">MathMore</A> library.
28
29
30@defgroup StatFunc Statistical functions
31
32@ingroup MathCore
33@ingroup MathMore
34
35*/
36
37#ifndef ROOT_Math_PdfFuncMathCore
38#define ROOT_Math_PdfFuncMathCore
39
40#include "Math/Math.h"
42
43#include <limits>
44
45namespace ROOT {
46namespace Math {
47
48
49
50 /** @defgroup PdfFunc Probability Density Functions (PDF)
51 * @ingroup StatFunc
52 * Probability density functions of various statistical distributions
53 * (continuous and discrete).
54 * The probability density function returns the probability that
55 * the variate has the value x.
56 * In statistics the PDF is also called the frequency function.
57 *
58 *
59 */
60
61 /** @name Probability Density Functions from MathCore
62 * Additional PDF's are provided in the MathMore library
63 * (see PDF functions from MathMore)
64 */
65
66 //@{
67
68 /**
69
70 Probability density function of the beta distribution.
71
72 \f[ p(x) = \frac{\Gamma (a + b) } {\Gamma(a)\Gamma(b) } x ^{a-1} (1 - x)^{b-1} \f]
73
74 for \f$0 \leq x \leq 1 \f$. For detailed description see
75 <A HREF="http://mathworld.wolfram.com/BetaDistribution.html">
76 Mathworld</A>.
77
78 @ingroup PdfFunc
79
80 */
81
82 inline double beta_pdf(double x, double a, double b) {
83 // Inlined to enable clad-auto-derivation for this function.
84
85 if (x < 0 || x > 1.0) return 0;
86 if (x == 0 ) {
87 // need this work Windows
88 if (a < 1) return std::numeric_limits<double>::infinity();
89 else if (a > 1) return 0;
90 else if ( a == 1) return b; // to avoid a nan from log(0)*0
91 }
92 if (x == 1 ) {
93 // need this work Windows
94 if (b < 1) return std::numeric_limits<double>::infinity();
95 else if (b > 1) return 0;
96 else if ( b == 1) return a; // to avoid a nan from log(0)*0
97 }
99 std::log(x) * (a -1.) + ROOT::Math::log1p(-x ) * (b - 1.) );
100 }
101
102
103
104 /**
105
106 Probability density function of the binomial distribution.
107
108 \f[ p(k) = \frac{n!}{k! (n-k)!} p^k (1-p)^{n-k} \f]
109
110 for \f$ 0 \leq k \leq n \f$. For detailed description see
111 <A HREF="http://mathworld.wolfram.com/BinomialDistribution.html">
112 Mathworld</A>.
113
114 @ingroup PdfFunc
115
116 */
117
118 inline double binomial_pdf(unsigned int k, double p, unsigned int n) {
119 // Inlined to enable clad-auto-derivation for this function.
120 if (k > n)
121 return 0.0;
122
123 double coeff = ROOT::Math::lgamma(n+1) - ROOT::Math::lgamma(k+1) - ROOT::Math::lgamma(n-k+1);
124 return std::exp(coeff + k * std::log(p) + (n - k) * ROOT::Math::log1p(-p));
125 }
126
127
128
129 /**
130
131 Probability density function of the negative binomial distribution.
132
133 \f[ p(k) = \frac{(k+n-1)!}{k! (n-1)!} p^{n} (1-p)^{k} \f]
134
135 For detailed description see
136 <A HREF="http://mathworld.wolfram.com/NegativeBinomialDistribution.html">
137 Mathworld</A> (where \f$k \to x\f$ and \f$n \to r\f$).
138 The distribution in <A HREF="http://en.wikipedia.org/wiki/Negative_binomial_distribution">
139 Wikipedia</A> is defined with a \f$p\f$ corresponding to \f$1-p\f$ in this case.
140
141
142 @ingroup PdfFunc
143
144 */
145
146 inline double negative_binomial_pdf(unsigned int k, double p, double n) {
147 // Inlined to enable clad-auto-derivation for this function.
148
149 // implement in term of gamma function
150
151 if (n < 0) return 0.0;
152 if (p < 0 || p > 1.0) return 0.0;
153
154 double coeff = ROOT::Math::lgamma(k+n) - ROOT::Math::lgamma(k+1.0) - ROOT::Math::lgamma(n);
155 return std::exp(coeff + n * std::log(p) + double(k) * ROOT::Math::log1p(-p));
156
157 }
158
159
160
161
162 /**
163
164 Probability density function of Breit-Wigner distribution, which is similar, just
165 a different definition of the parameters, to the Cauchy distribution
166 (see #cauchy_pdf )
167
168 \f[ p(x) = \frac{1}{\pi} \frac{\frac{1}{2} \Gamma}{x^2 + (\frac{1}{2} \Gamma)^2} \f]
169
170
171 @ingroup PdfFunc
172
173 */
174
175 inline double breitwigner_pdf(double x, double gamma, double x0 = 0) {
176 // Inlined to enable clad-auto-derivation for this function.
177 double gammahalf = gamma/2.0;
178 return gammahalf/(M_PI * ((x-x0)*(x-x0) + gammahalf*gammahalf));
179 }
180
181
182
183
184 /**
185
186 Probability density function of the Cauchy distribution which is also
187 called Lorentzian distribution.
188
189
190 \f[ p(x) = \frac{1}{\pi} \frac{ b }{ (x-m)^2 + b^2} \f]
191
192 For detailed description see
193 <A HREF="http://mathworld.wolfram.com/CauchyDistribution.html">
194 Mathworld</A>. It is also related to the #breitwigner_pdf which
195 will call the same implementation.
196
197 @ingroup PdfFunc
198
199 */
200
201 inline double cauchy_pdf(double x, double b = 1, double x0 = 0) {
202
203 return b/(M_PI * ((x-x0)*(x-x0) + b*b));
204
205 }
206
207
208
209
210 /**
211
212 Probability density function of the \f$\chi^2\f$ distribution with \f$r\f$
213 degrees of freedom.
214
215 \f[ p_r(x) = \frac{1}{\Gamma(r/2) 2^{r/2}} x^{r/2-1} e^{-x/2} \f]
216
217 for \f$x \geq 0\f$. For detailed description see
218 <A HREF="http://mathworld.wolfram.com/Chi-SquaredDistribution.html">
219 Mathworld</A>.
220
221 @ingroup PdfFunc
222
223 */
224
225 inline double chisquared_pdf(double x, double r, double x0 = 0) {
226 // Inlined to enable clad-auto-derivation for this function.
227
228 if ((x-x0) < 0) {
229 return 0.0;
230 }
231 double a = r/2 -1.;
232 // let return inf for case x = x0 and treat special case of r = 2 otherwise will return nan
233 if (x == x0 && a == 0) return 0.5;
234
235 return std::exp ((r/2 - 1) * std::log((x-x0)/2) - (x-x0)/2 - ROOT::Math::lgamma(r/2))/2;
236
237 }
238
239
240 /**
241
242 Crystal ball function
243
244 See the definition at
245 <A HREF="http://en.wikipedia.org/wiki/Crystal_Ball_function">
246 Wikipedia</A>.
247
248 It is not really a pdf since it is not normalized
249
250 @ingroup PdfFunc
251
252 */
253
254 inline double crystalball_function(double x, double alpha, double n, double sigma, double mean = 0) {
255 // Inlined to enable clad-auto-derivation for this function.
256
257 // evaluate the crystal ball function
258 if (sigma < 0.) return 0.;
259 double z = (x - mean)/sigma;
260 if (alpha < 0) z = -z;
261 double abs_alpha = std::abs(alpha);
262 // double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
263 // double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
264 // double N = 1./(sigma*(C+D));
265 if (z > - abs_alpha)
266 return std::exp(- 0.5 * z * z);
267 //double A = std::pow(n/abs_alpha,n) * std::exp(-0.5*abs_alpha*abs_alpha);
268 double nDivAlpha = n/abs_alpha;
269 double AA = std::exp(-0.5*abs_alpha*abs_alpha);
270 double B = nDivAlpha -abs_alpha;
271 double arg = nDivAlpha/(B-z);
272 return AA * std::pow(arg,n);
273 }
274
275 /**
276 pdf definition of the crystal_ball which is defined only for n > 1 otherwise integral is diverging
277 */
278 inline double crystalball_pdf(double x, double alpha, double n, double sigma, double mean = 0) {
279 // Inlined to enable clad-auto-derivation for this function.
280
281 // evaluation of the PDF ( is defined only for n >1)
282 if (sigma < 0.) return 0.;
283 if ( n <= 1) return std::numeric_limits<double>::quiet_NaN(); // pdf is not normalized for n <=1
284 double abs_alpha = std::abs(alpha);
285 double C = n/abs_alpha * 1./(n-1.) * std::exp(-alpha*alpha/2.);
286 double D = std::sqrt(M_PI/2.)*(1.+ROOT::Math::erf(abs_alpha/std::sqrt(2.)));
287 double N = 1./(sigma*(C+D));
288 return N * crystalball_function(x,alpha,n,sigma,mean);
289 }
290
291 /**
292
293 Probability density function of the exponential distribution.
294
295 \f[ p(x) = \lambda e^{-\lambda x} \f]
296
297 for x>0. For detailed description see
298 <A HREF="http://mathworld.wolfram.com/ExponentialDistribution.html">
299 Mathworld</A>.
300
301
302 @ingroup PdfFunc
303
304 */
305
306 inline double exponential_pdf(double x, double lambda, double x0 = 0) {
307 // Inlined to enable clad-auto-derivation for this function.
308
309 if ((x-x0) < 0)
310 return 0.0;
311 return lambda * std::exp (-lambda * (x-x0));
312 }
313
314
315
316
317 /**
318
319 Probability density function of the F-distribution.
320
321 \f[ p_{n,m}(x) = \frac{\Gamma(\frac{n+m}{2})}{\Gamma(\frac{n}{2}) \Gamma(\frac{m}{2})} n^{n/2} m^{m/2} x^{n/2 -1} (m+nx)^{-(n+m)/2} \f]
322
323 for x>=0. For detailed description see
324 <A HREF="http://mathworld.wolfram.com/F-Distribution.html">
325 Mathworld</A>.
326
327 @ingroup PdfFunc
328
329 */
330
331
332 inline double fdistribution_pdf(double x, double n, double m, double x0 = 0) {
333 // Inlined to enable clad-auto-derivation for this function.
334
335 // function is defined only for both n and m > 0
336 if (n < 0 || m < 0)
337 return std::numeric_limits<double>::quiet_NaN();
338 if ((x-x0) < 0)
339 return 0.0;
340
341 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)
342 + (n/2 -1) * std::log(x-x0) - ((n+m)/2) * std::log(m + n*(x-x0)) );
343
344 }
345
346
347
348
349 /**
350
351 Probability density function of the gamma distribution.
352
353 \f[ p(x) = {1 \over \Gamma(\alpha) \theta^{\alpha}} x^{\alpha-1} e^{-x/\theta} \f]
354
355 for x>0. For detailed description see
356 <A HREF="http://mathworld.wolfram.com/GammaDistribution.html">
357 Mathworld</A>.
358
359 @ingroup PdfFunc
360
361 */
362
363 inline double gamma_pdf(double x, double alpha, double theta, double x0 = 0) {
364 // Inlined to enable clad-auto-derivation for this function.
365
366 if ((x-x0) < 0) {
367 return 0.0;
368 } else if ((x-x0) == 0) {
369
370 if (alpha == 1) {
371 return 1.0/theta;
372 } else {
373 return 0.0;
374 }
375
376 } else if (alpha == 1) {
377 return std::exp(-(x-x0)/theta)/theta;
378 } else {
379 return std::exp((alpha - 1) * std::log((x-x0)/theta) - (x-x0)/theta - ROOT::Math::lgamma(alpha))/theta;
380 }
381
382 }
383
384
385
386
387 /**
388
389 Probability density function of the normal (Gaussian) distribution.
390
391 \f[ p(x) = {1 \over \sqrt{2 \pi \sigma^2}} e^{-x^2 / 2\sigma^2} \f]
392
393 For detailed description see
394 <A HREF="http://mathworld.wolfram.com/NormalDistribution.html">
395 Mathworld</A>. It can also be evaluated using #normal_pdf which will
396 call the same implementation.
397
398 @ingroup PdfFunc
399
400 */
401
402 inline double gaussian_pdf(double x, double sigma = 1, double x0 = 0) {
403
404 double tmp = (x-x0)/sigma;
405 return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
406 }
407
408 /**
409
410 Probability density function of the bi-dimensional (Gaussian) distribution.
411
412 \f[ p(x) = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp (-(x^2/\sigma_x^2 + y^2/\sigma_y^2 - 2 \rho x y/(\sigma_x\sigma_y))/2(1-\rho^2)) \f]
413
414 For detailed description see
415 <A HREF="http://mathworld.wolfram.com/BivariateNormalDistribution.html">
416 Mathworld</A>. It can also be evaluated using #normal_pdf which will
417 call the same implementation.
418
419 @param rho correlation , must be between -1,1
420
421 @ingroup PdfFunc
422
423 */
424
425 inline double bigaussian_pdf(double x, double y, double sigmax = 1, double sigmay = 1, double rho = 0, double x0 = 0, double y0 = 0) {
426 double u = (x-x0)/sigmax;
427 double v = (y-y0)/sigmay;
428 double c = 1. - rho*rho;
429 double z = u*u - 2.*rho*u*v + v*v;
430 return 1./(2 * M_PI * sigmax * sigmay * std::sqrt(c) ) * std::exp(- z / (2. * c) );
431 }
432
433 /**
434
435 Probability density function of the Landau distribution:
436 \f[ p(x) = \frac{1}{\xi} \phi (\lambda) \f]
437 with
438 \f[ \phi(\lambda) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{\lambda s + s \log{s}} ds\f]
439 where \f$\lambda = (x-x_0)/\xi\f$. For a detailed description see
440 K.S. K&ouml;lbig and B. Schorr, A program package for the Landau distribution,
441 <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
442 <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
443 The same algorithms as in
444 <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g110/top.html">
445 CERNLIB</A> (DENLAN) is used
446
447 @param x The argument \f$x\f$
448 @param xi The width parameter \f$\xi\f$
449 @param x0 The location parameter \f$x_0\f$
450
451 @ingroup PdfFunc
452
453 */
454
455 double landau_pdf(double x, double xi = 1, double x0 = 0);
456
457
458
459 /**
460
461 Probability density function of the lognormal distribution.
462
463 \f[ p(x) = {1 \over x \sqrt{2 \pi s^2} } e^{-(\ln{x} - m)^2/2 s^2} \f]
464
465 for x>0. For detailed description see
466 <A HREF="http://mathworld.wolfram.com/LogNormalDistribution.html">
467 Mathworld</A>.
468 @param s scale parameter (not the sigma of the distribution which is not even defined)
469 @param x0 location parameter, corresponds approximately to the most probable value. For x0 = 0, sigma = 1, the x_mpv = -0.22278
470
471 @ingroup PdfFunc
472
473 */
474
475 inline double lognormal_pdf(double x, double m, double s, double x0 = 0) {
476 // Inlined to enable clad-auto-derivation for this function.
477 if ((x-x0) <= 0)
478 return 0.0;
479 double tmp = (std::log((x-x0)) - m)/s;
480 return 1.0 / ((x-x0) * std::fabs(s) * std::sqrt(2 * M_PI)) * std::exp(-(tmp * tmp) /2);
481 }
482
483
484
485
486 /**
487
488 Probability density function of the normal (Gaussian) distribution.
489
490 \f[ p(x) = {1 \over \sqrt{2 \pi \sigma^2}} e^{-x^2 / 2\sigma^2} \f]
491
492 For detailed description see
493 <A HREF="http://mathworld.wolfram.com/NormalDistribution.html">
494 Mathworld</A>. It can also be evaluated using #gaussian_pdf which will call the same
495 implementation.
496
497 @ingroup PdfFunc
498
499 */
500
501 inline double normal_pdf(double x, double sigma =1, double x0 = 0) {
502 // Inlined to enable clad-auto-derivation for this function.
503
504 double tmp = (x-x0)/sigma;
505 return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
506
507 }
508
509
510 /**
511
512 Probability density function of the Poisson distribution.
513
514 \f[ p(n) = \frac{\mu^n}{n!} e^{- \mu} \f]
515
516 For detailed description see
517 <A HREF="http://mathworld.wolfram.com/PoissonDistribution.html">
518 Mathworld</A>.
519
520 @ingroup PdfFunc
521
522 */
523
524 inline double poisson_pdf(unsigned int n, double mu) {
525 // Inlined to enable clad-auto-derivation for this function.
526
527 if (n > 0 && mu >= 0)
528 return std::exp (n*std::log(mu) - ROOT::Math::lgamma(n+1) - mu);
529
530 // when n = 0 and mu = 0, 1 is returned
531 if (mu >= 0)
532 return std::exp(-mu);
533
534 // return a nan for mu < 0 since it does not make sense
535 return std::numeric_limits<double>::quiet_NaN();
536 }
537
538
539
540
541 /**
542
543 Probability density function of Student's t-distribution.
544
545 \f[ p_{r}(x) = \frac{\Gamma(\frac{r+1}{2})}{\sqrt{r \pi}\Gamma(\frac{r}{2})} \left( 1+\frac{x^2}{r}\right)^{-(r+1)/2} \f]
546
547 for \f$k \geq 0\f$. For detailed description see
548 <A HREF="http://mathworld.wolfram.com/Studentst-Distribution.html">
549 Mathworld</A>.
550
551 @ingroup PdfFunc
552
553 */
554
555 inline double tdistribution_pdf(double x, double r, double x0 = 0) {
556 // Inlined to enable clad-auto-derivation for this function.
557
558 return (std::exp (ROOT::Math::lgamma((r + 1.0)/2.0) - ROOT::Math::lgamma(r/2.0)) / std::sqrt (M_PI * r))
559 * std::pow ((1.0 + (x-x0)*(x-x0)/r), -(r + 1.0)/2.0);
560
561 }
562
563
564
565
566 /**
567
568 Probability density function of the uniform (flat) distribution.
569
570 \f[ p(x) = {1 \over (b-a)} \f]
571
572 if \f$a \leq x<b\f$ and 0 otherwise. For detailed description see
573 <A HREF="http://mathworld.wolfram.com/UniformDistribution.html">
574 Mathworld</A>.
575
576 @ingroup PdfFunc
577
578 */
579
580 inline double uniform_pdf(double x, double a, double b, double x0 = 0) {
581 // Inlined to enable clad-auto-derivation for this function.
582
583 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! when a=b
584
585 if ((x-x0) < b && (x-x0) >= a)
586 return 1.0/(b - a);
587 return 0.0;
588
589 }
590
591
592
593} // namespace Math
594} // namespace ROOT
595
596
597
598#endif // ROOT_Math_PdfFunc
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define M_PI
Definition: Rotated.cxx:105
#define N
double uniform_pdf(double x, double a, double b, double x0=0)
Probability density function of the uniform (flat) 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.
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
double binomial_pdf(unsigned int k, double p, unsigned int n)
Probability density function of the binomial distribution.
double fdistribution_pdf(double x, double n, double m, double x0=0)
Probability density function of the F-distribution.
double negative_binomial_pdf(unsigned int k, double p, double n)
Probability density function of the negative binomial distribution.
double gamma_pdf(double x, double alpha, double theta, double x0=0)
Probability density function of the gamma distribution.
double landau_pdf(double x, double xi=1, double x0=0)
Probability density function of the Landau distribution:
double exponential_pdf(double x, double lambda, double x0=0)
Probability density function of the exponential distribution.
double gaussian_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
double chisquared_pdf(double x, double r, double x0=0)
Probability density function of the distribution with degrees of freedom.
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 crystalball_function(double x, double alpha, double n, double sigma, double mean=0)
Crystal ball function.
double crystalball_pdf(double x, double alpha, double n, double sigma, double mean=0)
pdf definition of the crystal_ball which is defined only for n > 1 otherwise integral is diverging
double beta_pdf(double x, double a, double b)
Probability density function of the beta distribution.
double poisson_pdf(unsigned int n, double mu)
Probability density function of the Poisson distribution.
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 tdistribution_pdf(double x, double r, double x0=0)
Probability density function of Student's t-distribution.
double lgamma(double x)
Calculates the logarithm of the gamma function.
double erf(double x)
Error function encountered in integrating the normal distribution.
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
static double B[]
static double C[]
double gamma(double x)
double log1p(double x)
declarations for functions which are not implemented by some compilers
Definition: Math.h:98
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
static constexpr double s
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12