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: [EMAIL PROTECTED] 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(Γ(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 [EMAIL PROTECTED] * * 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 [EMAIL PROTECTED] * * 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 [EMAIL PROTECTED] * * 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 * [EMAIL PROTECTED] 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 < 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 * [EMAIL PROTECTED] #inverseCummulativeProbability(double)} to find critical values. * * @param p the desired probability for the critical value * @return domain value lower bound, i.e. * P(X < <i>lower bound</i>) < <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 * [EMAIL PROTECTED] #inverseCummulativeProbability(double)} to find critical values. * * @param p the desired probability for the critical value * @return domain value upper bound, i.e. * P(X < <i>upper bound</i>) > <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 * [EMAIL PROTECTED] #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 [EMAIL PROTECTED] * * 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 * [EMAIL PROTECTED] 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 [EMAIL PROTECTED] * * 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: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]