// CHI SQUARE DISTRIBUTION AND CHI SQUARE FUNCTIONS
// Chi-Square Cumulative Distribution Function
// probability that an observed chi-square value for a correct model should be less than chiSquare
// nu = the degrees of freedom
public static double chiSquareCDF(double chiSquare, int nu){
if(nu<=0)throw new IllegalArgumentException("The degrees of freedom [nu], " + nu + ", must be greater than zero");
return Stat.incompleteGamma((double)nu/2.0D, chiSquare/2.0D);
}
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt
// Retained for backward compatibility
public static double incompleteGamma(double a, double x){
return regularisedGammaFunction(a, x);
}
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt
public static double regularisedGammaFunction(double a, double x){
if(a<0.0D || x<0.0D)throw new IllegalArgumentException("\nFunction defined only for a >= 0 and x>=0");
boolean oldIgSupress = Stat.igSupress;
Stat.igSupress = true;
double igf = 0.0D;
if(x!=0){
if(x < a+1.0D){
// Series representation
igf = incompleteGammaSer(a, x);
}
else{
// Continued fraction representation
igf = incompleteGammaFract(a, x);
}
if(igf!=igf)igf = 1.0 - Stat.crigfGaussQuad(a, x);
}
if(igf<0.0)igf = 0.0;
Stat.igSupress = oldIgSupress;
return igf;
}
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt
// Series representation of the function - valid for x < a + 1
public static double incompleteGammaSer(double a, double x){
if(a<0.0D || x<0.0D)throw new IllegalArgumentException("\nFunction defined only for a >= 0 and x>=0");
if(x>=a+1) throw new IllegalArgumentException("\nx >= a+1 use Continued Fraction Representation");
double igf = 0.0D;
if(x!=0.0D){
int i = 0;
boolean check = true;
double acopy = a;
double sum = 1.0/a;
double incr = sum;
double loggamma = Stat.logGamma(a);
while(check){
++i;
++a;
incr *= x/a;
sum += incr;
if(Math.abs(incr) < Math.abs(sum)*Stat.igfeps){
igf = sum*Math.exp(-x+acopy*Math.log(x)- loggamma);
check = false;
}
if(i>=Stat.igfiter){
check=false;
igf = Double.NaN;
if(!Stat.igSupress){
System.out.println("\nMaximum number of iterations were exceeded in Stat.incompleteGammaSer().");
System.out.println("NaN returned.\nIncrement = "+String.valueOf(incr)+".");
System.out.println("Sum = "+String.valueOf(sum)+".\nTolerance = "+String.valueOf(igfeps));
}
}
}
}
return igf;
}
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt
// Continued Fraction representation of the function - valid for x >= a + 1
// This method follows the general procedure used in Numerical Recipes for C,
// The Art of Scientific Computing
// by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
// Cambridge University Press, http://www.nr.com/
public static double incompleteGammaFract(double a, double x){
if(a<0.0D || x<0.0D)throw new IllegalArgumentException("\nFunction defined only for a >= 0 and x>=0");
if(x<a+1) throw new IllegalArgumentException("\nx < a+1 Use Series Representation");
double igf = 0.0D;
if(x!=0.0D){
int i = 0;
double ii = 0;
boolean check = true;
double loggamma = Stat.logGamma(a);
double numer = 0.0D;
double incr = 0.0D;
double denom = x - a + 1.0D;
double first = 1.0D/denom;
double term = 1.0D/FPMIN;
double prod = first;
while(check){
++i;
ii = (double)i;
numer = -ii*(ii - a);
denom += 2.0D;
first = numer*first + denom;
if(Math.abs(first) < Stat.FPMIN){
first = Stat.FPMIN;
}
term = denom + numer/term;
if(Math.abs(term) < Stat.FPMIN){
term = Stat.FPMIN;
}
first = 1.0D/first;
incr = first*term;
prod *= incr;
if(Math.abs(incr - 1.0D) < igfeps)check = false;
if(i>=Stat.igfiter){
check=false;
igf = Double.NaN;
if(!Stat.igSupress){
System.out.println("\nMaximum number of iterations were exceeded in Stat.incompleteGammaFract().");
System.out.println("NaN returned.\nIncrement - 1 = "+String.valueOf(incr-1)+".");
System.out.println("Tolerance = "+String.valueOf(igfeps));
}
}
}
igf = 1.0D - Math.exp(-x+a*Math.log(x)-loggamma)*prod;
}
return igf;
}
// log to base e of the Gamma function
// Lanczos approximation (6 terms)
// Retained for backward compatibility
public static double logGamma(double x){
double xcopy = x;
double fg = 0.0D;
double first = x + lgfGamma + 0.5;
double second = lgfCoeff[0];
if(x>=0.0){
first -= (x + 0.5)*Math.log(first);
for(int i=1; i<=lgfN; i++)second += lgfCoeff[i]/++xcopy;
fg = Math.log(Math.sqrt(2.0*Math.PI)*second/x) - first;
}
else{
fg = Math.PI/(Stat.gamma(1.0D-x)*Math.sin(Math.PI*x));
if(fg!=1.0/0.0 && fg!=-1.0/0.0){
if(fg<0){
throw new IllegalArgumentException("\nThe gamma function is negative");
}
else{
fg = Math.log(fg);
}
}
}
return fg;
}
// Guassian quadrature estimation of the complementary regularised incomplete gamma function
private static double crigfGaussQuad(double a, double x){
double sum = 0.0;
// set increment details
double upper = 100.0*a;
double range = upper - x;
double incr = 0;
if(upper>x && range >100){
incr = range/1000;
}
else{
upper = x + 100.0;
range = 100.0;
incr = 0.1;
}
int nIncr = (int)Math.round(range/incr);
incr = range/nIncr;
// Instantiate integration function
CrigFunct f1 = new CrigFunct();
f1.setA(a);
f1.setB(Stat.logGammaFunction(a));
// Instantiate Integration
Integration intgn1 = new Integration(f1);
double xx = x;
double yy = x + incr;
intgn1.setLimits(xx, yy);
// Perform quadrature
sum = intgn1.gaussQuad(64);
boolean test2 = true;
for(int i=1; i<nIncr; i++){
xx = yy;
yy = xx + incr;
intgn1.setLimits(xx, yy);
sum += intgn1.gaussQuad(64);
}
return sum;
}
卡方累积分布函数
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt // Retained for backward compatibility public static double incompleteGamma(double a, double x){ return regularisedGammaFunction(a, x); }
// Regularised Incomplete Gamma Function P(a,x) = integral from zero to x of (exp(-t)t^(a-1))dt // Series representation of the function - valid for x < a + 1 public static double incompleteGammaSer(double a, double x){ if(a<0.0D || x<0.0D)throw new IllegalArgumentException("\nFunction defined only for a >= 0 and x>=0"); if(x>=a+1) throw new IllegalArgumentException("\nx >= a+1 use Continued Fraction Representation");