source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/fdistr.cs @ 4068

/*************************************************************************
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier

Contributors:
    * Sergey Bochkanov (ALGLIB project). Translation from C to
      pseudocode.

See subroutines comments for additional copyrights.

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/


namespace alglib {
  public class fdistr {
    /*************************************************************************
    F distribution

    Returns the area from zero to x under the F density
    function (also known as Snedcor's density or the
    variance ratio density).  This is the density
    of x = (u1/df1)/(u2/df2), where u1 and u2 are random
    variables having Chi square distributions with df1
    and df2 degrees of freedom, respectively.
    The incomplete beta integral is used, according to the
    formula

    P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ).


    The arguments a and b are greater than zero, and x is
    nonnegative.

    ACCURACY:

    Tested at random points (a,b,x).

                   x     a,b                     Relative error:
    arithmetic  domain  domain     # trials      peak         rms
       IEEE      0,1    0,100       100000      9.8e-15     1.7e-15
       IEEE      1,5    0,100       100000      6.5e-15     3.5e-16
       IEEE      0,1    1,10000     100000      2.2e-11     3.3e-12
       IEEE      1,5    1,10000     100000      1.1e-11     1.7e-13

    Cephes Math Library Release 2.8:  June, 2000
    Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    public static double fdistribution(int a,
        int b,
        double x) {
      double result = 0;
      double w = 0;

      System.Diagnostics.Debug.Assert(a >= 1 & b >= 1 & (double)(x) >= (double)(0), "Domain error in FDistribution");
      w = a * x;
      w = w / (b + w);
      result = ibetaf.incompletebeta(0.5 * a, 0.5 * b, w);
      return result;
    }


    /*************************************************************************
    Complemented F distribution

    Returns the area from x to infinity under the F density
    function (also known as Snedcor's density or the
    variance ratio density).


                         inf.
                          -
                 1       | |  a-1      b-1
    1-P(x)  =  ------    |   t    (1-t)    dt
               B(a,b)  | |
                        -
                         x


    The incomplete beta integral is used, according to the
    formula

    P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ).


    ACCURACY:

    Tested at random points (a,b,x) in the indicated intervals.
                   x     a,b                     Relative error:
    arithmetic  domain  domain     # trials      peak         rms
       IEEE      0,1    1,100       100000      3.7e-14     5.9e-16
       IEEE      1,5    1,100       100000      8.0e-15     1.6e-15
       IEEE      0,1    1,10000     100000      1.8e-11     3.5e-13
       IEEE      1,5    1,10000     100000      2.0e-11     3.0e-12

    Cephes Math Library Release 2.8:  June, 2000
    Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    public static double fcdistribution(int a,
        int b,
        double x) {
      double result = 0;
      double w = 0;

      System.Diagnostics.Debug.Assert(a >= 1 & b >= 1 & (double)(x) >= (double)(0), "Domain error in FCDistribution");
      w = b / (b + a * x);
      result = ibetaf.incompletebeta(0.5 * b, 0.5 * a, w);
      return result;
    }


    /*************************************************************************
    Inverse of complemented F distribution

    Finds the F density argument x such that the integral
    from x to infinity of the F density is equal to the
    given probability p.

    This is accomplished using the inverse beta integral
    function and the relations

         z = incbi( df2/2, df1/2, p )
         x = df2 (1-z) / (df1 z).

    Note: the following relations hold for the inverse of
    the uncomplemented F distribution:

         z = incbi( df1/2, df2/2, p )
         x = df2 z / (df1 (1-z)).

    ACCURACY:

    Tested at random points (a,b,p).

                 a,b                     Relative error:
    arithmetic  domain     # trials      peak         rms
     For p between .001 and 1:
       IEEE     1,100       100000      8.3e-15     4.7e-16
       IEEE     1,10000     100000      2.1e-11     1.4e-13
     For p between 10^-6 and 10^-3:
       IEEE     1,100        50000      1.3e-12     8.4e-15
       IEEE     1,10000      50000      3.0e-12     4.8e-14

    Cephes Math Library Release 2.8:  June, 2000
    Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
    *************************************************************************/
    public static double invfdistribution(int a,
        int b,
        double y) {
      double result = 0;
      double w = 0;

      System.Diagnostics.Debug.Assert(a >= 1 & b >= 1 & (double)(y) > (double)(0) & (double)(y) <= (double)(1), "Domain error in InvFDistribution");

      //
      // Compute probability for x = 0.5
      //
      w = ibetaf.incompletebeta(0.5 * b, 0.5 * a, 0.5);

      //
      // If that is greater than y, then the solution w < .5
      // Otherwise, solve at 1-y to remove cancellation in (b - b*w)
      //
      if ((double)(w) > (double)(y) | (double)(y) < (double)(0.001)) {
        w = ibetaf.invincompletebeta(0.5 * b, 0.5 * a, y);
        result = (b - b * w) / (a * w);
      } else {
        w = ibetaf.invincompletebeta(0.5 * a, 0.5 * b, 1.0 - y);
        result = b * w / (a * (1.0 - w));
      }
      return result;
    }
  }
}