X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=RALICE%2FAliMath.cxx;h=2c0864fa4ba51fbb23feaa2a005e8279eedfb09d;hb=cf3100fab50e161050cd0b4ea2023986e164a199;hp=fcf42eaecef6a394f4f5daea75ab0bf6a63d67c6;hpb=c72198f1e40c1ac6d987530ef7440da78e24eab8;p=u%2Fmrichter%2FAliRoot.git diff --git a/RALICE/AliMath.cxx b/RALICE/AliMath.cxx index fcf42eaecef..2c0864fa4ba 100644 --- a/RALICE/AliMath.cxx +++ b/RALICE/AliMath.cxx @@ -49,12 +49,12 @@ AliMath::~AliMath() // Destructor } /////////////////////////////////////////////////////////////////////////// -AliMath::AliMath(AliMath& m) : TObject(m) +AliMath::AliMath(const AliMath& m) : TObject(m) { // Copy constructor } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::Gamma(Double_t z) +Double_t AliMath::Gamma(Double_t z) const { // Computation of gamma(z) for all z>0. // @@ -75,7 +75,7 @@ Double_t AliMath::Gamma(Double_t z) return exp(v); } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::Gamma(Double_t a,Double_t x) +Double_t AliMath::Gamma(Double_t a,Double_t x) const { // Computation of the incomplete gamma function P(a,x) // @@ -106,7 +106,7 @@ Double_t AliMath::Gamma(Double_t a,Double_t x) } } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::LnGamma(Double_t z) +Double_t AliMath::LnGamma(Double_t z) const { // Computation of ln[gamma(z)] for all z>0. // @@ -149,7 +149,7 @@ Double_t AliMath::LnGamma(Double_t z) return v; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::GamSer(Double_t a,Double_t x) +Double_t AliMath::GamSer(Double_t a,Double_t x) const { // Computation of the incomplete gamma function P(a,x) // via its series representation. @@ -190,7 +190,7 @@ Double_t AliMath::GamSer(Double_t a,Double_t x) return v; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::GamCf(Double_t a,Double_t x) +Double_t AliMath::GamCf(Double_t a,Double_t x) const { // Computation of the incomplete gamma function P(a,x) // via its continued fraction representation. @@ -240,7 +240,7 @@ Double_t AliMath::GamCf(Double_t a,Double_t x) return (1.-v); } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::Erf(Double_t x) +Double_t AliMath::Erf(Double_t x) const { // Computation of the error function erf(x). // @@ -249,7 +249,7 @@ Double_t AliMath::Erf(Double_t x) return (1.-Erfc(x)); } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::Erfc(Double_t x) +Double_t AliMath::Erfc(Double_t x) const { // Computation of the complementary error function erfc(x). // @@ -261,11 +261,11 @@ Double_t AliMath::Erfc(Double_t x) //--- Nve 14-nov-1998 UU-SAP Utrecht // The parameters of the Chebyshev fit - const Double_t a1=-1.26551223, a2=1.00002368, - a3= 0.37409196, a4=0.09678418, - a5=-0.18628806, a6=0.27886807, - a7=-1.13520398, a8=1.48851587, - a9=-0.82215223, a10=0.17087277; + const Double_t ka1=-1.26551223, ka2=1.00002368, + ka3= 0.37409196, ka4=0.09678418, + ka5=-0.18628806, ka6=0.27886807, + ka7=-1.13520398, ka8=1.48851587, + ka9=-0.82215223, ka10=0.17087277; Double_t v=1.; // The return value @@ -276,20 +276,33 @@ Double_t AliMath::Erfc(Double_t x) Double_t t=1./(1.+0.5*z); v=t*exp((-z*z) - +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10))))))))); + +ka1+t*(ka2+t*(ka3+t*(ka4+t*(ka5+t*(ka6+t*(ka7+t*(ka8+t*(ka9+t*ka10))))))))); if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x) return v; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::Prob(Double_t chi2,Int_t ndf) +Double_t AliMath::Prob(Double_t chi2,Int_t ndf,Int_t mode) const { // Computation of the probability for a certain Chi-squared (chi2) // and number of degrees of freedom (ndf). // -// Calculations are based on the incomplete gamma function P(a,x), -// where a=ndf/2 and x=chi2/2. +// A more clear and flexible facility is offered by Chi2Pvalue. +// +// According to the value of the parameter "mode" various algorithms +// can be selected. +// +// mode = 0 : Calculations are based on the incomplete gamma function P(a,x), +// where a=ndf/2 and x=chi2/2. +// +// 1 : Same as for mode=0. However, in case ndf=1 an exact expression +// based on the error function Erf() is used. +// +// 2 : Same as for mode=0. However, in case ndf>30 a Gaussian approximation +// is used instead of the gamma function. +// +// When invoked as Prob(chi2,ndf) the default mode=1 is used. // // P(a,x) represents the probability that the observed Chi-squared // for a correct model should be less than the value chi2. @@ -313,31 +326,36 @@ Double_t AliMath::Prob(Double_t chi2,Int_t ndf) return 1; } } + + Double_t v=-1.; + + switch (mode) + { + case 1: // Exact expression for ndf=1 as alternative for the gamma function + if (ndf==1) v=1.-Erf(sqrt(chi2)/sqrt(2.)); + break; -// Alternative which is exact -// This code may be activated in case the gamma function gives problems -// if (ndf==1) -// { -// Double_t v=1.-Erf(sqrt(chi2)/sqrt(2.)); -// return v; -// } - -// Gaussian approximation for large ndf -// This code may be activated in case the gamma function shows a problem -// Double_t q=sqrt(2.*chi2)-sqrt(double(2*ndf-1)); -// if (n>30 && q>0.) -// { -// Double_t v=0.5*(1.-Erf(q/sqrt(2.))); -// return v; -// } + case 2: // Gaussian approximation for large ndf (i.e. ndf>30) as alternative for the gamma function + if (ndf>30) + { + Double_t q=sqrt(2.*chi2)-sqrt(double(2*ndf-1)); + if (q>0.) v=0.5*(1.-Erf(q/sqrt(2.))); + } + break; + } - // Evaluate the incomplete gamma function - Double_t a=double(ndf)/2.; - Double_t x=chi2/2.; - return (1.-Gamma(a,x)); + if (v<0.) + { + // Evaluate the incomplete gamma function + Double_t a=double(ndf)/2.; + Double_t x=chi2/2.; + v=1.-Gamma(a,x); + } + + return v; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselI0(Double_t x) +Double_t AliMath::BesselI0(Double_t x) const { // Computation of the modified Bessel function I_0(x) for any real x. // @@ -350,12 +368,12 @@ Double_t AliMath::BesselI0(Double_t x) //--- NvE 12-mar-2000 UU-SAP Utrecht // Parameters of the polynomial approximation - const Double_t p1=1.0, p2=3.5156229, p3=3.0899424, - p4=1.2067492, p5=0.2659732, p6=3.60768e-2, p7=4.5813e-3; + const Double_t kp1=1.0, kp2=3.5156229, kp3=3.0899424, + kp4=1.2067492, kp5=0.2659732, kp6=3.60768e-2, kp7=4.5813e-3; - const Double_t q1= 0.39894228, q2= 1.328592e-2, q3= 2.25319e-3, - q4=-1.57565e-3, q5= 9.16281e-3, q6=-2.057706e-2, - q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3; + const Double_t kq1= 0.39894228, kq2= 1.328592e-2, kq3= 2.25319e-3, + kq4=-1.57565e-3, kq5= 9.16281e-3, kq6=-2.057706e-2, + kq7= 2.635537e-2, kq8=-1.647633e-2, kq9= 3.92377e-3; Double_t ax=fabs(x); @@ -364,18 +382,19 @@ Double_t AliMath::BesselI0(Double_t x) if (ax < 3.75) { y=pow(x/3.75,2); - result=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))); + result=kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7))))); } else { y=3.75/ax; - result=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))); + result=(exp(ax)/sqrt(ax)) + *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*(kq7+y*(kq8+y*kq9)))))))); } return result; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselK0(Double_t x) +Double_t AliMath::BesselK0(Double_t x) const { // Computation of the modified Bessel function K_0(x) for positive real x. // @@ -388,11 +407,11 @@ Double_t AliMath::BesselK0(Double_t x) //--- NvE 12-mar-2000 UU-SAP Utrecht // Parameters of the polynomial approximation - const Double_t p1=-0.57721566, p2=0.42278420, p3=0.23069756, - p4= 3.488590e-2, p5=2.62698e-3, p6=1.0750e-4, p7=7.4e-5; + const Double_t kp1=-0.57721566, kp2=0.42278420, kp3=0.23069756, + kp4= 3.488590e-2, kp5=2.62698e-3, kp6=1.0750e-4, kp7=7.4e-6; - const Double_t q1= 1.25331414, q2=-7.832358e-2, q3= 2.189568e-2, - q4=-1.062446e-2, q5= 5.87872e-3, q6=-2.51540e-3, q7=5.3208e-4; + const Double_t kq1= 1.25331414, kq2=-7.832358e-2, kq3= 2.189568e-2, + kq4=-1.062446e-2, kq5= 5.87872e-3, kq6=-2.51540e-3, kq7=5.3208e-4; if (x <= 0) { @@ -405,18 +424,20 @@ Double_t AliMath::BesselK0(Double_t x) if (x <= 2) { y=x*x/4.; - result=(-log(x/2.)*BesselI0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))); + result=(-log(x/2.)*BesselI0(x)) + +(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7)))))); } else { y=2./x; - result=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7)))))); + result=(exp(-x)/sqrt(x)) + *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*kq7)))))); } return result; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselI1(Double_t x) +Double_t AliMath::BesselI1(Double_t x) const { // Computation of the modified Bessel function I_1(x) for any real x. // @@ -429,12 +450,12 @@ Double_t AliMath::BesselI1(Double_t x) //--- NvE 12-mar-2000 UU-SAP Utrecht // Parameters of the polynomial approximation - const Double_t p1=0.5, p2=0.87890594, p3=0.51498869, - p4=0.15084934, p5=2.658733e-2, p6=3.01532e-3, p7=3.2411e-4; + const Double_t kp1=0.5, kp2=0.87890594, kp3=0.51498869, + kp4=0.15084934, kp5=2.658733e-2, kp6=3.01532e-3, kp7=3.2411e-4; - const Double_t q1= 0.39894228, q2=-3.988024e-2, q3=-3.62018e-3, - q4= 1.63801e-3, q5=-1.031555e-2, q6= 2.282967e-2, - q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3; + const Double_t kq1= 0.39894228, kq2=-3.988024e-2, kq3=-3.62018e-3, + kq4= 1.63801e-3, kq5=-1.031555e-2, kq6= 2.282967e-2, + kq7=-2.895312e-2, kq8= 1.787654e-2, kq9=-4.20059e-3; Double_t ax=fabs(x); @@ -443,19 +464,20 @@ Double_t AliMath::BesselI1(Double_t x) if (ax < 3.75) { y=pow(x/3.75,2); - result=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))); + result=x*(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7)))))); } else { y=3.75/ax; - result=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))); + result=(exp(ax)/sqrt(ax)) + *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*(kq7+y*(kq8+y*kq9)))))))); if (x < 0) result=-result; } return result; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselK1(Double_t x) +Double_t AliMath::BesselK1(Double_t x) const { // Computation of the modified Bessel function K_1(x) for positive real x. // @@ -468,11 +490,11 @@ Double_t AliMath::BesselK1(Double_t x) //--- NvE 12-mar-2000 UU-SAP Utrecht // Parameters of the polynomial approximation - const Double_t p1= 1., p2= 0.15443144, p3=-0.67278579, - p4=-0.18156897, p5=-1.919402e-2, p6=-1.10404e-3, p7=-4.686e-5; + const Double_t kp1= 1., kp2= 0.15443144, kp3=-0.67278579, + kp4=-0.18156897, kp5=-1.919402e-2, kp6=-1.10404e-3, kp7=-4.686e-5; - const Double_t q1= 1.25331414, q2= 0.23498619, q3=-3.655620e-2, - q4= 1.504268e-2, q5=-7.80353e-3, q6= 3.25614e-3, q7=-6.8245e-4; + const Double_t kq1= 1.25331414, kq2= 0.23498619, kq3=-3.655620e-2, + kq4= 1.504268e-2, kq5=-7.80353e-3, kq6= 3.25614e-3, kq7=-6.8245e-4; if (x <= 0) { @@ -485,18 +507,20 @@ Double_t AliMath::BesselK1(Double_t x) if (x <= 2) { y=x*x/4.; - result=(log(x/2.)*BesselI1(x))+(1./x)*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))); + result=(log(x/2.)*BesselI1(x)) + +(1./x)*(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7)))))); } else { y=2./x; - result=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*q7)))))); + result=(exp(-x)/sqrt(x)) + *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*kq7)))))); } return result; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselK(Int_t n,Double_t x) +Double_t AliMath::BesselK(Int_t n,Double_t x) const { // Computation of the Integer Order Modified Bessel function K_n(x) // for n=0,1,2,... and positive real x. @@ -534,7 +558,7 @@ Double_t AliMath::BesselK(Int_t n,Double_t x) return bk; } /////////////////////////////////////////////////////////////////////////// -Double_t AliMath::BesselI(Int_t n,Double_t x) +Double_t AliMath::BesselI(Int_t n,Double_t x) const { // Computation of the Integer Order Modified Bessel function I_n(x) // for n=0,1,2,... and any real x. @@ -587,3 +611,817 @@ Double_t AliMath::BesselI(Int_t n,Double_t x) return result; } /////////////////////////////////////////////////////////////////////////// +TF1 AliMath::Chi2Dist(Int_t ndf) const +{ +// Provide the Chi-squared distribution function corresponding to the +// specified ndf degrees of freedom. +// +// Details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =ndf Var(chi2)=2*ndf + + TF1 chi2dist("chi2dist","1./(TMath::Gamma([0]/2.)*pow(2,[0]/2.))*pow(x,[0]/2.-1.)*exp(-x/2.)"); + TString title="#chi^{2} distribution for ndf = "; + title+=ndf; + chi2dist.SetTitle(title.Data()); + chi2dist.SetParName(0,"ndf"); + chi2dist.SetParameter(0,ndf); + + return chi2dist; +} +/////////////////////////////////////////////////////////////////////////// +TF1 AliMath::StudentDist(Double_t ndf) const +{ +// Provide the Student's T distribution function corresponding to the +// specified ndf degrees of freedom. +// +// In a frequentist approach, the Student's T distribution is particularly +// useful in making inferences about the mean of an underlying population +// based on the data from a random sample. +// +// In a Bayesian context it is used to characterise the posterior PDF +// for a particular state of information. +// +// Note : ndf is not restricted to integer values +// +// Details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =0 Var(T)=ndf/(ndf-2) + + TF1 tdist("tdist","(TMath::Gamma(([0]+1.)/2.)/(sqrt(pi*[0])*TMath::Gamma([0]/2.)))*pow(1.+x*x/[0],-([0]+1.)/2.)"); + TString title="Student's t distribution for ndf = "; + title+=ndf; + tdist.SetTitle(title.Data()); + tdist.SetParName(0,"ndf"); + tdist.SetParameter(0,ndf); + + return tdist; +} +/////////////////////////////////////////////////////////////////////////// +TF1 AliMath::FratioDist(Int_t ndf1,Int_t ndf2) const +{ +// Provide the F (ratio) distribution function corresponding to the +// specified ndf1 and ndf2 degrees of freedom of the two samples. +// +// In a frequentist approach, the F (ratio) distribution is particularly useful +// in making inferences about the ratio of the variances of two underlying +// populations based on a measurement of the variance of a random sample taken +// from each one of the two populations. +// So the F test provides a means to investigate the degree of equality of +// two populations. +// +// Details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =ndf2/(ndf2-2) Var(F)=2*ndf2*ndf2*(ndf2+ndf1-2)/(ndf1*(ndf2-1)*(ndf2-1)*(ndf2-4)) + + TF1 fdist("fdist", + "(TMath::Gamma(([0]+[1])/2.)/(TMath::Gamma([0]/2.)*TMath::Gamma([1]/2.)))*pow([0]/[1],[0]/2.)*pow(x,([0]-2.)/2.)/pow(1.+x*[0]/[1],([0]+[1])/2.)"); + TString title="F (ratio) distribution for ndf1 = "; + title+=ndf1; + title+=" ndf2 = "; + title+=ndf2; + fdist.SetTitle(title.Data()); + fdist.SetParName(0,"ndf1"); + fdist.SetParameter(0,ndf1); + fdist.SetParName(1,"ndf2"); + fdist.SetParameter(1,ndf2); + + return fdist; +} +/////////////////////////////////////////////////////////////////////////// +TF1 AliMath::BinomialDist(Int_t n,Double_t p) const +{ +// Provide the Binomial distribution function corresponding to the +// specified number of trials n and probability p of success. +// +// p(k|n,p) = probability to obtain exactly k successes in n trials. +// +// Details can be found in the excelent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =n*p Var(k)=n*p*(1-p) + + TF1 bindist("bindist","TMath::Binomial(int([0]),int(x))*pow([1],int(x))*pow(1.-[1],int([0])-int(x))"); + TString title="Binomial distribution for n = "; + title+=n; + title+=" p = "; + title+=p; + bindist.SetTitle(title.Data()); + bindist.SetParName(0,"n"); + bindist.SetParameter(0,n); + bindist.SetParName(1,"p"); + bindist.SetParameter(1,p); + + return bindist; +} +/////////////////////////////////////////////////////////////////////////// +TF1 AliMath::NegBinomialDist(Int_t k,Double_t p) const +{ +// Provide the Negative Binomial distribution function corresponding to the +// specified number of successes k and probability p of success. +// +// p(n|k,p) = probability to have reached k successes after n trials. +// +// Details can be found in the excelent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =k*(1-p)/p Var(n)=k*(1-p)/(p*p) + + TF1 negbindist("negbindist","TMath::Binomial(int(x)-1,int([0])-1)*pow([1],int([0]))*pow(1.-[1],int(x)-int([0]))"); + TString title="Negative Binomial distribution for k = "; + title+=k; + title+=" p = "; + title+=p; + negbindist.SetTitle(title.Data()); + negbindist.SetParName(0,"k"); + negbindist.SetParameter(0,k); + negbindist.SetParName(1,"p"); + negbindist.SetParameter(1,p); + + return negbindist; +} +/////////////////////////////////////////////////////////////////////////// +TF1 AliMath::PoissonDist(Double_t mu) const +{ +// Provide the Poisson distribution function corresponding to the +// specified average rate (in time or space) mu of occurrences. +// +// p(n|mu) = probability for n occurrences in a certain time or space interval. +// +// Details can be found in the excelent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =mu Var(n)=mu + + TF1 poissdist("poissdist","exp(-[0])*pow([0],int(x))/TMath::Factorial(int(x))"); + TString title="Poisson distribution for mu = "; + title+=mu; + poissdist.SetTitle(title.Data()); + poissdist.SetParName(0,"mu"); + poissdist.SetParameter(0,mu); + + return poissdist; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::Chi2Pvalue(Double_t chi2,Int_t ndf,Int_t sides,Int_t sigma,Int_t mode) const +{ +// Computation of the P-value for a certain specified Chi-squared (chi2) value +// for a Chi-squared distribution with ndf degrees of freedom. +// +// The P-value for a certain Chi-squared value chi2 corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield a Chi-squared value exceeding (less than) the value chi2 for an +// upper (lower) tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =ndf Var(Chi2)=2*ndf +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input chi2 value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference chi2- expressed in units of sigma +// Note : This difference may be negative. +// +// According to the value of the parameter "mode" various algorithms +// can be selected. +// +// mode = 0 : Calculations are based on the incomplete gamma function. +// +// 1 : Same as for mode=0. However, in case ndf=1 an exact expression +// based on the error function Erf() is used. +// +// 2 : Same as for mode=0. However, in case ndf>30 a Gaussian approximation +// is used instead of the gamma function. +// +// The default values are sides=0, sigma=0 and mode=1. +// +//--- NvE 21-may-2005 Utrecht University + + if (ndf<=0) return 0; + + Double_t mean=ndf; + + if (!sides) // Automatic one-sided test + { + sides=1; + if (chi20) // Upper tail + { + val=Prob(chi2,ndf,mode); + } + else // Lower tail + { + val=1.-Prob(chi2,ndf,mode); + } + } + + if (abs(sides)==2) val=val*2.; + + return val; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::StudentPvalue(Double_t t,Double_t ndf,Int_t sides,Int_t sigma) const +{ +// Computation of the P-value for a certain specified Student's t value +// for a Student's T distribution with ndf degrees of freedom. +// +// In a frequentist approach, the Student's T distribution is particularly +// useful in making inferences about the mean of an underlying population +// based on the data from a random sample. +// +// The P-value for a certain t value corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield a T value exceeding (less than) the value t for an upper (lower) +// tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =0 Var(T)=ndf/(ndf-2) +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input t value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference t- expressed in units of sigma +// Note : This difference may be negative and sigma +// is only defined for ndf>2. +// +// The default values are sides=0 and sigma=0. +// +//--- NvE 21-may-2005 Utrecht University + + if (ndf<=0) return 0; + + Double_t mean=0; + + if (!sides) // Automatic one-sided test + { + sides=1; + if (t2) // Sigma is only defined for ndf>2 + { + Double_t s=sqrt(ndf/(ndf-2.)); + val=t/s; + } + } + else // P-value from tail contents + { + if (sides>0) // Upper tail + { + val=1.-TMath::StudentI(t,ndf); + } + else // Lower tail + { + val=TMath::StudentI(t,ndf); + } + } + + if (abs(sides)==2) val=val*2.; + + return val; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::FratioPvalue(Double_t f,Int_t ndf1,Int_t ndf2,Int_t sides,Int_t sigma) const +{ +// Computation of the P-value for a certain specified F ratio f value +// for an F (ratio) distribution with ndf1 and ndf2 degrees of freedom +// for the two samples X,Y respectively to be compared in the ratio X/Y. +// +// In a frequentist approach, the F (ratio) distribution is particularly useful +// in making inferences about the ratio of the variances of two underlying +// populations based on a measurement of the variance of a random sample taken +// from each one of the two populations. +// So the F test provides a means to investigate the degree of equality of +// two populations. +// +// The P-value for a certain f value corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield an F value exceeding (less than) the value f for an upper (lower) +// tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =ndf2/(ndf2-2) Var(F)=2*ndf2*ndf2*(ndf2+ndf1-2)/(ndf1*(ndf2-1)*(ndf2-1)*(ndf2-4)) +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input f value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference f- expressed in units of sigma +// Note : This difference may be negative and sigma +// is only defined for ndf2>4. +// +// The default values are sides=0 and sigma=0. +// +//--- NvE 21-may-2005 Utrecht University + + if (ndf1<=0 || ndf2<=0 || f<=0) return 0; + + Double_t mean=double(ndf2/(ndf2-2)); + + if (!sides) // Automatic one-sided test + { + sides=1; + if (f4 + if (ndf2>4) + { + Double_t s=sqrt(double(ndf2*ndf2*(2*ndf2+2*ndf1-4))/double(ndf1*pow(double(ndf2-1),2)*(ndf2-4))); + val=(f-mean)/s; + } + } + else // P-value from tail contents + { + if (sides>0) // Upper tail + { + val=1.-TMath::FDistI(f,ndf1,ndf2); + } + else // Lower tail + { + val=TMath::FDistI(f,ndf1,ndf2); + } + } + + if (abs(sides)==2) val=val*2.; + + return val; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::BinomialPvalue(Int_t k,Int_t n,Double_t p,Int_t sides,Int_t sigma,Int_t mode) const +{ +// Computation of the P-value for a certain specified number of successes k +// for a Binomial distribution with n trials and success probability p. +// +// The P-value for a certain number of successes k corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield a number of successes exceeding (less than) the value k for an +// upper (lower) tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =n*p Var(K)=n*p*(1-p) +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input k value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference k- expressed in units of sigma +// Note : This difference may be negative. +// +// mode = 0 : Incomplete Beta function will be used to calculate the tail content. +// 1 : Straightforward summation of the Binomial terms is used. +// +// The Incomplete Beta function in general provides the most accurate values. +// +// The default values are sides=0, sigma=0 and mode=0. +// +//--- NvE 24-may-2005 Utrecht University + + Double_t mean=double(n)*p; + + if (!sides) // Automatic one-sided test + { + sides=1; + if (k0) + { + if (!mode) // Use Incomplete Beta function + { + val=TMath::BetaIncomplete(p,k+1,n-k); + } + else // Use straightforward summation + { + for (Int_t i=k+1; i<=n; i++) + { + val+=TMath::Binomial(n,i)*pow(p,i)*pow(1.-p,n-i); + } + } + } + else + { + if (!mode) // Use Incomplete Beta function + { + val=1.-TMath::BetaIncomplete(p,k+1,n-k); + } + else + { + for (Int_t j=0; j<=k; j++) + { + val+=TMath::Binomial(n,j)*pow(p,j)*pow(1.-p,n-j); + } + } + } + } + + if (abs(sides)==2) val=val*2.; + + return val; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::PoissonPvalue(Int_t k,Double_t mu,Int_t sides,Int_t sigma) const +{ +// Computation of the P-value for a certain specified number of occurrences k +// for a Poisson distribution with a specified average rate (in time or space) +// mu of occurrences. +// +// The P-value for a certain number of occurrences k corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield a number of occurrences exceeding (less than) the value k for an +// upper (lower) tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =mu Var(K)=mu +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input k value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference k- expressed in units of sigma +// Note : This difference may be negative. +// +// The default values are sides=0 and sigma=0. +// +// Note : The tail contents are given by the incomplete Gamma function P(a,x). +// Lower tail contents = 1-P(k,mu) +// Upper tail contents = P(k,mu) +// +//--- NvE 24-may-2005 Utrecht University + + Double_t mean=mu; + + if (!sides) // Automatic one-sided test + { + sides=1; + if (k0) // Upper tail + { + val=Gamma(k,mu); + } + else // Lower tail + { + val=1.-Gamma(k,mu); + } + } + + if (abs(sides)==2) val=val*2.; + + return val; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::NegBinomialPvalue(Int_t n,Int_t k,Double_t p,Int_t sides,Int_t sigma,Int_t mode) const +{ +// Computation of the P-value for a certain specified number of trials n +// for a Negative Binomial distribution where exactly k successes are to +// be reached which have each a probability p. +// +// The P-value for a certain number of trials n corresponds to the fraction of +// repeatedly drawn equivalent samples from a certain population, which is expected +// to yield a number of trials exceeding (less than) the value n for an +// upper (lower) tail test in case a certain hypothesis is true. +// +// Further details can be found in the excellent textbook of Phil Gregory +// "Bayesian Logical Data Analysis for the Physical Sciences". +// +// Note : =k*(1-p)/p Var(N)=k*(1-p)/(p*p) +// +// With the "sides" parameter a one-sided or two-sided test can be selected +// using either the upper or lower tail contents. +// In case of automatic upper/lower selection the decision is made on basis +// of the location of the input n value w.r.t. of the distribution. +// +// sides = 1 : One-sided test using the upper tail contents +// 2 : Two-sided test using the upper tail contents +// -1 : One-sided test using the lower tail contents +// -2 : Two-sided test using the lower tail contents +// 0 : One-sided test using the auto-selected upper or lower tail contents +// 3 : Two-sided test using the auto-selected upper or lower tail contents +// +// The argument "sigma" allows for the following return values : +// +// sigma = 0 : P-value is returned as the above specified fraction +// 1 : The difference n- expressed in units of sigma +// Note : This difference may be negative. +// +// mode = 0 : Incomplete Beta function will be used to calculate the tail content. +// 1 : Straightforward summation of the Negative Binomial terms is used. +// +// The Incomplete Beta function in general provides the most accurate values. +// +// The default values are sides=0, sigma=0 and mode=0. +// +//--- NvE 24-may-2005 Utrecht University + + Double_t mean=double(k)*(1.-p)/p; + + if (!sides) // Automatic one-sided test + { + sides=1; + if (n0) // Upper tail + { + if (!mode) // Use Incomplete Beta function + { + val=1.-TMath::BetaIncomplete(p,k,n-k); + } + else // Use straightforward summation + { + for (Int_t i=1; i Calculation by means of straightforward multiplication +// : 1 ==> Calculation by means of Stirling's approximation +// : 2 ==> Calculation by means of n!=Gamma(n+1) +// +// For large n the calculation modes 1 and 2 will in general be faster. +// By default mode=0 is used. +// For n<0 the value 0 will be returned. +// +// Note : Because of Double_t value overflow the maximum value is n=170. +// +//--- NvE 20-jan-2007 Utrecht University + + if (n<0) return 0; + if (n==0) return 1; + + Double_t twopi=2.*acos(-1.); + Double_t z=0; + Double_t nfac=1; + Int_t i=n; + + switch (mode) + { + case 0: // Straightforward multiplication + while (i>1) + { + nfac*=Double_t(i); + i--; + } + break; + + case 1: // Stirling's approximation + z=n; + nfac=sqrt(twopi)*pow(z,z+0.5)*exp(-z)*(1.+1./(12.*z)); + break; + + case 2: // Use of Gamma(n+1) + z=n+1; + nfac=Gamma(z); + break; + + default: + nfac=0; + break; + } + + return nfac; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::LnNfac(Int_t n,Int_t mode) const +{ +// Compute ln(n!). +// The algorithm can be selected by the "mode" input argument. +// +// mode : 0 ==> Calculation via evaluation of n! followed by taking ln(n!) +// : 1 ==> Calculation via Stirling's approximation ln(n!)=0.5*ln(2*pi)+(n+0.5)*ln(n)-n+1/(12*n) +// : 2 ==> Calculation by means of ln(n!)=LnGamma(n+1) +// +// Note : Because of Double_t value overflow the maximum value is n=170 for mode=0. +// +// For mode=2 rather accurate results are obtained for both small and large n. +// By default mode=2 is used. +// For n<1 the value 0 will be returned. +// +//--- NvE 20-jan-2007 Utrecht University + + if (n<=1) return 0; + + Double_t twopi=2.*acos(-1.); + Double_t z=0; + Double_t lognfac=0; + + switch (mode) + { + case 0: // Straightforward ln(n!) + z=Nfac(n); + lognfac=log(z); + break; + + case 1: // Stirling's approximation + z=n; + lognfac=0.5*log(twopi)+(z+0.5)*log(z)-z+1./(12.*z); + break; + + case 2: // Use of LnGamma(n+1) + z=n+1; + lognfac=LnGamma(z); + break; + + default: + lognfac=0; + break; + } + + return lognfac; +} +/////////////////////////////////////////////////////////////////////////// +Double_t AliMath::LogNfac(Int_t n,Int_t mode) const +{ +// Compute log_10(n!). +// First ln(n!) is evaluated via invokation of LnNfac(n,mode). +// Then the algorithm log_10(z)=ln(z)*log_10(e) is used. +// +//--- NvE 20-jan-2007 Utrecht University + + Double_t e=exp(1.); + + Double_t val=LnNfac(n,mode); + val*=log10(e); + + return val; +} +///////////////////////////////////////////////////////////////////////////