Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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-x0)^2/\sigma_x^2 + (y-y0)^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 x x variable
420 @param y y variable
421 @param sigmax the stdev in x
422 @param sigmay the stdev in y
423 @param rho correlation, must be between -1,1
424 @param x0 the offset in x
425 @param y0 the offset in y
426
427 @ingroup PdfFunc
428
429 */
430
431 inline double bigaussian_pdf(double x, double y, double sigmax = 1, double sigmay = 1, double rho = 0, double x0 = 0, double y0 = 0) {
432 double u = (x-x0)/sigmax;
433 double v = (y-y0)/sigmay;
434 double c = 1. - rho*rho;
435 double z = u*u - 2.*rho*u*v + v*v;
436 return 1./(2 * M_PI * sigmax * sigmay * std::sqrt(c) ) * std::exp(- z / (2. * c) );
437 }
438
439 /**
440
441 Probability density function of the Landau distribution:
442 \f[ p(x) = \frac{1}{\xi} \phi (\lambda) \f]
443 with
444 \f[ \phi(\lambda) = \frac{1}{2 \pi i}\int_{c-i\infty}^{c+i\infty} e^{\lambda s + s \log{s}} ds\f]
445 where \f$\lambda = (x-x_0)/\xi\f$. For a detailed description see
446 K.S. K&ouml;lbig and B. Schorr, A program package for the Landau distribution,
447 <A HREF="http://dx.doi.org/10.1016/0010-4655(84)90085-7">Computer Phys. Comm. 31 (1984) 97-111</A>
448 <A HREF="http://dx.doi.org/10.1016/j.cpc.2008.03.002">[Erratum-ibid. 178 (2008) 972]</A>.
449 The same algorithms as in
450 <A HREF="https://cern-tex.web.cern.ch/cern-tex/shortwrupsdir/g110/top.html">
451 CERNLIB</A> (DENLAN) is used
452
453 @param x The argument \f$x\f$
454 @param xi The width parameter \f$\xi\f$
455 @param x0 The location parameter \f$x_0\f$
456
457 @ingroup PdfFunc
458
459 */
460
461 double landau_pdf(double x, double xi = 1, double x0 = 0);
462
463
464
465 /**
466
467 Probability density function of the lognormal distribution.
468
469 \f[ p(x) = {1 \over x \sqrt{2 \pi s^2} } e^{-(\ln{x} - m)^2/2 s^2} \f]
470
471 for x>0. For detailed description see
472 <A HREF="http://mathworld.wolfram.com/LogNormalDistribution.html">
473 Mathworld</A>.
474 @param x x variable
475 @param m M = 0 for lognormal
476 @param s scale parameter (not the sigma of the distribution which is not even defined)
477 @param x0 location parameter, corresponds approximately to the most probable value. For x0 = 0, sigma = 1, the x_mpv = -0.22278
478
479 @ingroup PdfFunc
480
481 */
482
483 inline double lognormal_pdf(double x, double m, double s, double x0 = 0) {
484 // Inlined to enable clad-auto-derivation for this function.
485 if ((x-x0) <= 0)
486 return 0.0;
487 double tmp = (std::log((x-x0)) - m)/s;
488 return 1.0 / ((x-x0) * std::fabs(s) * std::sqrt(2 * M_PI)) * std::exp(-(tmp * tmp) /2);
489 }
490
491
492
493
494 /**
495
496 Probability density function of the normal (Gaussian) distribution.
497
498 \f[ p(x) = {1 \over \sqrt{2 \pi \sigma^2}} e^{-x^2 / 2\sigma^2} \f]
499
500 For detailed description see
501 <A HREF="http://mathworld.wolfram.com/NormalDistribution.html">
502 Mathworld</A>. It can also be evaluated using #gaussian_pdf which will call the same
503 implementation.
504
505 @ingroup PdfFunc
506
507 */
508
509 inline double normal_pdf(double x, double sigma =1, double x0 = 0) {
510 // Inlined to enable clad-auto-derivation for this function.
511
512 double tmp = (x-x0)/sigma;
513 return (1.0/(std::sqrt(2 * M_PI) * std::fabs(sigma))) * std::exp(-tmp*tmp/2);
514
515 }
516
517
518 /**
519
520 Probability density function of the Poisson distribution.
521
522 \f[ p(n) = \frac{\mu^n}{n!} e^{- \mu} \f]
523
524 For detailed description see
525 <A HREF="http://mathworld.wolfram.com/PoissonDistribution.html">
526 Mathworld</A>.
527
528 @ingroup PdfFunc
529
530 */
531
532 inline double poisson_pdf(unsigned int n, double mu) {
533 // Inlined to enable clad-auto-derivation for this function.
534
535 if (n > 0 && mu >= 0)
536 return std::exp (n*std::log(mu) - ROOT::Math::lgamma(n+1) - mu);
537
538 // when n = 0 and mu = 0, 1 is returned
539 if (mu >= 0)
540 return std::exp(-mu);
541
542 // return a nan for mu < 0 since it does not make sense
543 return std::numeric_limits<double>::quiet_NaN();
544 }
545
546
547
548
549 /**
550
551 Probability density function of Student's t-distribution.
552
553 \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]
554
555 for \f$k \geq 0\f$. For detailed description see
556 <A HREF="http://mathworld.wolfram.com/Studentst-Distribution.html">
557 Mathworld</A>.
558
559 @ingroup PdfFunc
560
561 */
562
563 inline double tdistribution_pdf(double x, double r, double x0 = 0) {
564 // Inlined to enable clad-auto-derivation for this function.
565
566 return (std::exp (ROOT::Math::lgamma((r + 1.0)/2.0) - ROOT::Math::lgamma(r/2.0)) / std::sqrt (M_PI * r))
567 * std::pow ((1.0 + (x-x0)*(x-x0)/r), -(r + 1.0)/2.0);
568
569 }
570
571
572
573
574 /**
575
576 Probability density function of the uniform (flat) distribution.
577
578 \f[ p(x) = {1 \over (b-a)} \f]
579
580 if \f$a \leq x<b\f$ and 0 otherwise. For detailed description see
581 <A HREF="http://mathworld.wolfram.com/UniformDistribution.html">
582 Mathworld</A>.
583
584 @ingroup PdfFunc
585
586 */
587
588 inline double uniform_pdf(double x, double a, double b, double x0 = 0) {
589 // Inlined to enable clad-auto-derivation for this function.
590
591 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! when a=b
592
593 if ((x-x0) < b && (x-x0) >= a)
594 return 1.0/(b - a);
595 return 0.0;
596
597 }
598
599
600
601} // namespace Math
602} // namespace ROOT
603
604
605
606#endif // ROOT_Math_PdfFunc
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define M_PI
Definition Rotated.cxx:105
#define N
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
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.
double log1p(double x)
declarations for functions which are not implemented by some compilers
Definition Math.h:98
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
TMarker m
Definition textangle.C:8