source: src/main/java/weka/core/Statistics.java @ 6

Last change on this file since 6 was 4, checked in by gnappo, 14 years ago

Import di weka.
File size: 26.7 KB
Line 
1package weka.core;
2
3/**
4 * Class implementing some distributions, tests, etc. The code is mostly adapted from the CERN
5 * Jet Java libraries:
6 *
7 * Copyright 2001 University of Waikato
8 * Copyright 1999 CERN - European Organization for Nuclear Research.
9 * Permission to use, copy, modify, distribute and sell this software and its documentation for
10 * any purpose is hereby granted without fee, provided that the above copyright notice appear
11 * in all copies and that both that copyright notice and this permission notice appear in
12 * supporting documentation.
13 * CERN and the University of Waikato make no representations about the suitability of this
14 * software for any purpose. It is provided "as is" without expressed or implied warranty.
15 *
16 * @author peter.gedeck@pharma.Novartis.com
17 * @author wolfgang.hoschek@cern.ch
18 * @author Eibe Frank (eibe@cs.waikato.ac.nz)
19 * @author Richard Kirkby (rkirkby@cs.waikato.ac.nz)
20 * @version $Revision: 5616 $
21 */
22public class Statistics
23  implements RevisionHandler {
24
25  /** Some constants */
26  protected static final double MACHEP =  1.11022302462515654042E-16;
27  protected static final double MAXLOG =  7.09782712893383996732E2;
28  protected static final double MINLOG = -7.451332191019412076235E2;
29  protected static final double MAXGAM = 171.624376956302725;
30  protected static final double SQTPI  =  2.50662827463100050242E0;
31  protected static final double SQRTH  =  7.07106781186547524401E-1;
32  protected static final double LOGPI  =  1.14472988584940017414;
33 
34  protected static final double big    =  4.503599627370496e15;
35  protected static final double biginv =  2.22044604925031308085e-16;
36
37  /*************************************************
38   *    COEFFICIENTS FOR METHOD  normalInverse()   *
39   *************************************************/
40  /* approximation for 0 <= |y - 0.5| <= 3/8 */
41  protected static final double P0[] = {
42    -5.99633501014107895267E1,
43    9.80010754185999661536E1,
44    -5.66762857469070293439E1,
45    1.39312609387279679503E1,
46    -1.23916583867381258016E0,
47  };
48  protected static final double Q0[] = {
49    /* 1.00000000000000000000E0,*/
50    1.95448858338141759834E0,
51    4.67627912898881538453E0,
52    8.63602421390890590575E1,
53    -2.25462687854119370527E2,
54    2.00260212380060660359E2,
55    -8.20372256168333339912E1,
56    1.59056225126211695515E1,
57    -1.18331621121330003142E0,
58  };
59 
60  /* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
61   * i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
62   */
63  protected static final double P1[] = {
64    4.05544892305962419923E0,
65    3.15251094599893866154E1,
66    5.71628192246421288162E1,
67    4.40805073893200834700E1,
68    1.46849561928858024014E1,
69    2.18663306850790267539E0,
70    -1.40256079171354495875E-1,
71    -3.50424626827848203418E-2,
72    -8.57456785154685413611E-4,
73  };
74  protected static final double Q1[] = {
75    /*  1.00000000000000000000E0,*/
76    1.57799883256466749731E1,
77    4.53907635128879210584E1,
78    4.13172038254672030440E1,
79    1.50425385692907503408E1,
80    2.50464946208309415979E0,
81    -1.42182922854787788574E-1,
82    -3.80806407691578277194E-2,
83    -9.33259480895457427372E-4,
84  };
85 
86  /* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
87   * i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
88   */
89  protected static final double  P2[] = {
90    3.23774891776946035970E0,
91    6.91522889068984211695E0,
92    3.93881025292474443415E0,
93    1.33303460815807542389E0,
94    2.01485389549179081538E-1,
95    1.23716634817820021358E-2,
96    3.01581553508235416007E-4,
97    2.65806974686737550832E-6,
98    6.23974539184983293730E-9,
99  };
100  protected static final double  Q2[] = {
101    /*  1.00000000000000000000E0,*/
102    6.02427039364742014255E0,
103    3.67983563856160859403E0,
104    1.37702099489081330271E0,
105    2.16236993594496635890E-1,
106    1.34204006088543189037E-2,
107    3.28014464682127739104E-4,
108    2.89247864745380683936E-6,
109    6.79019408009981274425E-9,
110  };
111 
112  /**
113   * Computes standard error for observed values of a binomial
114   * random variable.
115   *
116   * @param p the probability of success
117   * @param n the size of the sample
118   * @return the standard error
119   */
120  public static double binomialStandardError(double p, int n) {
121   
122    if (n == 0) {
123      return 0; 
124    }
125    return Math.sqrt((p*(1-p))/(double) n);
126  }
127 
128  /**
129   * Returns chi-squared probability for given value and degrees
130   * of freedom. (The probability that the chi-squared variate
131   * will be greater than x for the given degrees of freedom.)
132   *
133   * @param x the value
134   * @param v the number of degrees of freedom
135   * @return the chi-squared probability
136   */
137  public static double chiSquaredProbability(double x, double v) { 
138
139    if( x < 0.0 || v < 1.0 ) return 0.0;
140    return incompleteGammaComplement( v/2.0, x/2.0 );
141  }
142
143  /**
144   * Computes probability of F-ratio.
145   *
146   * @param F the F-ratio
147   * @param df1 the first number of degrees of freedom
148   * @param df2 the second number of degrees of freedom
149   * @return the probability of the F-ratio.
150   */
151  public static double FProbability(double F, int df1, int df2) {
152   
153    return incompleteBeta( df2/2.0, df1/2.0, df2/(df2+df1*F) );
154  }
155
156  /**
157   * Returns the area under the Normal (Gaussian) probability density
158   * function, integrated from minus infinity to <tt>x</tt>
159   * (assumes mean is zero, variance is one).
160   * <pre>
161   *                            x
162   *                             -
163   *                   1        | |          2
164   *  normal(x)  = ---------    |    exp( - t /2 ) dt
165   *               sqrt(2pi)  | |
166   *                           -
167   *                          -inf.
168   *
169   *             =  ( 1 + erf(z) ) / 2
170   *             =  erfc(z) / 2
171   * </pre>
172   * where <tt>z = x/sqrt(2)</tt>.
173   * Computation is via the functions <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
174   *
175   * @param a the z-value
176   * @return the probability of the z value according to the normal pdf
177   */
178  public static double normalProbability(double a) { 
179
180    double x, y, z;
181 
182    x = a * SQRTH;
183    z = Math.abs(x);
184 
185    if( z < SQRTH ) y = 0.5 + 0.5 * errorFunction(x);
186    else {
187      y = 0.5 * errorFunctionComplemented(z);
188      if( x > 0 )  y = 1.0 - y;
189    } 
190    return y;
191  }
192
193  /**
194   * Returns the value, <tt>x</tt>, for which the area under the
195   * Normal (Gaussian) probability density function (integrated from
196   * minus infinity to <tt>x</tt>) is equal to the argument <tt>y</tt>
197   * (assumes mean is zero, variance is one).
198   * <p>
199   * For small arguments <tt>0 < y < exp(-2)</tt>, the program computes
200   * <tt>z = sqrt( -2.0 * log(y) )</tt>;  then the approximation is
201   * <tt>x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z)</tt>.
202   * There are two rational functions P/Q, one for <tt>0 < y < exp(-32)</tt>
203   * and the other for <tt>y</tt> up to <tt>exp(-2)</tt>.
204   * For larger arguments,
205   * <tt>w = y - 0.5</tt>, and  <tt>x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2))</tt>.
206   *
207   * @param y0 the area under the normal pdf
208   * @return the z-value
209   */
210  public static double normalInverse(double y0) { 
211
212    double x, y, z, y2, x0, x1;
213    int code;
214
215    final double s2pi = Math.sqrt(2.0*Math.PI);
216
217    if( y0 <= 0.0 ) throw new IllegalArgumentException();
218    if( y0 >= 1.0 ) throw new IllegalArgumentException();
219    code = 1;
220    y = y0;
221    if( y > (1.0 - 0.13533528323661269189) ) { /* 0.135... = exp(-2) */
222      y = 1.0 - y;
223      code = 0;
224    }
225
226    if( y > 0.13533528323661269189 ) {
227      y = y - 0.5;
228      y2 = y * y;
229      x = y + y * (y2 * polevl( y2, P0, 4)/p1evl( y2, Q0, 8 ));
230      x = x * s2pi; 
231      return(x);
232    }
233
234    x = Math.sqrt( -2.0 * Math.log(y) );
235    x0 = x - Math.log(x)/x;
236
237    z = 1.0/x;
238    if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
239      x1 = z * polevl( z, P1, 8 )/p1evl( z, Q1, 8 );
240    else
241      x1 = z * polevl( z, P2, 8 )/p1evl( z, Q2, 8 );
242    x = x0 - x1;
243    if( code != 0 )
244      x = -x;
245    return( x );
246  }
247
248  /**
249   * Returns natural logarithm of gamma function.
250   *
251   * @param x the value
252   * @return natural logarithm of gamma function
253   */
254  public static double lnGamma(double x) {
255
256    double p, q, w, z;
257 
258    double A[] = {
259      8.11614167470508450300E-4,
260      -5.95061904284301438324E-4,
261      7.93650340457716943945E-4,
262      -2.77777777730099687205E-3,
263      8.33333333333331927722E-2
264    };
265    double B[] = {
266      -1.37825152569120859100E3,
267      -3.88016315134637840924E4,
268      -3.31612992738871184744E5,
269      -1.16237097492762307383E6,
270      -1.72173700820839662146E6,
271      -8.53555664245765465627E5
272    };
273    double C[] = {
274      /* 1.00000000000000000000E0, */
275      -3.51815701436523470549E2,
276      -1.70642106651881159223E4,
277      -2.20528590553854454839E5,
278      -1.13933444367982507207E6,
279      -2.53252307177582951285E6,
280      -2.01889141433532773231E6
281    };
282 
283    if( x < -34.0 ) {
284      q = -x;
285      w = lnGamma(q);
286      p = Math.floor(q);
287      if( p == q ) throw new ArithmeticException("lnGamma: Overflow");
288      z = q - p;
289      if( z > 0.5 ) {
290        p += 1.0;
291        z = p - q;
292      }
293      z = q * Math.sin( Math.PI * z );
294      if( z == 0.0 ) throw new 
295        ArithmeticException("lnGamma: Overflow");
296      z = LOGPI - Math.log( z ) - w;
297      return z;
298    }
299 
300    if( x < 13.0 ) {
301      z = 1.0;
302      while( x >= 3.0 ) {
303        x -= 1.0;
304        z *= x;
305      }
306      while( x < 2.0 ) {
307        if( x == 0.0 ) throw new 
308          ArithmeticException("lnGamma: Overflow");
309        z /= x;
310        x += 1.0;
311      }
312      if( z < 0.0 ) z = -z;
313      if( x == 2.0 ) return Math.log(z);
314      x -= 2.0;
315      p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
316      return( Math.log(z) + p );
317    }
318 
319    if( x > 2.556348e305 ) throw new ArithmeticException("lnGamma: Overflow");
320 
321    q = ( x - 0.5 ) * Math.log(x) - x + 0.91893853320467274178;
322 
323    if( x > 1.0e8 ) return( q );
324 
325    p = 1.0/(x*x);
326    if( x >= 1000.0 )
327      q += ((   7.9365079365079365079365e-4 * p
328                - 2.7777777777777777777778e-3) *p
329            + 0.0833333333333333333333) / x;
330    else
331      q += polevl( p, A, 4 ) / x;
332    return q;
333  }
334
335  /**
336   * Returns the error function of the normal distribution.
337   * The integral is
338   * <pre>
339   *                           x
340   *                            -
341   *                 2         | |          2
342   *   erf(x)  =  --------     |    exp( - t  ) dt.
343   *              sqrt(pi)   | |
344   *                          -
345   *                           0
346   * </pre>
347   * <b>Implementation:</b>
348   * For <tt>0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2)</tt>; otherwise
349   * <tt>erf(x) = 1 - erfc(x)</tt>.
350   * <p>
351   * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
352   * Java 2D Graph Package 2.4</A>,
353   * which in turn is a port from the
354   * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
355   * Math Library (C).
356   *
357   * @param a the argument to the function.
358   */
359  public static double errorFunction(double x) { 
360    double y, z;
361    final double T[] = {
362      9.60497373987051638749E0,
363      9.00260197203842689217E1,
364      2.23200534594684319226E3,
365      7.00332514112805075473E3,
366      5.55923013010394962768E4
367    };
368    final double U[] = {
369      //1.00000000000000000000E0,
370      3.35617141647503099647E1,
371      5.21357949780152679795E2,
372      4.59432382970980127987E3,
373      2.26290000613890934246E4,
374      4.92673942608635921086E4
375    };
376 
377    if( Math.abs(x) > 1.0 ) return( 1.0 - errorFunctionComplemented(x) );
378    z = x * x;
379    y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
380    return y;
381  }
382
383  /**
384   * Returns the complementary Error function of the normal distribution.
385   * <pre>
386   *  1 - erf(x) =
387   *
388   *                           inf.
389   *                             -
390   *                  2         | |          2
391   *   erfc(x)  =  --------     |    exp( - t  ) dt
392   *               sqrt(pi)   | |
393   *                           -
394   *                            x
395   * </pre>
396   * <b>Implementation:</b>
397   * For small x, <tt>erfc(x) = 1 - erf(x)</tt>; otherwise rational
398   * approximations are computed.
399   * <p>
400   * Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">
401   * Java 2D Graph Package 2.4</A>,
402   * which in turn is a port from the
403   * <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A>
404   * Math Library (C).
405   *
406   * @param a the argument to the function.
407   */
408  public static double errorFunctionComplemented(double a) { 
409    double x,y,z,p,q;
410 
411    double P[] = {
412      2.46196981473530512524E-10,
413      5.64189564831068821977E-1,
414      7.46321056442269912687E0,
415      4.86371970985681366614E1,
416      1.96520832956077098242E2,
417      5.26445194995477358631E2,
418      9.34528527171957607540E2,
419      1.02755188689515710272E3,
420      5.57535335369399327526E2
421    };
422    double Q[] = {
423      //1.0
424      1.32281951154744992508E1,
425      8.67072140885989742329E1,
426      3.54937778887819891062E2,
427      9.75708501743205489753E2,
428      1.82390916687909736289E3,
429      2.24633760818710981792E3,
430      1.65666309194161350182E3,
431      5.57535340817727675546E2
432    };
433 
434    double R[] = {
435      5.64189583547755073984E-1,
436      1.27536670759978104416E0,
437      5.01905042251180477414E0,
438      6.16021097993053585195E0,
439      7.40974269950448939160E0,
440      2.97886665372100240670E0
441    };
442    double S[] = {
443      //1.00000000000000000000E0,
444      2.26052863220117276590E0,
445      9.39603524938001434673E0,
446      1.20489539808096656605E1,
447      1.70814450747565897222E1,
448      9.60896809063285878198E0,
449      3.36907645100081516050E0
450    };
451 
452    if( a < 0.0 )   x = -a;
453    else            x = a;
454 
455    if( x < 1.0 )   return 1.0 - errorFunction(a);
456 
457    z = -a * a;
458 
459    if( z < -MAXLOG ) {
460      if( a < 0 )  return( 2.0 );
461      else         return( 0.0 );
462    }
463 
464    z = Math.exp(z);
465 
466    if( x < 8.0 ) {
467      p = polevl( x, P, 8 );
468      q = p1evl( x, Q, 8 );
469    } else {
470      p = polevl( x, R, 5 );
471      q = p1evl( x, S, 6 );
472    }
473 
474    y = (z * p)/q;
475 
476    if( a < 0 ) y = 2.0 - y;
477 
478    if( y == 0.0 ) {
479      if( a < 0 ) return 2.0;
480      else        return( 0.0 );
481    }
482    return y;
483  }
484 
485  /**
486   * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
487   * Evaluates polynomial when coefficient of N is 1.0.
488   * Otherwise same as <tt>polevl()</tt>.
489   * <pre>
490   *                     2          N
491   * y  =  C  + C x + C x  +...+ C x
492   *        0    1     2          N
493   *
494   * Coefficients are stored in reverse order:
495   *
496   * coef[0] = C  , ..., coef[N] = C  .
497   *            N                   0
498   * </pre>
499   * The function <tt>p1evl()</tt> assumes that <tt>coef[N] = 1.0</tt> and is
500   * omitted from the array.  Its calling arguments are
501   * otherwise the same as <tt>polevl()</tt>.
502   * <p>
503   * In the interest of speed, there are no checks for out of bounds arithmetic.
504   *
505   * @param x argument to the polynomial.
506   * @param coef the coefficients of the polynomial.
507   * @param N the degree of the polynomial.
508   */
509  public static double p1evl( double x, double coef[], int N ) {
510 
511    double ans;
512    ans = x + coef[0];
513 
514    for(int i=1; i<N; i++) ans = ans*x+coef[i];
515 
516    return ans;
517  }
518
519  /**
520   * Evaluates the given polynomial of degree <tt>N</tt> at <tt>x</tt>.
521   * <pre>
522   *                     2          N
523   * y  =  C  + C x + C x  +...+ C x
524   *        0    1     2          N
525   *
526   * Coefficients are stored in reverse order:
527   *
528   * coef[0] = C  , ..., coef[N] = C  .
529   *            N                   0
530   * </pre>
531   * In the interest of speed, there are no checks for out of bounds arithmetic.
532   *
533   * @param x argument to the polynomial.
534   * @param coef the coefficients of the polynomial.
535   * @param N the degree of the polynomial.
536   */
537  public static double polevl( double x, double coef[], int N ) {
538
539    double ans;
540    ans = coef[0];
541 
542    for(int i=1; i<=N; i++) ans = ans*x+coef[i];
543 
544    return ans;
545  }
546
547  /**
548   * Returns the Incomplete Gamma function.
549   * @param a the parameter of the gamma distribution.
550   * @param x the integration end point.
551   */
552  public static double incompleteGamma(double a, double x) 
553    { 
554 
555    double ans, ax, c, r;
556 
557    if( x <= 0 || a <= 0 ) return 0.0;
558 
559    if( x > 1.0 && x > a ) return 1.0 - incompleteGammaComplement(a,x);
560
561    /* Compute  x**a * exp(-x) / gamma(a)  */
562    ax = a * Math.log(x) - x - lnGamma(a);
563    if( ax < -MAXLOG ) return( 0.0 );
564
565    ax = Math.exp(ax);
566
567    /* power series */
568    r = a;
569    c = 1.0;
570    ans = 1.0;
571
572    do {
573      r += 1.0;
574      c *= x/r;
575      ans += c;
576    }
577    while( c/ans > MACHEP );
578 
579    return( ans * ax/a );
580  }
581
582  /**
583   * Returns the Complemented Incomplete Gamma function.
584   * @param a the parameter of the gamma distribution.
585   * @param x the integration start point.
586   */
587  public static double incompleteGammaComplement( double a, double x ) {
588
589    double ans, ax, c, yc, r, t, y, z;
590    double pk, pkm1, pkm2, qk, qkm1, qkm2;
591
592    if( x <= 0 || a <= 0 ) return 1.0;
593 
594    if( x < 1.0 || x < a ) return 1.0 - incompleteGamma(a,x);
595 
596    ax = a * Math.log(x) - x - lnGamma(a);
597    if( ax < -MAXLOG ) return 0.0;
598 
599    ax = Math.exp(ax);
600 
601    /* continued fraction */
602    y = 1.0 - a;
603    z = x + y + 1.0;
604    c = 0.0;
605    pkm2 = 1.0;
606    qkm2 = x;
607    pkm1 = x + 1.0;
608    qkm1 = z * x;
609    ans = pkm1/qkm1;
610 
611    do {
612      c += 1.0;
613      y += 1.0;
614      z += 2.0;
615      yc = y * c;
616      pk = pkm1 * z  -  pkm2 * yc;
617      qk = qkm1 * z  -  qkm2 * yc;
618      if( qk != 0 ) {
619        r = pk/qk;
620        t = Math.abs( (ans - r)/r );
621        ans = r;
622      } else
623        t = 1.0;
624
625      pkm2 = pkm1;
626      pkm1 = pk;
627      qkm2 = qkm1;
628      qkm1 = qk;
629      if( Math.abs(pk) > big ) {
630        pkm2 *= biginv;
631        pkm1 *= biginv;
632        qkm2 *= biginv;
633        qkm1 *= biginv;
634      }
635    } while( t > MACHEP );
636 
637    return ans * ax;
638  }
639
640  /**
641   * Returns the Gamma function of the argument.
642   */
643  public static double gamma(double x) {
644
645    double P[] = {
646      1.60119522476751861407E-4,
647      1.19135147006586384913E-3,
648      1.04213797561761569935E-2,
649      4.76367800457137231464E-2,
650      2.07448227648435975150E-1,
651      4.94214826801497100753E-1,
652      9.99999999999999996796E-1
653    };
654    double Q[] = {
655      -2.31581873324120129819E-5,
656      5.39605580493303397842E-4,
657      -4.45641913851797240494E-3,
658      1.18139785222060435552E-2,
659      3.58236398605498653373E-2,
660      -2.34591795718243348568E-1,
661      7.14304917030273074085E-2,
662      1.00000000000000000320E0
663    };
664
665    double p, z;
666    int i;
667
668    double q = Math.abs(x);
669
670    if( q > 33.0 ) {
671      if( x < 0.0 ) {
672        p = Math.floor(q);
673        if( p == q ) throw new ArithmeticException("gamma: overflow");
674        i = (int)p;
675        z = q - p;
676        if( z > 0.5 ) {
677          p += 1.0;
678          z = q - p;
679        }
680        z = q * Math.sin( Math.PI * z );
681        if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");
682        z = Math.abs(z);
683        z = Math.PI/(z * stirlingFormula(q) );
684
685        return -z;
686      } else {
687        return stirlingFormula(x);
688      }
689    }
690
691    z = 1.0;
692    while( x >= 3.0 ) {
693      x -= 1.0;
694      z *= x;
695    }
696
697    while( x < 0.0 ) {
698      if( x == 0.0 ) {
699        throw new ArithmeticException("gamma: singular");
700      } else
701        if( x > -1.E-9 ) {
702          return( z/((1.0 + 0.5772156649015329 * x) * x) );
703        }
704      z /= x;
705      x += 1.0;
706    }
707
708    while( x < 2.0 ) {
709      if( x == 0.0 ) {
710        throw new ArithmeticException("gamma: singular");
711      } else
712        if( x < 1.e-9 ) {
713          return( z/((1.0 + 0.5772156649015329 * x) * x) );
714        }
715      z /= x;
716      x += 1.0;
717    }
718
719    if( (x == 2.0) || (x == 3.0) )      return z;
720
721    x -= 2.0;
722    p = polevl( x, P, 6 );
723    q = polevl( x, Q, 7 );
724    return  z * p / q;
725  }
726
727  /**
728   * Returns the Gamma function computed by Stirling's formula.
729   * The polynomial STIR is valid for 33 <= x <= 172.
730   */
731  public static double stirlingFormula(double x) {
732
733    double STIR[] = {
734      7.87311395793093628397E-4,
735      -2.29549961613378126380E-4,
736      -2.68132617805781232825E-3,
737      3.47222221605458667310E-3,
738      8.33333333333482257126E-2,
739    };
740    double MAXSTIR = 143.01608;
741
742    double w = 1.0/x;
743    double  y = Math.exp(x);
744
745    w = 1.0 + w * polevl( w, STIR, 4 );
746
747    if( x > MAXSTIR ) {
748      /* Avoid overflow in Math.pow() */
749      double v = Math.pow( x, 0.5 * x - 0.25 );
750      y = v * (v / y);
751    } else {
752      y = Math.pow( x, x - 0.5 ) / y;
753    }
754    y = SQTPI * y * w;
755    return y;
756  }
757
758  /**
759   * Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>.
760   *
761   * @param aa the alpha parameter of the beta distribution.
762   * @param bb the beta parameter of the beta distribution.
763   * @param xx the integration end point.
764   */
765  public static double incompleteBeta( double aa, double bb, double xx ) {
766
767    double a, b, t, x, xc, w, y;
768    boolean flag;
769
770    if( aa <= 0.0 || bb <= 0.0 ) throw new 
771      ArithmeticException("ibeta: Domain error!");
772
773    if( (xx <= 0.0) || ( xx >= 1.0) ) {
774      if( xx == 0.0 ) return 0.0;
775      if( xx == 1.0 ) return 1.0;
776      throw new ArithmeticException("ibeta: Domain error!");
777    }
778
779    flag = false;
780    if( (bb * xx) <= 1.0 && xx <= 0.95) {
781      t = powerSeries(aa, bb, xx);
782      return t;
783    }
784
785    w = 1.0 - xx;
786
787    /* Reverse a and b if x is greater than the mean. */
788    if( xx > (aa/(aa+bb)) ) {
789      flag = true;
790      a = bb;
791      b = aa;
792      xc = xx;
793      x = w;
794    } else {
795      a = aa;
796      b = bb;
797      xc = w;
798      x = xx;
799    }
800
801    if( flag  && (b * x) <= 1.0 && x <= 0.95) {
802      t = powerSeries(a, b, x);
803      if( t <= MACHEP )         t = 1.0 - MACHEP;
804      else                      t = 1.0 - t;
805      return t;
806    }
807
808    /* Choose expansion for better convergence. */
809    y = x * (a+b-2.0) - (a-1.0);
810    if( y < 0.0 )
811      w = incompleteBetaFraction1( a, b, x );
812    else
813      w = incompleteBetaFraction2( a, b, x ) / xc;
814
815    /* Multiply w by the factor
816       a      b   _             _     _
817       x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */
818
819    y = a * Math.log(x);
820    t = b * Math.log(xc);
821    if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {
822      t = Math.pow(xc,b);
823      t *= Math.pow(x,a);
824      t /= a;
825      t *= w;
826      t *= gamma(a+b) / (gamma(a) * gamma(b));
827      if( flag ) {
828        if( t <= MACHEP )       t = 1.0 - MACHEP;
829        else                    t = 1.0 - t;
830      }
831      return t;
832    }
833    /* Resort to logarithms.  */
834    y += t + lnGamma(a+b) - lnGamma(a) - lnGamma(b);
835    y += Math.log(w/a);
836    if( y < MINLOG )
837      t = 0.0;
838    else
839      t = Math.exp(y);
840
841    if( flag ) {
842      if( t <= MACHEP )         t = 1.0 - MACHEP;
843      else                      t = 1.0 - t;
844    }
845    return t;
846  }   
847
848  /**
849   * Continued fraction expansion #1 for incomplete beta integral.
850   */
851  public static double incompleteBetaFraction1( double a, double b, double x ) {
852
853    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
854    double k1, k2, k3, k4, k5, k6, k7, k8;
855    double r, t, ans, thresh;
856    int n;
857
858    k1 = a;
859    k2 = a + b;
860    k3 = a;
861    k4 = a + 1.0;
862    k5 = 1.0;
863    k6 = b - 1.0;
864    k7 = k4;
865    k8 = a + 2.0;
866
867    pkm2 = 0.0;
868    qkm2 = 1.0;
869    pkm1 = 1.0;
870    qkm1 = 1.0;
871    ans = 1.0;
872    r = 1.0;
873    n = 0;
874    thresh = 3.0 * MACHEP;
875    do {
876      xk = -( x * k1 * k2 )/( k3 * k4 );
877      pk = pkm1 +  pkm2 * xk;
878      qk = qkm1 +  qkm2 * xk;
879      pkm2 = pkm1;
880      pkm1 = pk;
881      qkm2 = qkm1;
882      qkm1 = qk;
883
884      xk = ( x * k5 * k6 )/( k7 * k8 );
885      pk = pkm1 +  pkm2 * xk;
886      qk = qkm1 +  qkm2 * xk;
887      pkm2 = pkm1;
888      pkm1 = pk;
889      qkm2 = qkm1;
890      qkm1 = qk;
891
892      if( qk != 0 )             r = pk/qk;
893      if( r != 0 ) {
894        t = Math.abs( (ans - r)/r );
895        ans = r;
896      } else
897        t = 1.0;
898
899      if( t < thresh ) return ans;
900
901      k1 += 1.0;
902      k2 += 1.0;
903      k3 += 2.0;
904      k4 += 2.0;
905      k5 += 1.0;
906      k6 -= 1.0;
907      k7 += 2.0;
908      k8 += 2.0;
909
910      if( (Math.abs(qk) + Math.abs(pk)) > big ) {
911        pkm2 *= biginv;
912        pkm1 *= biginv;
913        qkm2 *= biginv;
914        qkm1 *= biginv;
915      }
916      if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
917        pkm2 *= big;
918        pkm1 *= big;
919        qkm2 *= big;
920        qkm1 *= big;
921      }
922    } while( ++n < 300 );
923
924    return ans;
925  }   
926
927  /**
928   * Continued fraction expansion #2 for incomplete beta integral.
929   */
930  public static double incompleteBetaFraction2( double a, double b, double x ) {
931
932    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
933    double k1, k2, k3, k4, k5, k6, k7, k8;
934    double r, t, ans, z, thresh;
935    int n;
936
937    k1 = a;
938    k2 = b - 1.0;
939    k3 = a;
940    k4 = a + 1.0;
941    k5 = 1.0;
942    k6 = a + b;
943    k7 = a + 1.0;;
944    k8 = a + 2.0;
945
946    pkm2 = 0.0;
947    qkm2 = 1.0;
948    pkm1 = 1.0;
949    qkm1 = 1.0;
950    z = x / (1.0-x);
951    ans = 1.0;
952    r = 1.0;
953    n = 0;
954    thresh = 3.0 * MACHEP;
955    do {
956      xk = -( z * k1 * k2 )/( k3 * k4 );
957      pk = pkm1 +  pkm2 * xk;
958      qk = qkm1 +  qkm2 * xk;
959      pkm2 = pkm1;
960      pkm1 = pk;
961      qkm2 = qkm1;
962      qkm1 = qk;
963
964      xk = ( z * k5 * k6 )/( k7 * k8 );
965      pk = pkm1 +  pkm2 * xk;
966      qk = qkm1 +  qkm2 * xk;
967      pkm2 = pkm1;
968      pkm1 = pk;
969      qkm2 = qkm1;
970      qkm1 = qk;
971
972      if( qk != 0 )  r = pk/qk;
973      if( r != 0 ) {
974        t = Math.abs( (ans - r)/r );
975        ans = r;
976      } else
977        t = 1.0;
978
979      if( t < thresh ) return ans;
980
981      k1 += 1.0;
982      k2 -= 1.0;
983      k3 += 2.0;
984      k4 += 2.0;
985      k5 += 1.0;
986      k6 += 1.0;
987      k7 += 2.0;
988      k8 += 2.0;
989
990      if( (Math.abs(qk) + Math.abs(pk)) > big ) {
991        pkm2 *= biginv;
992        pkm1 *= biginv;
993        qkm2 *= biginv;
994        qkm1 *= biginv;
995      }
996      if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
997        pkm2 *= big;
998        pkm1 *= big;
999        qkm2 *= big;
1000        qkm1 *= big;
1001      }
1002    } while( ++n < 300 );
1003
1004    return ans;
1005  }
1006
1007  /**
1008   * Power series for incomplete beta integral.
1009   * Use when b*x is small and x not too close to 1. 
1010   */
1011  public static double powerSeries( double a, double b, double x ) {
1012
1013    double s, t, u, v, n, t1, z, ai;
1014   
1015    ai = 1.0 / a;
1016    u = (1.0 - b) * x;
1017    v = u / (a + 1.0);
1018    t1 = v;
1019    t = u;
1020    n = 2.0;
1021    s = 0.0;
1022    z = MACHEP * ai;
1023    while( Math.abs(v) > z ) {
1024      u = (n - b) * x / n;
1025      t *= u;
1026      v = t / (a + n);
1027      s += v; 
1028      n += 1.0;
1029    }
1030    s += t1;
1031    s += ai;
1032
1033    u = a * Math.log(x);
1034    if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {
1035      t = gamma(a+b)/(gamma(a)*gamma(b));
1036      s = s * t * Math.pow(x,a);
1037    } else {
1038      t = lnGamma(a+b) - lnGamma(a) - lnGamma(b) + u + Math.log(s);
1039      if( t < MINLOG )  s = 0.0;
1040      else                  s = Math.exp(t);
1041    }
1042    return s;
1043  }
1044 
1045  /**
1046   * Returns the revision string.
1047   *
1048   * @return            the revision
1049   */
1050  public String getRevision() {
1051    return RevisionUtils.extract("$Revision: 5616 $");
1052  }
1053
1054  /**
1055   * Main method for testing this class.
1056   */
1057  public static void main(String[] ops) {
1058
1059    System.out.println("Binomial standard error (0.5, 100): " + 
1060                       Statistics.binomialStandardError(0.5, 100));
1061    System.out.println("Chi-squared probability (2.558, 10): " +
1062                       Statistics.chiSquaredProbability(2.558, 10));
1063    System.out.println("Normal probability (0.2): " +
1064                       Statistics.normalProbability(0.2));
1065    System.out.println("F probability (5.1922, 4, 5): " +
1066                       Statistics.FProbability(5.1922, 4, 5));
1067    System.out.println("lnGamma(6): "+ Statistics.lnGamma(6));
1068  }
1069}
Note: See TracBrowser for help on using the repository browser.