commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From mdigg...@apache.org
Subject cvs commit: jakarta-commons-sandbox/math/src/java/org/apache/commons/math ContinuedFraction.java
Date Wed, 11 Jun 2003 01:19:19 GMT
mdiggory    2003/06/10 18:19:19

  Modified:    math/src/java/org/apache/commons/math/special Beta.java
                        Gamma.java
               math/src/test/org/apache/commons/math/stat/distribution
                        DistributionFactoryImplTest.java
               math/src/java/org/apache/commons/math/stat/distribution
                        DistributionFactory.java
                        DistributionFactoryImpl.java
  Added:       math/src/test/org/apache/commons/math/stat/distribution
                        FDistributionTest.java
               math/src/test/org/apache/commons/math
                        ContinuedFractionTest.java
               math/src/java/org/apache/commons/math/stat/distribution
                        FDistributionImpl.java FDistribution.java
               math/src/java/org/apache/commons/math ContinuedFraction.java
  Log:
  PR: http://nagoya.apache.org/bugzilla/show_bug.cgi?id=20601
  Submitted by:	brent@worden.org
  
  Revision  Changes    Path
  1.2       +62 -16    jakarta-commons-sandbox/math/src/java/org/apache/commons/math/special/Beta.java
  
  Index: Beta.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons-sandbox/math/src/java/org/apache/commons/math/special/Beta.java,v
  retrieving revision 1.1
  retrieving revision 1.2
  diff -u -r1.1 -r1.2
  --- Beta.java	7 Jun 2003 13:57:54 -0000	1.1
  +++ Beta.java	11 Jun 2003 01:19:18 -0000	1.2
  @@ -53,9 +53,11 @@
    */
   package org.apache.commons.math.special;
   
  +import org.apache.commons.math.ContinuedFraction;
  +
   /**
    * This is a utility class that provides computation methods related to the
  - * Gamma family of functions.
  + * Beta family of functions.
    * 
    * @author Brent Worden
    */
  @@ -124,8 +126,8 @@
        * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
        * Regularized Beta Function</a>.</li>
        * <li>
  -     * <a href="http://mathworld.wolfram.com/IncompleteBetaFunction.html">
  -     * Incomplete Beta Function</a>.</li>
  +     * <a href="http://functions.wolfram.com/06.21.10.0001.01">
  +     * Regularized Beta Function</a>.</li>
        * </ul>
        * </p>
        * 
  @@ -149,15 +151,45 @@
           } else if(x == 1.0){
               ret = 1.0;
           } else {
  -            double n = 0.0;
  -            double an = 1.0 / a;
  -            double s = an;
  -            while(Math.abs(an) > epsilon && n < maxIterations){
  -                n = n + 1.0;
  -                an = an * (n - b) / n * x / (a + n) * (a + n - 1);
  -                s = s + an;
  -            }
  -            ret = Math.exp(a * Math.log(x) - logBeta(a, b)) * s;
  +            ContinuedFraction fraction = new ContinuedFraction() {
  +                protected double getB(int n, double x) {
  +                    double ret;
  +                    double m;
  +                    switch (n) {
  +                        case 1 :
  +                            ret = 1.0;
  +                            break;
  +                        default :
  +                            if (n % 2 == 0) { // even
  +                                m = (n - 2.0) / 2.0;
  +                                ret =
  +                                    - ((a + m) * (a + b + m) * x)
  +                                        / ((a + (2 * m)) * (a + (2 * m) + 1.0));
  +                            } else {
  +                                m = (n - 1.0) / 2.0;
  +                                ret =
  +                                    (m * (b - m) * x)
  +                                        / ((a + (2 * m) - 1) * (a + (2 * m)));
  +                            }
  +                            break;
  +                    }
  +                    return ret;
  +                }
  +
  +                protected double getA(int n, double x) {
  +                    double ret;
  +                    switch (n) {
  +                        case 0 :
  +                            ret = 0.0;
  +                            break;
  +                        default :
  +                            ret = 1.0;
  +                            break;
  +                    }
  +                    return ret;
  +                }
  +			};
  +            ret = Math.exp((a * Math.log(x)) + (b * Math.log(1.0 - x)) - Math.log(a) -
logBeta(a, b, epsilon, maxIterations)) * fraction.evaluate(x, epsilon, maxIterations);
           }
   
           return ret;
  @@ -167,6 +199,19 @@
        * <p>
        * Returns the natural logarithm of the beta function B(a, b).
        * </p>
  +     * 
  +     * @param a ???
  +     * @param b ???
  +     * @return log(B(a, b))
  +     */
  +    public static double logBeta(double a, double b) {
  +        return logBeta(a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
  +    }
  +    
  +    /**
  +     * <p>
  +     * Returns the natural logarithm of the beta function B(a, b).
  +     * </p>
        *
        * <p> 
        * The implementation of this method is based on:
  @@ -176,10 +221,11 @@
        * </ul>
        * </p>
        * 
  -     * @param x ???
  +     * @param a ???
  +     * @param b ???
        * @return log(B(a, b))
        */
  -    public static double logBeta(double a, double b) {
  +    public static double logBeta(double a, double b, double epsilon, int maxIterations)
{
           double ret;
   
           if (a <= 0.0) {
  @@ -187,8 +233,8 @@
           } else if (b <= 0.0) {
               throw new IllegalArgumentException("b must be positive");
           } else {
  -            ret = Gamma.logGamma(a) + Gamma.logGamma(b)
  -                - Gamma.logGamma(a + b);
  +            ret = Gamma.logGamma(a, epsilon, maxIterations) + Gamma.logGamma(b, epsilon,
maxIterations)
  +                - Gamma.logGamma(a + b, epsilon, maxIterations);
           }
   
           return ret;
  
  
  
  1.4       +20 -12    jakarta-commons-sandbox/math/src/java/org/apache/commons/math/special/Gamma.java
  
  Index: Gamma.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons-sandbox/math/src/java/org/apache/commons/math/special/Gamma.java,v
  retrieving revision 1.3
  retrieving revision 1.4
  diff -u -r1.3 -r1.4
  --- Gamma.java	5 Jun 2003 14:03:53 -0000	1.3
  +++ Gamma.java	11 Jun 2003 01:19:18 -0000	1.4
  @@ -62,14 +62,9 @@
    * @author Brent Worden
    */
   public class Gamma {
  -    /** Maximum number of iteration allowed for iterative methods. */
  -    // TODO: try to reduce this.  regularizedGammaP doesn't converge very
  -    // fast for large values of x.
  -    private static final int MAXIMUM_ITERATIONS = 100;
  -
       /** Maximum allowed numerical error. */
  -    private static final double EPSILON = 10e-9;
  -
  +    private static final double DEFAULT_EPSILON = 10e-9;
  +    
       /**
        * Default constructor.  Prohibit instantiation.
        */
  @@ -82,6 +77,19 @@
        * Returns the regularized gamma function P(a, x).
        * </p>
        * 
  +     * @param a ???
  +     * @param x ???
  +     * @return the regularized gamma function P(a, x)
  +     */
  +    public static double regularizedGammaP(double a, double x) {
  +        return regularizedGammaP(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
  +    }
  +    
  +    /**
  +     * <p>
  +     * Returns the regularized gamma function P(a, x).
  +     * </p>
  +     * 
        * <p>
        * The implementation of this method is based on:
        * <ul>
  @@ -102,7 +110,7 @@
        * @param x ???
        * @return the regularized gamma function P(a, x)
        */
  -    public static double regularizedGammaP(double a, double x) {
  +    public static double regularizedGammaP(double a, double x, double epsilon, int maxIterations)
{
           double ret;
   
           if (a <= 0.0) {
  @@ -114,7 +122,7 @@
               double n = 0.0; // current element index
               double an = 1.0 / a; // n-th element in the series
               double sum = an; // partial sum
  -            while (Math.abs(an) > EPSILON && n < MAXIMUM_ITERATIONS) {
  +            while (Math.abs(an) > epsilon && n < maxIterations) {
                   // compute next element in the series
                   n = n + 1.0;
                   an = an * (x / (a + n));
  @@ -122,11 +130,11 @@
                   // update partial sum
                   sum = sum + an;
               }
  -            if (n >= MAXIMUM_ITERATIONS) {
  +            if (n >= maxIterations) {
                   throw new ConvergenceException(
                       "maximum number of iterations reached");
               } else {
  -                ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
  +                ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a, epsilon, maxIterations))
* sum;
               }
           }
   
  @@ -154,7 +162,7 @@
        * @param x ???
        * @return log(&#915;(x))
        */
  -    public static double logGamma(double x) {
  +    public static double logGamma(double x, double epsilon, int maxIterations) {
           double ret;
   
           if (x <= 0.0) {
  
  
  
  1.5       +44 -0     jakarta-commons-sandbox/math/src/test/org/apache/commons/math/stat/distribution/DistributionFactoryImplTest.java
  
  Index: DistributionFactoryImplTest.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons-sandbox/math/src/test/org/apache/commons/math/stat/distribution/DistributionFactoryImplTest.java,v
  retrieving revision 1.4
  retrieving revision 1.5
  diff -u -r1.4 -r1.5
  --- DistributionFactoryImplTest.java	7 Jun 2003 13:57:54 -0000	1.4
  +++ DistributionFactoryImplTest.java	11 Jun 2003 01:19:18 -0000	1.5
  @@ -58,6 +58,50 @@
           }
       }
       
  +    public void testCreateFDistributionNegativePositive(){
  +        try {
  +            factory.createFDistribution(-1.0, 1.0);
  +            fail("negative degrees of freedom.  IllegalArgumentException expected");
  +        } catch (IllegalArgumentException ex) {
  +            ;
  +        }
  +    }
  +    
  +    public void testCreateFDistributionZeroPositive(){
  +        try {
  +            factory.createFDistribution(0.0, 1.0);
  +            fail("zero degrees of freedom.  IllegalArgumentException expected");
  +        } catch (IllegalArgumentException ex) {
  +            ;
  +        }
  +    }
  +    
  +    public void testCreateFDistributionPositiveNegative(){
  +        try {
  +            factory.createFDistribution(1.0, -1.0);
  +            fail("negative degrees of freedom.  IllegalArgumentException expected");
  +        } catch (IllegalArgumentException ex) {
  +            ;
  +        }
  +    }
  +    
  +    public void testCreateFDistributionPositiveZero(){
  +        try {
  +            factory.createFDistribution(1.0, 0.0);
  +            fail("zero degrees of freedom.  IllegalArgumentException expected");
  +        } catch (IllegalArgumentException ex) {
  +            ;
  +        }
  +    }
  +    
  +    public void testCreateFDistributionPositivePositive(){
  +        try {
  +            factory.createFDistribution(1.0, 1.0);
  +        } catch (IllegalArgumentException ex) {
  +            fail("positive degrees of freedom.  IllegalArgumentException is not expected");
  +        }
  +    }
  +    
       public void testCreateGammaDistributionNegativePositive(){
           try {
               factory.createGammaDistribution(-1.0, 1.0);
  
  
  
  1.1                  jakarta-commons-sandbox/math/src/test/org/apache/commons/math/stat/distribution/FDistributionTest.java
  
  Index: FDistributionTest.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2003 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowlegement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowlegement may appear in the software itself,
   *    if and wherever such third-party acknowlegements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their names without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  package org.apache.commons.math.stat.distribution;
  
  import junit.framework.TestCase;
  
  /**
   * @author Brent Worden
   */
  public class FDistributionTest extends TestCase {
      private FDistribution f;
      
      /**
       * Constructor for ChiSquareDistributionTest.
       * @param name
       */
      public FDistributionTest(String name) {
          super(name);
      }
  
      /*
       * @see TestCase#setUp()
       */
      protected void setUp() throws Exception {
          super.setUp();
          f = DistributionFactory.newInstance().createFDistribution(5.0, 6.0);
      }
  
      /*
       * @see TestCase#tearDown()
       */
      protected void tearDown() throws Exception {
          f = null;
          super.tearDown();
      }
  
      public void testLowerTailProbability(){
          testProbability(1.0 / 10.67, .010);
          testProbability(1.0 /  6.98, .025);
          testProbability(1.0 /  4.95, .050);
          testProbability(1.0 /  3.40, .100);
      }
  
      public void testUpperTailProbability(){
          testProbability(8.75, .990);
          testProbability(5.99, .975);
          testProbability(4.39, .950);
          testProbability(3.11, .900);
      }
      
      public void testLowerTailValues(){
          testValue(1.0 / 10.67, .010);
          testValue(1.0 /  6.98, .025);
          testValue(1.0 /  4.95, .050);
          testValue(1.0 /  3.40, .100);
      }
      
      public void testUpperTailValues(){
          testValue(8.75, .990);
          testValue(5.99, .975);
          testValue(4.39, .950);
          testValue(3.11, .900);
      }
      
      private void testProbability(double x, double expected){
          double actual = f.cummulativeProbability(x);
          assertEquals("probability for " + x, expected, actual, 1e-3);
      }
      
      private void testValue(double expected, double p){
          double actual = f.inverseCummulativeProbability(p);
          assertEquals("value for " + p, expected, actual, 1e-2);
      }
  }
  
  
  
  1.1                  jakarta-commons-sandbox/math/src/test/org/apache/commons/math/ContinuedFractionTest.java
  
  Index: ContinuedFractionTest.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2003 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowlegement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowlegement may appear in the software itself,
   *    if and wherever such third-party acknowlegements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their names without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  package org.apache.commons.math;
  
  import junit.framework.TestCase;
  
  /**
   * @author Brent Worden
   */
  public class ContinuedFractionTest extends TestCase {
  	/**
  	 * Constructor for ContinuedFractionTest.
  	 * @param name
  	 */
  	public ContinuedFractionTest(String name) {
  		super(name);
  	}
  
  	public void testGoldenRation(){
          ContinuedFraction cf = new ContinuedFraction() {
  			public double getA(int n, double x) {
  				return 1.0;
  			}
  
  			public double getB(int n, double x) {
  				return 1.0;
  			}
  		};
          double gr = cf.evaluate(0.0, 10e-9);
          assertEquals(1.61803399, gr, 10e-9);
  	}
  }
  
  
  
  1.5       +12 -2     jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/DistributionFactory.java
  
  Index: DistributionFactory.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/DistributionFactory.java,v
  retrieving revision 1.4
  retrieving revision 1.5
  diff -u -r1.4 -r1.5
  --- DistributionFactory.java	7 Jun 2003 13:57:54 -0000	1.4
  +++ DistributionFactory.java	11 Jun 2003 01:19:19 -0000	1.5
  @@ -59,7 +59,9 @@
    * The following distributions are supported:
    * <ul>
    * <li>Chi-Squared</li>
  + * <li>F</li>
    * <li>Gamma</li>
  + * <li>Student's t</li>
    * </ul>
    * </p>
    * 
  @@ -99,8 +101,16 @@
        * @return a new chi-square distribution.  
        */
       public abstract ChiSquaredDistribution createChiSquareDistribution(
  -        double degreesOfFreedom
  -    );
  +        double degreesOfFreedom);
  +    
  +    /**
  +     * Create a new F-distribution with the given degrees of freedom.
  +     * @param numeratorDegreesOfFreedom numerator degrees of freedom.
  +     * @param denominatorDegreesOfFreedom denominator degrees of freedom.
  +     * @return a new F-distribution.  
  +     */
  +    public abstract FDistribution createFDistribution(
  +        double numeratorDegreesOfFreedom, double denominatorDegreesOfFreedom);
       
       /**
        * Create a new gamma distribution with the given alpha and beta values.
  
  
  
  1.5       +14 -0     jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/DistributionFactoryImpl.java
  
  Index: DistributionFactoryImpl.java
  ===================================================================
  RCS file: /home/cvs/jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/DistributionFactoryImpl.java,v
  retrieving revision 1.4
  retrieving revision 1.5
  diff -u -r1.4 -r1.5
  --- DistributionFactoryImpl.java	7 Jun 2003 13:57:54 -0000	1.4
  +++ DistributionFactoryImpl.java	11 Jun 2003 01:19:19 -0000	1.5
  @@ -99,4 +99,18 @@
       public TDistribution createTDistribution(double degreesOfFreedom) {
           return new TDistributionImpl(degreesOfFreedom);
       }
  +    
  +    /**
  +     * Create a new F-distribution with the given degrees of freedom.
  +     * @param numeratorDegreesOfFreedom numerator degrees of freedom.
  +     * @param denominatorDegreesOfFreedom denominator degrees of freedom.
  +     * @return a new F-distribution.  
  +     */
  +	public FDistribution createFDistribution(
  +		double numeratorDegreesOfFreedom,
  +		double denominatorDegreesOfFreedom) {
  +		return new FDistributionImpl(numeratorDegreesOfFreedom,
  +            denominatorDegreesOfFreedom);
  +	}
  +
   }
  
  
  
  1.1                  jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/FDistributionImpl.java
  
  Index: FDistributionImpl.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2003 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowlegement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowlegement may appear in the software itself,
   *    if and wherever such third-party acknowlegements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their names without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  package org.apache.commons.math.stat.distribution;
  
  import org.apache.commons.math.special.Beta;
  
  /**
   * Default implementation of
   * {@link org.apache.commons.math.stat.distribution.TDistribution}.
   * 
   * @author Brent Worden
   */
  public class FDistributionImpl
      extends AbstractContinuousDistribution
      implements FDistribution {
  
      /** The numerator degrees of freedom*/
      private double numeratorDegreesOfFreedom;
  
      /** The numerator degrees of freedom*/
      private double denominatorDegreesOfFreedom;
      
      /**
       * Create a F distribution using the given degrees of freedom.
       * @param degreesOfFreedom the degrees of freedom.
       */
      public FDistributionImpl(double numeratorDegreesOfFreedom,
              double denominatorDegreesOfFreedom){
          super();
          setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
          setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
      }
      
      /**
       * <p>
       * For this disbution, X, this method returns P(X &lt; x).
       * </p>
       * 
       * <p>
       * The implementation of this method is based on:
       * <ul>
       * <li>
       * <a href="http://mathworld.wolfram.com/F-Distribution.html">
       * F-Distribution</a>, equation (4).</li>
       * </p>
       * 
       * @param x the value at which the CDF is evaluated.
       * @return CDF for this distribution. 
       */
      public double cummulativeProbability(double x) {
          double ret;
          if(x <= 0.0){
              ret = 0.0;
          } else {
              double n = getNumeratorDegreesOfFreedom();
              double m = getDenominatorDegreesOfFreedom();
              
              ret = Beta.regularizedBeta((n * x) / (m + n * x),
                  0.5 * n,
                  0.5 * m);
          }
          return ret;
      }
          
      /**
       * Access the domain value lower bound, based on <code>p</code>, used to
       * bracket a CDF root.  This method is used by
       * {@link #inverseCummulativeProbability(double)} to find critical values.
       * 
       * @param p the desired probability for the critical value
       * @return domain value lower bound, i.e.
       *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>

       */
      protected double getDomainLowerBound(double p){
          return 0.0;
      }
  
      /**
       * Access the domain value upper bound, based on <code>p</code>, used to
       * bracket a CDF root.  This method is used by
       * {@link #inverseCummulativeProbability(double)} to find critical values.
       * 
       * @param p the desired probability for the critical value
       * @return domain value upper bound, i.e.
       *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>

       */
      protected double getDomainUpperBound(double p){
          return Double.MAX_VALUE;
      }
  
      /**
       * Access the initial domain value, based on <code>p</code>, used to
       * bracket a CDF root.  This method is used by
       * {@link #inverseCummulativeProbability(double)} to find critical values.
       * 
       * @param p the desired probability for the critical value
       * @return initial domain value
       */
      protected double getInitialDomain(double p){
          return getDenominatorDegreesOfFreedom() / (getDenominatorDegreesOfFreedom() - 2.0);
      }
      
      /**
       * Modify the numerator degrees of freedom.
       * @param degreesOfFreedom the new numerator degrees of freedom.
       */
      public void setNumeratorDegreesOfFreedom(double degreesOfFreedom){
          if(degreesOfFreedom <= 0.0){
              throw new IllegalArgumentException(
                  "degrees of freedom must be positive.");
          }
          this.numeratorDegreesOfFreedom = degreesOfFreedom;
      }
      
      /**
       * Access the numerator degrees of freedom.
       * @return the numerator degrees of freedom.
       */
      public double getNumeratorDegreesOfFreedom(){
          return numeratorDegreesOfFreedom;
      }
      
      /**
       * Modify the denominator degrees of freedom.
       * @param degreesOfFreedom the new denominator degrees of freedom.
       */
      public void setDenominatorDegreesOfFreedom(double degreesOfFreedom){
          if(degreesOfFreedom <= 0.0){
              throw new IllegalArgumentException(
                  "degrees of freedom must be positive.");
          }
          this.denominatorDegreesOfFreedom = degreesOfFreedom;
      }
      
      /**
       * Access the denominator degrees of freedom.
       * @return the denominator degrees of freedom.
       */
      public double getDenominatorDegreesOfFreedom(){
          return denominatorDegreesOfFreedom;
      }
  }
  
  
  
  1.1                  jakarta-commons-sandbox/math/src/java/org/apache/commons/math/stat/distribution/FDistribution.java
  
  Index: FDistribution.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2003 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowlegement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowlegement may appear in the software itself,
   *    if and wherever such third-party acknowlegements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their names without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  package org.apache.commons.math.stat.distribution;
  
  /**
   * <p>
   * F-Distribution.
   * </p>
   * 
   * <p>
   * Instances of FDistribution objects should be created using
   * {@link DistributionFactory#createFDistribution(double)}
   * </p>
   * 
   * <p>
   * Reference:<br/>
   * <a href="http://mathworld.wolfram.com/F-Distribution.html">
   * F-Distribution</a>
   * </p>
   * 
   * @author Brent Worden
   */
  public interface FDistribution extends ContinuousDistribution {
      /**
       * Modify the numerator degrees of freedom.
       * @param degreesOfFreedom the new numerator degrees of freedom.
       */
      void setNumeratorDegreesOfFreedom(double degreesOfFreedom);
      
      /**
       * Access the numerator degrees of freedom.
       * @return the numerator degrees of freedom.
       */
      double getNumeratorDegreesOfFreedom();
      
      /**
       * Modify the denominator degrees of freedom.
       * @param degreesOfFreedom the new denominator degrees of freedom.
       */
      void setDenominatorDegreesOfFreedom(double degreesOfFreedom);
      
      /**
       * Access the denominator degrees of freedom.
       * @return the denominator degrees of freedom.
       */
      double getDenominatorDegreesOfFreedom();
  }
  
  
  
  1.1                  jakarta-commons-sandbox/math/src/java/org/apache/commons/math/ContinuedFraction.java
  
  Index: ContinuedFraction.java
  ===================================================================
  /* ====================================================================
   * The Apache Software License, Version 1.1
   *
   * Copyright (c) 2003 The Apache Software Foundation.  All rights
   * reserved.
   *
   * Redistribution and use in source and binary forms, with or without
   * modification, are permitted provided that the following conditions
   * are met:
   *
   * 1. Redistributions of source code must retain the above copyright
   *    notice, this list of conditions and the following disclaimer.
   *
   * 2. Redistributions in binary form must reproduce the above copyright
   *    notice, this list of conditions and the following disclaimer in
   *    the documentation and/or other materials provided with the
   *    distribution.
   *
   * 3. The end-user documentation included with the redistribution, if
   *    any, must include the following acknowlegement:
   *       "This product includes software developed by the
   *        Apache Software Foundation (http://www.apache.org/)."
   *    Alternately, this acknowlegement may appear in the software itself,
   *    if and wherever such third-party acknowlegements normally appear.
   *
   * 4. The names "The Jakarta Project", "Commons", and "Apache Software
   *    Foundation" must not be used to endorse or promote products derived
   *    from this software without prior written permission. For written
   *    permission, please contact apache@apache.org.
   *
   * 5. Products derived from this software may not be called "Apache"
   *    nor may "Apache" appear in their names without prior written
   *    permission of the Apache Software Foundation.
   *
   * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
   * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
   * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   * DISCLAIMED.  IN NO EVENT SHALL THE APACHE SOFTWARE FOUNDATION OR
   * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
   * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
   * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
   * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
   * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
   * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
   * SUCH DAMAGE.
   * ====================================================================
   *
   * This software consists of voluntary contributions made by many
   * individuals on behalf of the Apache Software Foundation.  For more
   * information on the Apache Software Foundation, please see
   * <http://www.apache.org/>.
   */
  package org.apache.commons.math;
  
  /**
   * <p>
   * Provides a generic means to evaluate continued fractions.  Subclasses simply
   * provided the a and b coefficients to evaluate the continued fraction.
   * </p>
   * 
   * <p>
   * Reference:<br/>
   * <a href="http://mathworld.wolfram.com/ContinuedFraction.html">
   * Continued Fraction</a>
   * </p>
   * 
   * @author Brent Worden
   */
  public abstract class ContinuedFraction {
      /** Maximum allowed numerical error. */
      private static final double DEFAULT_EPSILON = 10e-9;
      
  	/**
  	 * Default constructor.
  	 */
  	protected ContinuedFraction() {
  		super();
  	}
      
      /**
       * Access the n-th a coefficient of the continued fraction.  Since a can be
       * a function of the evaluation point, x, that is passed in as well.
       * @param n the coefficient index to retrieve.
       * @param x the evaluation point.
       * @return the n-th a coefficient.
       */
      protected abstract double getA(int n, double x);
      
      /**
       * Access the n-th b coefficient of the continued fraction.  Since b can be
       * a function of the evaluation point, x, that is passed in as well.
       * @param n the coefficient index to retrieve.
       * @param x the evaluation point.
       * @return the n-th b coefficient.
       */
      protected abstract double getB(int n, double x);
      
      /**
       * Evaluates the continued fraction at the value x.
       * @param x the evaluation point.
       * @return the value of the continued fraction evaluated at x. 
       */
      public double evaluate(double x){
          return evaluate(x, DEFAULT_EPSILON, Integer.MAX_VALUE);
      }    
      
      /**
       * Evaluates the continued fraction at the value x.
       * @param x the evaluation point.
       * @param epsilon maximum error allowed.
       * @return the value of the continued fraction evaluated at x. 
       */
      public double evaluate(double x, double epsilon){
          return evaluate(x, epsilon, Integer.MAX_VALUE);
      }    
      
      /**
       * Evaluates the continued fraction at the value x.
       * @param x the evaluation point.
       * @param maxIterations maximum number of convergents
       * @return the value of the continued fraction evaluated at x. 
       */
      public double evaluate(double x, int maxIterations){
          return evaluate(x, DEFAULT_EPSILON, maxIterations);
      }    
      
      /**
       * <p>
       * Evaluates the continued fraction at the value x.
       * </p>
       * 
       * <p>
       * The implementation of this method is based on:
       * <ul>
       * <li>O. E-gecio-glu, C . K. Koc, J. Rifa i Coma,
       * <a href="http://citeseer.nj.nec.com/egecioglu91fast.html">
       * Fast Computation of Continued Fractions</a>, Computers Math. Applic.,
       * 21(2--3), 1991, 167--169.</li>
       * </ul>
       * </p>
       * 
       * @param x the evaluation point.
       * @param epsilon maximum error allowed.
       * @param maxIterations maximum number of convergents
       * @return the value of the continued fraction evaluated at x. 
       */
      public double evaluate(double x, double epsilon, int maxIterations) {
          double[][] f = new double[2][2];
          double[][] a = new double[2][2];
          double[][] an = new double[2][2];
          
          a[0][0] = getA(0, x);
          a[0][1] = 1.0;
          a[1][0] = 1.0;
          a[1][1] = 0.0;
  
          return evaluate(1, x, a, an, f, epsilon, maxIterations);
      }
      
      /**
       * Evaluates the n-th convergent, fn = pn / qn, for this continued fraction at the value
x.
       * @param n the convergent to compute.
       * @param x the evaluation point.
       * @param a (n-1)-th convergent matrix.  (Input)
       * @param an the n-th coefficient matrix. (Output)
       * @param f the n-th convergent matrix. (Output)
       * @param epsilon maximum error allowed.
       * @param maxIterations maximum number of convergents
       * @return the value of the the n-th convergent for this continued fraction evaluated
at x. 
       */
      private double evaluate(int n, double x, double[][] a, double[][] an, double[][] f,
double epsilon, int maxIterations) {
          double ret;
          
          // create next matrix
          an[0][0] = getA(n, x);
          an[0][1] = 1.0;
          an[1][0] = getB(n, x);
          an[1][1] = 0.0;
          
          // multiply a and an, save as f
          f[0][0] = (a[0][0] * an[0][0]) + (a[0][1] * an[1][0]);
          f[0][1] = (a[0][0] * an[0][1]) + (a[0][1] * an[1][1]);
          f[1][0] = (a[1][0] * an[0][0]) + (a[1][1] * an[1][0]);
          f[1][1] = (a[1][0] * an[0][1]) + (a[1][1] * an[1][1]);
          
          // determine if we're close enough
          if(Math.abs((f[0][0] * f[1][1]) - (f[1][0] * f[0][1])) < Math.abs(epsilon * f[1][0]
* f[1][1])){
              ret = f[0][0] / f[1][0];
          } else {
              if(n >= maxIterations){
                  throw new ConvergenceException("Continued fraction convergents failed to
converge.");
              }
              // compute next
              ret = evaluate(n + 1, x, f /* new a */, an /* reuse an */, a /* new f */, epsilon,
maxIterations);
          }
          
          return ret;
      }
  }
  
  
  

---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message