// 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.
//
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)
//
}
}
///////////////////////////////////////////////////////////////////////////
-Double_t AliMath::LnGamma(Double_t z)
+Double_t AliMath::LnGamma(Double_t z) const
{
// Computation of ln[gamma(z)] for all z>0.
//
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.
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.
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).
//
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).
//
//--- 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
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.
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.
//
//--- 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);
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.
//
//--- 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)
{
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.
//
//--- 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);
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.
//
//--- 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)
{
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.
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.
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 : <chi2>=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 : <T>=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 : <F>=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 : <k>=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 : <n>=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 : <n>=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 : <Chi2>=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. <Chi2> 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-<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 (chi2<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (chi2<mean) sides=-2;
+ }
+
+ Double_t val=0;
+ if (sigma) // P-value in units of sigma
+ {
+ Double_t s=sqrt(double(2*ndf));
+ val=(chi2-mean)/s;
+ }
+ else // P-value from tail contents
+ {
+ if (sides>0) // 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 : <T>=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. <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-<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 (t<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (t<mean) sides=-2;
+ }
+
+ Double_t val=0;
+ if (sigma) // P-value in units of sigma
+ {
+ if (ndf>2) // 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 : <F>=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. <F> 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-<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 (f<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (f<mean) sides=-2;
+ }
+
+ Double_t val=0;
+ if (sigma) // P-value in units of sigma
+ {
+ // Sigma is only defined for ndf2>4
+ 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 : <K>=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. <K> 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-<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 (k<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (k<mean) sides=-2;
+ }
+
+ Double_t val=0;
+
+ if (sigma) // P-value in units of sigma
+ {
+ Double_t s=sqrt(double(n)*p*(1.-p));
+ val=(double(k)-mean)/s;
+ }
+ else // P-value from tail contents
+ {
+ if (sides>0)
+ {
+ 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 : <K>=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. <K> 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-<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 (k<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (k<mean) sides=-2;
+ }
+
+ Double_t val=0;
+
+ if (sigma) // P-value in units of sigma
+ {
+ Double_t s=sqrt(mu);
+ val=(double(k)-mean)/s;
+ }
+ else // P-value from tail contents
+ {
+ if (sides>0) // 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 : <N>=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. <N> 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-<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 (n<mean) sides=-1;
+ }
+
+ if (sides==3) // Automatic two-sided test
+ {
+ sides=2;
+ if (n<mean) sides=-2;
+ }
+
+ Double_t val=0;
+
+ if (sigma) // P-value in units of sigma
+ {
+ Double_t s=sqrt(double(k)*(1.-p)/(p*p));
+ val=(double(n)-mean)/s;
+ }
+ else // P-value from tail contents
+ {
+ if (sides>0) // 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<n; i++)
+ {
+ val+=TMath::Binomial(i-1,k-1)*pow(p,k)*pow(1.-p,i-k);
+ }
+ val=1.-val;
+ }
+ }
+ else // Lower tail
+ {
+ if (!mode) // Use Incomplete Beta function
+ {
+ val=TMath::BetaIncomplete(p,k,n-k);
+ }
+ else
+ {
+ for (Int_t j=1; j<n; j++)
+ {
+ val+=TMath::Binomial(j-1,k-1)*pow(p,k)*pow(1.-p,j-k);
+ }
+ }
+ }
+ }
+
+ if (abs(sides)==2) val=val*2.;
+
+ return val;
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t AliMath::Nfac(Int_t n,Int_t mode) const
+{
+// Compute n!.
+// The algorithm can be selected by the "mode" input argument.
+//
+// mode : 0 ==> 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;
+}
+///////////////////////////////////////////////////////////////////////////
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff