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