0001: package weka.core;
0002:
0003: /**
0004: * Class implementing some distributions, tests, etc. The code is mostly adapted from the CERN
0005: * Jet Java libraries:
0006: *
0007: * Copyright 2001 University of Waikato
0008: * Copyright 1999 CERN - European Organization for Nuclear Research.
0009: * Permission to use, copy, modify, distribute and sell this software and its documentation for
0010: * any purpose is hereby granted without fee, provided that the above copyright notice appear
0011: * in all copies and that both that copyright notice and this permission notice appear in
0012: * supporting documentation.
0013: * CERN and the University of Waikato make no representations about the suitability of this
0014: * software for any purpose. It is provided "as is" without expressed or implied warranty.
0015: *
0016: * @author peter.gedeck@pharma.Novartis.com
0017: * @author wolfgang.hoschek@cern.ch
0018: * @author Eibe Frank (eibe@cs.waikato.ac.nz)
0019: * @author Richard Kirkby (rkirkby@cs.waikato.ac.nz)
0020: * @version $Revision: 1.9 $
0021: */
0022: public class Statistics {
0023:
0024: /** Some constants */
0025: protected static final double MACHEP = 1.11022302462515654042E-16;
0026: protected static final double MAXLOG = 7.09782712893383996732E2;
0027: protected static final double MINLOG = -7.451332191019412076235E2;
0028: protected static final double MAXGAM = 171.624376956302725;
0029: protected static final double SQTPI = 2.50662827463100050242E0;
0030: protected static final double SQRTH = 7.07106781186547524401E-1;
0031: protected static final double LOGPI = 1.14472988584940017414;
0032:
0033: protected static final double big = 4.503599627370496e15;
0034: protected static final double biginv = 2.22044604925031308085e-16;
0035:
0036: /*************************************************
0037: * COEFFICIENTS FOR METHOD normalInverse() *
0038: *************************************************/
0039: /* approximation for 0 <= |y - 0.5| <= 3/8 */
0040: protected static final double P0[] = { -5.99633501014107895267E1,
0041: 9.80010754185999661536E1, -5.66762857469070293439E1,
0042: 1.39312609387279679503E1, -1.23916583867381258016E0, };
0043: protected static final double Q0[] = {
0044: /* 1.00000000000000000000E0,*/
0045: 1.95448858338141759834E0, 4.67627912898881538453E0,
0046: 8.63602421390890590575E1, -2.25462687854119370527E2,
0047: 2.00260212380060660359E2, -8.20372256168333339912E1,
0048: 1.59056225126211695515E1, -1.18331621121330003142E0, };
0049:
0050: /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
0051: * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
0052: */
0053: protected static final double P1[] = { 4.05544892305962419923E0,
0054: 3.15251094599893866154E1, 5.71628192246421288162E1,
0055: 4.40805073893200834700E1, 1.46849561928858024014E1,
0056: 2.18663306850790267539E0, -1.40256079171354495875E-1,
0057: -3.50424626827848203418E-2, -8.57456785154685413611E-4, };
0058: protected static final double Q1[] = {
0059: /* 1.00000000000000000000E0,*/
0060: 1.57799883256466749731E1, 4.53907635128879210584E1,
0061: 4.13172038254672030440E1, 1.50425385692907503408E1,
0062: 2.50464946208309415979E0, -1.42182922854787788574E-1,
0063: -3.80806407691578277194E-2, -9.33259480895457427372E-4, };
0064:
0065: /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
0066: * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
0067: */
0068: protected static final double P2[] = { 3.23774891776946035970E0,
0069: 6.91522889068984211695E0, 3.93881025292474443415E0,
0070: 1.33303460815807542389E0, 2.01485389549179081538E-1,
0071: 1.23716634817820021358E-2, 3.01581553508235416007E-4,
0072: 2.65806974686737550832E-6, 6.23974539184983293730E-9, };
0073: protected static final double Q2[] = {
0074: /* 1.00000000000000000000E0,*/
0075: 6.02427039364742014255E0, 3.67983563856160859403E0,
0076: 1.37702099489081330271E0, 2.16236993594496635890E-1,
0077: 1.34204006088543189037E-2, 3.28014464682127739104E-4,
0078: 2.89247864745380683936E-6, 6.79019408009981274425E-9, };
0079:
0080: /**
0081: * Computes standard error for observed values of a binomial
0082: * random variable.
0083: *
0084: * @param p the probability of success
0085: * @param n the size of the sample
0086: * @return the standard error
0087: */
0088: public static double binomialStandardError(double p, int n) {
0089:
0090: if (n == 0) {
0091: return 0;
0092: }
0093: return Math.sqrt((p * (1 - p)) / (double) n);
0094: }
0095:
0096: /**
0097: * Returns chi-squared probability for given value and degrees
0098: * of freedom. (The probability that the chi-squared variate
0099: * will be greater than x for the given degrees of freedom.)
0100: *
0101: * @param x the value
0102: * @param v the number of degrees of freedom
0103: * @return the chi-squared probability
0104: */
0105: public static double chiSquaredProbability(double x, double v) {
0106:
0107: if (x < 0.0 || v < 1.0)
0108: return 0.0;
0109: return incompleteGammaComplement(v / 2.0, x / 2.0);
0110: }
0111:
0112: /**
0113: * Computes probability of F-ratio.
0114: *
0115: * @param F the F-ratio
0116: * @param df1 the first number of degrees of freedom
0117: * @param df2 the second number of degrees of freedom
0118: * @return the probability of the F-ratio.
0119: */
0120: public static double FProbability(double F, int df1, int df2) {
0121:
0122: return incompleteBeta(df2 / 2.0, df1 / 2.0, df2
0123: / (df2 + df1 * F));
0124: }
0125:
0126: /**
0127: * Returns the area under the Normal (Gaussian) probability density
0128: * function, integrated from minus infinity to <tt>x</tt>
0129: * (assumes mean is zero, variance is one).
0130: * <pre>
0131: * x
0132: * -
0133: * 1 | | 2
0134: * normal(x) = --------- | exp( - t /2 ) dt
0135: * sqrt(2pi) | |
0136: * -
0137: * -inf.
0138: *
0139: * = ( 1 + erf(z) ) / 2
0140: * = erfc(z) / 2
0141: * </pre>
0142: * where <tt>z = x/sqrt(2)</tt>.
0143: * Computation is via the functions <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
0144: *
0145: * @param a the z-value
0146: * @return the probability of the z value according to the normal pdf
0147: */
0148: public static double normalProbability(double a) {
0149:
0150: double x, y, z;
0151:
0152: x = a * SQRTH;
0153: z = Math.abs(x);
0154:
0155: if (z < SQRTH)
0156: y = 0.5 + 0.5 * errorFunction(x);
0157: else {
0158: y = 0.5 * errorFunctionComplemented(z);
0159: if (x > 0)
0160: y = 1.0 - y;
0161: }
0162: return y;
0163: }
0164:
0165: /**
0166: * Returns the value, <tt>x</tt>, for which the area under the
0167: * Normal (Gaussian) probability density function (integrated from
0168: * minus infinity to <tt>x</tt>) is equal to the argument <tt>y</tt>
0169: * (assumes mean is zero, variance is one).
0170: * <p>
0171: * For small arguments <tt>0 < y < exp(-2)</tt>, the program computes
0172: * <tt>z = sqrt( -2.0 * log(y) )</tt>; then the approximation is
0173: * <tt>x = z - log(z)/z - (1/z) P(1/z) / Q(1/z)</tt>.
0174: * There are two rational functions P/Q, one for <tt>0 < y < exp(-32)</tt>
0175: * and the other for <tt>y</tt> up to <tt>exp(-2)</tt>.
0176: * For larger arguments,
0177: * <tt>w = y - 0.5</tt>, and <tt>x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2))</tt>.
0178: *
0179: * @param y0 the area under the normal pdf
0180: * @return the z-value
0181: */
0182: public static double normalInverse(double y0) {
0183:
0184: double x, y, z, y2, x0, x1;
0185: int code;
0186:
0187: final double s2pi = Math.sqrt(2.0 * Math.PI);
0188:
0189: if (y0 <= 0.0)
0190: throw new IllegalArgumentException();
0191: if (y0 >= 1.0)
0192: throw new IllegalArgumentException();
0193: code = 1;
0194: y = y0;
0195: if (y > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */
0196: y = 1.0 - y;
0197: code = 0;
0198: }
0199:
0200: if (y > 0.13533528323661269189) {
0201: y = y - 0.5;
0202: y2 = y * y;
0203: x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8));
0204: x = x * s2pi;
0205: return (x);
0206: }
0207:
0208: x = Math.sqrt(-2.0 * Math.log(y));
0209: x0 = x - Math.log(x) / x;
0210:
0211: z = 1.0 / x;
0212: if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */
0213: x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8);
0214: else
0215: x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8);
0216: x = x0 - x1;
0217: if (code != 0)
0218: x = -x;
0219: return (x);
0220: }
0221:
0222: /**
0223: * Returns natural logarithm of gamma function.
0224: *
0225: * @param x the value
0226: * @return natural logarithm of gamma function
0227: */
0228: public static double lnGamma(double x) {
0229:
0230: double p, q, w, z;
0231:
0232: double A[] = { 8.11614167470508450300E-4,
0233: -5.95061904284301438324E-4, 7.93650340457716943945E-4,
0234: -2.77777777730099687205E-3, 8.33333333333331927722E-2 };
0235: double B[] = { -1.37825152569120859100E3,
0236: -3.88016315134637840924E4, -3.31612992738871184744E5,
0237: -1.16237097492762307383E6, -1.72173700820839662146E6,
0238: -8.53555664245765465627E5 };
0239: double C[] = {
0240: /* 1.00000000000000000000E0, */
0241: -3.51815701436523470549E2, -1.70642106651881159223E4,
0242: -2.20528590553854454839E5, -1.13933444367982507207E6,
0243: -2.53252307177582951285E6, -2.01889141433532773231E6 };
0244:
0245: if (x < -34.0) {
0246: q = -x;
0247: w = lnGamma(q);
0248: p = Math.floor(q);
0249: if (p == q)
0250: throw new ArithmeticException("lnGamma: Overflow");
0251: z = q - p;
0252: if (z > 0.5) {
0253: p += 1.0;
0254: z = p - q;
0255: }
0256: z = q * Math.sin(Math.PI * z);
0257: if (z == 0.0)
0258: throw new ArithmeticException("lnGamma: Overflow");
0259: z = LOGPI - Math.log(z) - w;
0260: return z;
0261: }
0262:
0263: if (x < 13.0) {
0264: z = 1.0;
0265: while (x >= 3.0) {
0266: x -= 1.0;
0267: z *= x;
0268: }
0269: while (x < 2.0) {
0270: if (x == 0.0)
0271: throw new ArithmeticException("lnGamma: Overflow");
0272: z /= x;
0273: x += 1.0;
0274: }
0275: if (z < 0.0)
0276: z = -z;
0277: if (x == 2.0)
0278: return Math.log(z);
0279: x -= 2.0;
0280: p = x * polevl(x, B, 5) / p1evl(x, C, 6);
0281: return (Math.log(z) + p);
0282: }
0283:
0284: if (x > 2.556348e305)
0285: throw new ArithmeticException("lnGamma: Overflow");
0286:
0287: q = (x - 0.5) * Math.log(x) - x + 0.91893853320467274178;
0288:
0289: if (x > 1.0e8)
0290: return (q);
0291:
0292: p = 1.0 / (x * x);
0293: if (x >= 1000.0)
0294: q += ((7.9365079365079365079365e-4 * p - 2.7777777777777777777778e-3)
0295: * p + 0.0833333333333333333333)
0296: / x;
0297: else
0298: q += polevl(p, A, 4) / x;
0299: return q;
0300: }
0301:
0302: /**
0303: * Returns the error function of the normal distribution.
0304: * The integral is
0305: * <pre>
0306: * x
0307: * -
0308: * 2 | | 2
0309: * erf(x) = -------- | exp( - t ) dt.
0310: * sqrt(pi) | |
0311: * -
0312: * 0
0313: * </pre>
0314: * <b>Implementation:</b>
0315: * For <tt>0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2)</tt>; otherwise
0316: * <tt>erf(x) = 1 - erfc(x)</tt>.
0317: * <p>
0318: * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
0319: * Java 2D Graph Package 2.4</A>,
0320: * which in turn is a port from the
0321: * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
0322: * Math Library (C).
0323: *
0324: * @param a the argument to the function.
0325: */
0326: static double errorFunction(double x) {
0327: double y, z;
0328: final double T[] = { 9.60497373987051638749E0,
0329: 9.00260197203842689217E1, 2.23200534594684319226E3,
0330: 7.00332514112805075473E3, 5.55923013010394962768E4 };
0331: final double U[] = {
0332: //1.00000000000000000000E0,
0333: 3.35617141647503099647E1, 5.21357949780152679795E2,
0334: 4.59432382970980127987E3, 2.26290000613890934246E4,
0335: 4.92673942608635921086E4 };
0336:
0337: if (Math.abs(x) > 1.0)
0338: return (1.0 - errorFunctionComplemented(x));
0339: z = x * x;
0340: y = x * polevl(z, T, 4) / p1evl(z, U, 5);
0341: return y;
0342: }
0343:
0344: /**
0345: * Returns the complementary Error function of the normal distribution.
0346: * <pre>
0347: * 1 - erf(x) =
0348: *
0349: * inf.
0350: * -
0351: * 2 | | 2
0352: * erfc(x) = -------- | exp( - t ) dt
0353: * sqrt(pi) | |
0354: * -
0355: * x
0356: * </pre>
0357: * <b>Implementation:</b>
0358: * For small x, <tt>erfc(x) = 1 - erf(x)</tt>; otherwise rational
0359: * approximations are computed.
0360: * <p>
0361: * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
0362: * Java 2D Graph Package 2.4</A>,
0363: * which in turn is a port from the
0364: * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
0365: * Math Library (C).
0366: *
0367: * @param a the argument to the function.
0368: */
0369: static double errorFunctionComplemented(double a) {
0370: double x, y, z, p, q;
0371:
0372: double P[] = { 2.46196981473530512524E-10,
0373: 5.64189564831068821977E-1, 7.46321056442269912687E0,
0374: 4.86371970985681366614E1, 1.96520832956077098242E2,
0375: 5.26445194995477358631E2, 9.34528527171957607540E2,
0376: 1.02755188689515710272E3, 5.57535335369399327526E2 };
0377: double Q[] = {
0378: //1.0
0379: 1.32281951154744992508E1, 8.67072140885989742329E1,
0380: 3.54937778887819891062E2, 9.75708501743205489753E2,
0381: 1.82390916687909736289E3, 2.24633760818710981792E3,
0382: 1.65666309194161350182E3, 5.57535340817727675546E2 };
0383:
0384: double R[] = { 5.64189583547755073984E-1,
0385: 1.27536670759978104416E0, 5.01905042251180477414E0,
0386: 6.16021097993053585195E0, 7.40974269950448939160E0,
0387: 2.97886665372100240670E0 };
0388: double S[] = {
0389: //1.00000000000000000000E0,
0390: 2.26052863220117276590E0, 9.39603524938001434673E0,
0391: 1.20489539808096656605E1, 1.70814450747565897222E1,
0392: 9.60896809063285878198E0, 3.36907645100081516050E0 };
0393:
0394: if (a < 0.0)
0395: x = -a;
0396: else
0397: x = a;
0398:
0399: if (x < 1.0)
0400: return 1.0 - errorFunction(a);
0401:
0402: z = -a * a;
0403:
0404: if (z < -MAXLOG) {
0405: if (a < 0)
0406: return (2.0);
0407: else
0408: return (0.0);
0409: }
0410:
0411: z = Math.exp(z);
0412:
0413: if (x < 8.0) {
0414: p = polevl(x, P, 8);
0415: q = p1evl(x, Q, 8);
0416: } else {
0417: p = polevl(x, R, 5);
0418: q = p1evl(x, S, 6);
0419: }
0420:
0421: y = (z * p) / q;
0422:
0423: if (a < 0)
0424: y = 2.0 - y;
0425:
0426: if (y == 0.0) {
0427: if (a < 0)
0428: return 2.0;
0429: else
0430: return (0.0);
0431: }
0432: return y;
0433: }
0434:
0435: /**
0436: * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
0437: * Evaluates polynomial when coefficient of N is 1.0.
0438: * Otherwise same as <tt>polevl()</tt>.
0439: * <pre>
0440: * 2 N
0441: * y = C + C x + C x +...+ C x
0442: * 0 1 2 N
0443: *
0444: * Coefficients are stored in reverse order:
0445: *
0446: * coef[0] = C , ..., coef[N] = C .
0447: * N 0
0448: * </pre>
0449: * The function <tt>p1evl()</tt> assumes that <tt>coef[N] = 1.0</tt> and is
0450: * omitted from the array. Its calling arguments are
0451: * otherwise the same as <tt>polevl()</tt>.
0452: * <p>
0453: * In the interest of speed, there are no checks for out of bounds arithmetic.
0454: *
0455: * @param x argument to the polynomial.
0456: * @param coef the coefficients of the polynomial.
0457: * @param N the degree of the polynomial.
0458: */
0459: static double p1evl(double x, double coef[], int N) {
0460:
0461: double ans;
0462: ans = x + coef[0];
0463:
0464: for (int i = 1; i < N; i++)
0465: ans = ans * x + coef[i];
0466:
0467: return ans;
0468: }
0469:
0470: /**
0471: * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
0472: * <pre>
0473: * 2 N
0474: * y = C + C x + C x +...+ C x
0475: * 0 1 2 N
0476: *
0477: * Coefficients are stored in reverse order:
0478: *
0479: * coef[0] = C , ..., coef[N] = C .
0480: * N 0
0481: * </pre>
0482: * In the interest of speed, there are no checks for out of bounds arithmetic.
0483: *
0484: * @param x argument to the polynomial.
0485: * @param coef the coefficients of the polynomial.
0486: * @param N the degree of the polynomial.
0487: */
0488: static double polevl(double x, double coef[], int N) {
0489:
0490: double ans;
0491: ans = coef[0];
0492:
0493: for (int i = 1; i <= N; i++)
0494: ans = ans * x + coef[i];
0495:
0496: return ans;
0497: }
0498:
0499: /**
0500: * Returns the Incomplete Gamma function.
0501: * @param a the parameter of the gamma distribution.
0502: * @param x the integration end point.
0503: */
0504: static double incompleteGamma(double a, double x) {
0505:
0506: double ans, ax, c, r;
0507:
0508: if (x <= 0 || a <= 0)
0509: return 0.0;
0510:
0511: if (x > 1.0 && x > a)
0512: return 1.0 - incompleteGammaComplement(a, x);
0513:
0514: /* Compute x**a * exp(-x) / gamma(a) */
0515: ax = a * Math.log(x) - x - lnGamma(a);
0516: if (ax < -MAXLOG)
0517: return (0.0);
0518:
0519: ax = Math.exp(ax);
0520:
0521: /* power series */
0522: r = a;
0523: c = 1.0;
0524: ans = 1.0;
0525:
0526: do {
0527: r += 1.0;
0528: c *= x / r;
0529: ans += c;
0530: } while (c / ans > MACHEP);
0531:
0532: return (ans * ax / a);
0533: }
0534:
0535: /**
0536: * Returns the Complemented Incomplete Gamma function.
0537: * @param a the parameter of the gamma distribution.
0538: * @param x the integration start point.
0539: */
0540: static double incompleteGammaComplement(double a, double x) {
0541:
0542: double ans, ax, c, yc, r, t, y, z;
0543: double pk, pkm1, pkm2, qk, qkm1, qkm2;
0544:
0545: if (x <= 0 || a <= 0)
0546: return 1.0;
0547:
0548: if (x < 1.0 || x < a)
0549: return 1.0 - incompleteGamma(a, x);
0550:
0551: ax = a * Math.log(x) - x - lnGamma(a);
0552: if (ax < -MAXLOG)
0553: return 0.0;
0554:
0555: ax = Math.exp(ax);
0556:
0557: /* continued fraction */
0558: y = 1.0 - a;
0559: z = x + y + 1.0;
0560: c = 0.0;
0561: pkm2 = 1.0;
0562: qkm2 = x;
0563: pkm1 = x + 1.0;
0564: qkm1 = z * x;
0565: ans = pkm1 / qkm1;
0566:
0567: do {
0568: c += 1.0;
0569: y += 1.0;
0570: z += 2.0;
0571: yc = y * c;
0572: pk = pkm1 * z - pkm2 * yc;
0573: qk = qkm1 * z - qkm2 * yc;
0574: if (qk != 0) {
0575: r = pk / qk;
0576: t = Math.abs((ans - r) / r);
0577: ans = r;
0578: } else
0579: t = 1.0;
0580:
0581: pkm2 = pkm1;
0582: pkm1 = pk;
0583: qkm2 = qkm1;
0584: qkm1 = qk;
0585: if (Math.abs(pk) > big) {
0586: pkm2 *= biginv;
0587: pkm1 *= biginv;
0588: qkm2 *= biginv;
0589: qkm1 *= biginv;
0590: }
0591: } while (t > MACHEP);
0592:
0593: return ans * ax;
0594: }
0595:
0596: /**
0597: * Returns the Gamma function of the argument.
0598: */
0599: static double gamma(double x) {
0600:
0601: double P[] = { 1.60119522476751861407E-4,
0602: 1.19135147006586384913E-3, 1.04213797561761569935E-2,
0603: 4.76367800457137231464E-2, 2.07448227648435975150E-1,
0604: 4.94214826801497100753E-1, 9.99999999999999996796E-1 };
0605: double Q[] = { -2.31581873324120129819E-5,
0606: 5.39605580493303397842E-4, -4.45641913851797240494E-3,
0607: 1.18139785222060435552E-2, 3.58236398605498653373E-2,
0608: -2.34591795718243348568E-1, 7.14304917030273074085E-2,
0609: 1.00000000000000000320E0 };
0610:
0611: double p, z;
0612: int i;
0613:
0614: double q = Math.abs(x);
0615:
0616: if (q > 33.0) {
0617: if (x < 0.0) {
0618: p = Math.floor(q);
0619: if (p == q)
0620: throw new ArithmeticException("gamma: overflow");
0621: i = (int) p;
0622: z = q - p;
0623: if (z > 0.5) {
0624: p += 1.0;
0625: z = q - p;
0626: }
0627: z = q * Math.sin(Math.PI * z);
0628: if (z == 0.0)
0629: throw new ArithmeticException("gamma: overflow");
0630: z = Math.abs(z);
0631: z = Math.PI / (z * stirlingFormula(q));
0632:
0633: return -z;
0634: } else {
0635: return stirlingFormula(x);
0636: }
0637: }
0638:
0639: z = 1.0;
0640: while (x >= 3.0) {
0641: x -= 1.0;
0642: z *= x;
0643: }
0644:
0645: while (x < 0.0) {
0646: if (x == 0.0) {
0647: throw new ArithmeticException("gamma: singular");
0648: } else if (x > -1.E-9) {
0649: return (z / ((1.0 + 0.5772156649015329 * x) * x));
0650: }
0651: z /= x;
0652: x += 1.0;
0653: }
0654:
0655: while (x < 2.0) {
0656: if (x == 0.0) {
0657: throw new ArithmeticException("gamma: singular");
0658: } else if (x < 1.e-9) {
0659: return (z / ((1.0 + 0.5772156649015329 * x) * x));
0660: }
0661: z /= x;
0662: x += 1.0;
0663: }
0664:
0665: if ((x == 2.0) || (x == 3.0))
0666: return z;
0667:
0668: x -= 2.0;
0669: p = polevl(x, P, 6);
0670: q = polevl(x, Q, 7);
0671: return z * p / q;
0672: }
0673:
0674: /**
0675: * Returns the Gamma function computed by Stirling's formula.
0676: * The polynomial STIR is valid for 33 <= x <= 172.
0677: */
0678: static double stirlingFormula(double x) {
0679:
0680: double STIR[] = { 7.87311395793093628397E-4,
0681: -2.29549961613378126380E-4, -2.68132617805781232825E-3,
0682: 3.47222221605458667310E-3, 8.33333333333482257126E-2, };
0683: double MAXSTIR = 143.01608;
0684:
0685: double w = 1.0 / x;
0686: double y = Math.exp(x);
0687:
0688: w = 1.0 + w * polevl(w, STIR, 4);
0689:
0690: if (x > MAXSTIR) {
0691: /* Avoid overflow in Math.pow() */
0692: double v = Math.pow(x, 0.5 * x - 0.25);
0693: y = v * (v / y);
0694: } else {
0695: y = Math.pow(x, x - 0.5) / y;
0696: }
0697: y = SQTPI * y * w;
0698: return y;
0699: }
0700:
0701: /**
0702: * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>.
0703: *
0704: * @param aa the alpha parameter of the beta distribution.
0705: * @param bb the beta parameter of the beta distribution.
0706: * @param xx the integration end point.
0707: */
0708: public static double incompleteBeta(double aa, double bb, double xx) {
0709:
0710: double a, b, t, x, xc, w, y;
0711: boolean flag;
0712:
0713: if (aa <= 0.0 || bb <= 0.0)
0714: throw new ArithmeticException("ibeta: Domain error!");
0715:
0716: if ((xx <= 0.0) || (xx >= 1.0)) {
0717: if (xx == 0.0)
0718: return 0.0;
0719: if (xx == 1.0)
0720: return 1.0;
0721: throw new ArithmeticException("ibeta: Domain error!");
0722: }
0723:
0724: flag = false;
0725: if ((bb * xx) <= 1.0 && xx <= 0.95) {
0726: t = powerSeries(aa, bb, xx);
0727: return t;
0728: }
0729:
0730: w = 1.0 - xx;
0731:
0732: /* Reverse a and b if x is greater than the mean. */
0733: if (xx > (aa / (aa + bb))) {
0734: flag = true;
0735: a = bb;
0736: b = aa;
0737: xc = xx;
0738: x = w;
0739: } else {
0740: a = aa;
0741: b = bb;
0742: xc = w;
0743: x = xx;
0744: }
0745:
0746: if (flag && (b * x) <= 1.0 && x <= 0.95) {
0747: t = powerSeries(a, b, x);
0748: if (t <= MACHEP)
0749: t = 1.0 - MACHEP;
0750: else
0751: t = 1.0 - t;
0752: return t;
0753: }
0754:
0755: /* Choose expansion for better convergence. */
0756: y = x * (a + b - 2.0) - (a - 1.0);
0757: if (y < 0.0)
0758: w = incompleteBetaFraction1(a, b, x);
0759: else
0760: w = incompleteBetaFraction2(a, b, x) / xc;
0761:
0762: /* Multiply w by the factor
0763: a b _ _ _
0764: x (1-x) | (a+b) / ( a | (a) | (b) ) . */
0765:
0766: y = a * Math.log(x);
0767: t = b * Math.log(xc);
0768: if ((a + b) < MAXGAM && Math.abs(y) < MAXLOG
0769: && Math.abs(t) < MAXLOG) {
0770: t = Math.pow(xc, b);
0771: t *= Math.pow(x, a);
0772: t /= a;
0773: t *= w;
0774: t *= gamma(a + b) / (gamma(a) * gamma(b));
0775: if (flag) {
0776: if (t <= MACHEP)
0777: t = 1.0 - MACHEP;
0778: else
0779: t = 1.0 - t;
0780: }
0781: return t;
0782: }
0783: /* Resort to logarithms. */
0784: y += t + lnGamma(a + b) - lnGamma(a) - lnGamma(b);
0785: y += Math.log(w / a);
0786: if (y < MINLOG)
0787: t = 0.0;
0788: else
0789: t = Math.exp(y);
0790:
0791: if (flag) {
0792: if (t <= MACHEP)
0793: t = 1.0 - MACHEP;
0794: else
0795: t = 1.0 - t;
0796: }
0797: return t;
0798: }
0799:
0800: /**
0801: * Continued fraction expansion #1 for incomplete beta integral.
0802: */
0803: static double incompleteBetaFraction1(double a, double b, double x) {
0804:
0805: double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
0806: double k1, k2, k3, k4, k5, k6, k7, k8;
0807: double r, t, ans, thresh;
0808: int n;
0809:
0810: k1 = a;
0811: k2 = a + b;
0812: k3 = a;
0813: k4 = a + 1.0;
0814: k5 = 1.0;
0815: k6 = b - 1.0;
0816: k7 = k4;
0817: k8 = a + 2.0;
0818:
0819: pkm2 = 0.0;
0820: qkm2 = 1.0;
0821: pkm1 = 1.0;
0822: qkm1 = 1.0;
0823: ans = 1.0;
0824: r = 1.0;
0825: n = 0;
0826: thresh = 3.0 * MACHEP;
0827: do {
0828: xk = -(x * k1 * k2) / (k3 * k4);
0829: pk = pkm1 + pkm2 * xk;
0830: qk = qkm1 + qkm2 * xk;
0831: pkm2 = pkm1;
0832: pkm1 = pk;
0833: qkm2 = qkm1;
0834: qkm1 = qk;
0835:
0836: xk = (x * k5 * k6) / (k7 * k8);
0837: pk = pkm1 + pkm2 * xk;
0838: qk = qkm1 + qkm2 * xk;
0839: pkm2 = pkm1;
0840: pkm1 = pk;
0841: qkm2 = qkm1;
0842: qkm1 = qk;
0843:
0844: if (qk != 0)
0845: r = pk / qk;
0846: if (r != 0) {
0847: t = Math.abs((ans - r) / r);
0848: ans = r;
0849: } else
0850: t = 1.0;
0851:
0852: if (t < thresh)
0853: return ans;
0854:
0855: k1 += 1.0;
0856: k2 += 1.0;
0857: k3 += 2.0;
0858: k4 += 2.0;
0859: k5 += 1.0;
0860: k6 -= 1.0;
0861: k7 += 2.0;
0862: k8 += 2.0;
0863:
0864: if ((Math.abs(qk) + Math.abs(pk)) > big) {
0865: pkm2 *= biginv;
0866: pkm1 *= biginv;
0867: qkm2 *= biginv;
0868: qkm1 *= biginv;
0869: }
0870: if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
0871: pkm2 *= big;
0872: pkm1 *= big;
0873: qkm2 *= big;
0874: qkm1 *= big;
0875: }
0876: } while (++n < 300);
0877:
0878: return ans;
0879: }
0880:
0881: /**
0882: * Continued fraction expansion #2 for incomplete beta integral.
0883: */
0884: static double incompleteBetaFraction2(double a, double b, double x) {
0885:
0886: double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
0887: double k1, k2, k3, k4, k5, k6, k7, k8;
0888: double r, t, ans, z, thresh;
0889: int n;
0890:
0891: k1 = a;
0892: k2 = b - 1.0;
0893: k3 = a;
0894: k4 = a + 1.0;
0895: k5 = 1.0;
0896: k6 = a + b;
0897: k7 = a + 1.0;
0898: ;
0899: k8 = a + 2.0;
0900:
0901: pkm2 = 0.0;
0902: qkm2 = 1.0;
0903: pkm1 = 1.0;
0904: qkm1 = 1.0;
0905: z = x / (1.0 - x);
0906: ans = 1.0;
0907: r = 1.0;
0908: n = 0;
0909: thresh = 3.0 * MACHEP;
0910: do {
0911: xk = -(z * k1 * k2) / (k3 * k4);
0912: pk = pkm1 + pkm2 * xk;
0913: qk = qkm1 + qkm2 * xk;
0914: pkm2 = pkm1;
0915: pkm1 = pk;
0916: qkm2 = qkm1;
0917: qkm1 = qk;
0918:
0919: xk = (z * k5 * k6) / (k7 * k8);
0920: pk = pkm1 + pkm2 * xk;
0921: qk = qkm1 + qkm2 * xk;
0922: pkm2 = pkm1;
0923: pkm1 = pk;
0924: qkm2 = qkm1;
0925: qkm1 = qk;
0926:
0927: if (qk != 0)
0928: r = pk / qk;
0929: if (r != 0) {
0930: t = Math.abs((ans - r) / r);
0931: ans = r;
0932: } else
0933: t = 1.0;
0934:
0935: if (t < thresh)
0936: return ans;
0937:
0938: k1 += 1.0;
0939: k2 -= 1.0;
0940: k3 += 2.0;
0941: k4 += 2.0;
0942: k5 += 1.0;
0943: k6 += 1.0;
0944: k7 += 2.0;
0945: k8 += 2.0;
0946:
0947: if ((Math.abs(qk) + Math.abs(pk)) > big) {
0948: pkm2 *= biginv;
0949: pkm1 *= biginv;
0950: qkm2 *= biginv;
0951: qkm1 *= biginv;
0952: }
0953: if ((Math.abs(qk) < biginv) || (Math.abs(pk) < biginv)) {
0954: pkm2 *= big;
0955: pkm1 *= big;
0956: qkm2 *= big;
0957: qkm1 *= big;
0958: }
0959: } while (++n < 300);
0960:
0961: return ans;
0962: }
0963:
0964: /**
0965: * Power series for incomplete beta integral.
0966: * Use when b*x is small and x not too close to 1.
0967: */
0968: static double powerSeries(double a, double b, double x) {
0969:
0970: double s, t, u, v, n, t1, z, ai;
0971:
0972: ai = 1.0 / a;
0973: u = (1.0 - b) * x;
0974: v = u / (a + 1.0);
0975: t1 = v;
0976: t = u;
0977: n = 2.0;
0978: s = 0.0;
0979: z = MACHEP * ai;
0980: while (Math.abs(v) > z) {
0981: u = (n - b) * x / n;
0982: t *= u;
0983: v = t / (a + n);
0984: s += v;
0985: n += 1.0;
0986: }
0987: s += t1;
0988: s += ai;
0989:
0990: u = a * Math.log(x);
0991: if ((a + b) < MAXGAM && Math.abs(u) < MAXLOG) {
0992: t = gamma(a + b) / (gamma(a) * gamma(b));
0993: s = s * t * Math.pow(x, a);
0994: } else {
0995: t = lnGamma(a + b) - lnGamma(a) - lnGamma(b) + u
0996: + Math.log(s);
0997: if (t < MINLOG)
0998: s = 0.0;
0999: else
1000: s = Math.exp(t);
1001: }
1002: return s;
1003: }
1004:
1005: /**
1006: * Main method for testing this class.
1007: */
1008: public static void main(String[] ops) {
1009:
1010: System.out.println("Binomial standard error (0.5, 100): "
1011: + Statistics.binomialStandardError(0.5, 100));
1012: System.out.println("Chi-squared probability (2.558, 10): "
1013: + Statistics.chiSquaredProbability(2.558, 10));
1014: System.out.println("Normal probability (0.2): "
1015: + Statistics.normalProbability(0.2));
1016: System.out.println("F probability (5.1922, 4, 5): "
1017: + Statistics.FProbability(5.1922, 4, 5));
1018: System.out.println("lnGamma(6): " + Statistics.lnGamma(6));
1019: }
1020: }
|