*   * @param x argument to the polynomial.   * @param coef the coefficients of the polynomial.   * @param N the degree of the polynomial.   */  static double polevl( double x, double coef[], int N ) {    double ans;    ans = coef[0];      for(int i=1; i<=N; i++) ans = ans*x+coef[i];      return ans;  }  /**   * Returns the Incomplete Gamma function.   * @param a the parameter of the gamma distribution.   * @param x the integration end point.   */  static double incompleteGamma(double a, double x)     {      double ans, ax, c, r;     if( x <= 0 || a <= 0 ) return 0.0;     if( x > 1.0 && x > a ) return 1.0 - incompleteGammaComplement(a,x);    /* Compute  x**a * exp(-x) / gamma(a)  */    ax = a * Math.log(x) - x - lnGamma(a);    if( ax < -MAXLOG ) return( 0.0 );    ax = Math.exp(ax);    /* power series */    r = a;    c = 1.0;    ans = 1.0;    do {      r += 1.0;      c *= x/r;      ans += c;    }    while( c/ans > MACHEP );     return( ans * ax/a );  }  /**   * Returns the Complemented Incomplete Gamma function.   * @param a the parameter of the gamma distribution.   * @param x the integration start point.   */  static double incompleteGammaComplement( double a, double x ) {    double ans, ax, c, yc, r, t, y, z;    double pk, pkm1, pkm2, qk, qkm1, qkm2;    if( x <= 0 || a <= 0 ) return 1.0;     if( x < 1.0 || x < a ) return 1.0 - incompleteGamma(a,x);     ax = a * Math.log(x) - x - lnGamma(a);    if( ax < -MAXLOG ) return 0.0;     ax = Math.exp(ax);     /* continued fraction */    y = 1.0 - a;    z = x + y + 1.0;    c = 0.0;    pkm2 = 1.0;    qkm2 = x;    pkm1 = x + 1.0;    qkm1 = z * x;    ans = pkm1/qkm1;     do {      c += 1.0;      y += 1.0;      z += 2.0;      yc = y * c;      pk = pkm1 * z  -  pkm2 * yc;      qk = qkm1 * z  -  qkm2 * yc;      if( qk != 0 ) {	r = pk/qk;        t = Math.abs( (ans - r)/r );	ans = r;      } else	t = 1.0;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      if( Math.abs(pk) > big ) {	pkm2 *= biginv;        pkm1 *= biginv;	qkm2 *= biginv;	qkm1 *= biginv;      }    } while( t > MACHEP );     return ans * ax;  }  /**   * Returns the Gamma function of the argument.   */  static double gamma(double x) {    double P[] = {      1.60119522476751861407E-4,      1.19135147006586384913E-3,      1.04213797561761569935E-2,      4.76367800457137231464E-2,      2.07448227648435975150E-1,      4.94214826801497100753E-1,      9.99999999999999996796E-1    };    double Q[] = {      -2.31581873324120129819E-5,      5.39605580493303397842E-4,      -4.45641913851797240494E-3,      1.18139785222060435552E-2,      3.58236398605498653373E-2,      -2.34591795718243348568E-1,      7.14304917030273074085E-2,      1.00000000000000000320E0    };    double p, z;    int i;    double q = Math.abs(x);    if( q > 33.0 ) {      if( x < 0.0 ) {	p = Math.floor(q);	if( p == q ) throw new ArithmeticException("gamma: overflow");	i = (int)p;	z = q - p;	if( z > 0.5 ) {	  p += 1.0;	  z = q - p;	}	z = q * Math.sin( Math.PI * z );	if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");	z = Math.abs(z);	z = Math.PI/(z * stirlingFormula(q) );	return -z;      } else {	return stirlingFormula(x);      }    }    z = 1.0;    while( x >= 3.0 ) {      x -= 1.0;      z *= x;    }    while( x < 0.0 ) {      if( x == 0.0 ) {	throw new ArithmeticException("gamma: singular");      } else	if( x > -1.E-9 ) {	  return( z/((1.0 + 0.5772156649015329 * x) * x) );	}      z /= x;      x += 1.0;    }    while( x < 2.0 ) {      if( x == 0.0 ) {	throw new ArithmeticException("gamma: singular");      } else	if( x < 1.e-9 ) {	  return( z/((1.0 + 0.5772156649015329 * x) * x) );	}      z /= x;      x += 1.0;    }    if( (x == 2.0) || (x == 3.0) ) 	return z;    x -= 2.0;    p = polevl( x, P, 6 );    q = polevl( x, Q, 7 );    return  z * p / q;  }  /**   * Returns the Gamma function computed by Stirling's formula.   * The polynomial STIR is valid for 33 <= x <= 172.   */  static double stirlingFormula(double x) {    double STIR[] = {      7.87311395793093628397E-4,      -2.29549961613378126380E-4,      -2.68132617805781232825E-3,      3.47222221605458667310E-3,      8.33333333333482257126E-2,    };    double MAXSTIR = 143.01608;    double w = 1.0/x;    double  y = Math.exp(x);    w = 1.0 + w * polevl( w, STIR, 4 );    if( x > MAXSTIR ) {      /* Avoid overflow in Math.pow() */      double v = Math.pow( x, 0.5 * x - 0.25 );      y = v * (v / y);    } else {      y = Math.pow( x, x - 0.5 ) / y;    }    y = SQTPI * y * w;    return y;  }  /**   * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>.   *   * @param aa the alpha parameter of the beta distribution.   * @param bb the beta parameter of the beta distribution.   * @param xx the integration end point.   */  public static double incompleteBeta( double aa, double bb, double xx ) {    double a, b, t, x, xc, w, y;    boolean flag;    if( aa <= 0.0 || bb <= 0.0 ) throw new       ArithmeticException("ibeta: Domain error!");    if( (xx <= 0.0) || ( xx >= 1.0) ) {      if( xx == 0.0 ) return 0.0;      if( xx == 1.0 ) return 1.0;      throw new ArithmeticException("ibeta: Domain error!");    }    flag = false;    if( (bb * xx) <= 1.0 && xx <= 0.95) {      t = powerSeries(aa, bb, xx);      return t;    }    w = 1.0 - xx;    /* Reverse a and b if x is greater than the mean. */    if( xx > (aa/(aa+bb)) ) {      flag = true;      a = bb;      b = aa;      xc = xx;      x = w;    } else {      a = aa;      b = bb;      xc = w;      x = xx;    }    if( flag  && (b * x) <= 1.0 && x <= 0.95) {      t = powerSeries(a, b, x);      if( t <= MACHEP ) 	t = 1.0 - MACHEP;      else  		        t = 1.0 - t;      return t;    }    /* Choose expansion for better convergence. */    y = x * (a+b-2.0) - (a-1.0);    if( y < 0.0 )      w = incompleteBetaFraction1( a, b, x );    else      w = incompleteBetaFraction2( a, b, x ) / xc;    /* Multiply w by the factor       a      b   _             _     _       x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */    y = a * Math.log(x);    t = b * Math.log(xc);    if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {      t = Math.pow(xc,b);      t *= Math.pow(x,a);      t /= a;      t *= w;      t *= gamma(a+b) / (gamma(a) * gamma(b));      if( flag ) {	if( t <= MACHEP ) 	t = 1.0 - MACHEP;	else  		        t = 1.0 - t;      }      return t;    }    /* Resort to logarithms.  */    y += t + lnGamma(a+b) - lnGamma(a) - lnGamma(b);    y += Math.log(w/a);    if( y < MINLOG )      t = 0.0;    else      t = Math.exp(y);    if( flag ) {      if( t <= MACHEP ) 	t = 1.0 - MACHEP;      else  		        t = 1.0 - t;    }    return t;  }     /**   * Continued fraction expansion #1 for incomplete beta integral.   */  static double incompleteBetaFraction1( double a, double b, double x ) {    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;    double k1, k2, k3, k4, k5, k6, k7, k8;    double r, t, ans, thresh;    int n;    k1 = a;    k2 = a + b;    k3 = a;    k4 = a + 1.0;    k5 = 1.0;    k6 = b - 1.0;    k7 = k4;    k8 = a + 2.0;    pkm2 = 0.0;    qkm2 = 1.0;    pkm1 = 1.0;    qkm1 = 1.0;    ans = 1.0;    r = 1.0;    n = 0;    thresh = 3.0 * MACHEP;    do {      xk = -( x * k1 * k2 )/( k3 * k4 );      pk = pkm1 +  pkm2 * xk;      qk = qkm1 +  qkm2 * xk;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      xk = ( x * k5 * k6 )/( k7 * k8 );      pk = pkm1 +  pkm2 * xk;      qk = qkm1 +  qkm2 * xk;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      if( qk != 0 )		r = pk/qk;      if( r != 0 ) {	t = Math.abs( (ans - r)/r );	ans = r;      }	else	t = 1.0;      if( t < thresh ) return ans;      k1 += 1.0;      k2 += 1.0;      k3 += 2.0;      k4 += 2.0;      k5 += 1.0;      k6 -= 1.0;      k7 += 2.0;      k8 += 2.0;      if( (Math.abs(qk) + Math.abs(pk)) > big ) {	pkm2 *= biginv;	pkm1 *= biginv;	qkm2 *= biginv;	qkm1 *= biginv;      }      if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {	pkm2 *= big;	pkm1 *= big;	qkm2 *= big;	qkm1 *= big;      }    } while( ++n < 300 );    return ans;  }     /**   * Continued fraction expansion #2 for incomplete beta integral.   */  static double incompleteBetaFraction2( double a, double b, double x ) {    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;    double k1, k2, k3, k4, k5, k6, k7, k8;    double r, t, ans, z, thresh;    int n;    k1 = a;    k2 = b - 1.0;    k3 = a;    k4 = a + 1.0;    k5 = 1.0;    k6 = a + b;    k7 = a + 1.0;;    k8 = a + 2.0;    pkm2 = 0.0;    qkm2 = 1.0;    pkm1 = 1.0;    qkm1 = 1.0;    z = x / (1.0-x);    ans = 1.0;    r = 1.0;    n = 0;    thresh = 3.0 * MACHEP;    do {      xk = -( z * k1 * k2 )/( k3 * k4 );      pk = pkm1 +  pkm2 * xk;      qk = qkm1 +  qkm2 * xk;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      xk = ( z * k5 * k6 )/( k7 * k8 );      pk = pkm1 +  pkm2 * xk;      qk = qkm1 +  qkm2 * xk;      pkm2 = pkm1;      pkm1 = pk;      qkm2 = qkm1;      qkm1 = qk;      if( qk != 0 )  r = pk/qk;      if( r != 0 ) {	t = Math.abs( (ans - r)/r );	ans = r;      } else	t = 1.0;      if( t < thresh ) return ans;      k1 += 1.0;      k2 -= 1.0;      k3 += 2.0;      k4 += 2.0;      k5 += 1.0;      k6 += 1.0;      k7 += 2.0;      k8 += 2.0;      if( (Math.abs(qk) + Math.abs(pk)) > big ) {	pkm2 *= biginv;	pkm1 *= biginv;	qkm2 *= biginv;	qkm1 *= biginv;      }      if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {	pkm2 *= big;	pkm1 *= big;	qkm2 *= big;	qkm1 *= big;      }    } while( ++n < 300 );    return ans;  }  /**   * Power series for incomplete beta integral.   * Use when b*x is small and x not too close to 1.     */  static double powerSeries( double a, double b, double x ) {    double s, t, u, v, n, t1, z, ai;        ai = 1.0 / a;    u = (1.0 - b) * x;    v = u / (a + 1.0);    t1 = v;    t = u;    n = 2.0;    s = 0.0;    z = MACHEP * ai;    while( Math.abs(v) > z ) {      u = (n - b) * x / n;      t *= u;      v = t / (a + n);      s += v;       n += 1.0;    }    s += t1;    s += ai;    u = a * Math.log(x);    if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {      t = gamma(a+b)/(gamma(a)*gamma(b));      s = s * t * Math.pow(x,a);    } else {      t = lnGamma(a+b) - lnGamma(a) - lnGamma(b) + u + Math.log(s);      if( t < MINLOG ) 	s = 0.0;      else  	            s = Math.exp(t);    }    return s;  }  /**   * Main method for testing this class.   */  public static void main(String[] ops) {    System.out.println("Binomial standard error (0.5, 100): " + 		       Statistics.binomialStandardError(0.5, 100));    System.out.println("Chi-squared probability (2.558, 10): " +		       Statistics.chiSquaredProbability(2.558, 10));    System.out.println("Normal probability (0.2): " +		       Statistics.normalProbability(0.2));    System.out.println("F probability (5.1922, 4, 5): " +		       Statistics.FProbability(5.1922, 4, 5));    System.out.println("lnGamma(6): "+ Statistics.lnGamma(6));  }}