From b58585fb8def076b355699c8ee496b889534d744 Mon Sep 17 00:00:00 2001 From: "Mark R. Diggory" Date: Wed, 11 Jun 2003 01:19:19 +0000 Subject: [PATCH] PR: http://nagoya.apache.org/bugzilla/show_bug.cgi?id=20601 Submitted by: brent@worden.org git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@140899 13f79535-47bb-0310-9956-ffa450edef68 --- .../commons/math/ContinuedFraction.java | 200 ++++++++++++++++++ .../distribution/DistributionFactory.java | 14 +- .../distribution/DistributionFactoryImpl.java | 14 ++ .../math/distribution/FDistribution.java | 98 +++++++++ .../math/distribution/FDistributionImpl.java | 193 +++++++++++++++++ .../org/apache/commons/math/special/Beta.java | 78 +++++-- .../apache/commons/math/special/Gamma.java | 32 +-- .../commons/math/ContinuedFractionTest.java | 83 ++++++++ .../DistributionFactoryImplTest.java | 44 ++++ .../math/distribution/FDistributionTest.java | 125 +++++++++++ 10 files changed, 851 insertions(+), 30 deletions(-) create mode 100644 src/java/org/apache/commons/math/ContinuedFraction.java create mode 100644 src/java/org/apache/commons/math/distribution/FDistribution.java create mode 100644 src/java/org/apache/commons/math/distribution/FDistributionImpl.java create mode 100644 src/test/org/apache/commons/math/ContinuedFractionTest.java create mode 100644 src/test/org/apache/commons/math/distribution/FDistributionTest.java diff --git a/src/java/org/apache/commons/math/ContinuedFraction.java b/src/java/org/apache/commons/math/ContinuedFraction.java new file mode 100644 index 000000000..b2c54c420 --- /dev/null +++ b/src/java/org/apache/commons/math/ContinuedFraction.java @@ -0,0 +1,200 @@ +/* ==================================================================== + * 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 + * . + */ +package org.apache.commons.math; + +/** + *

+ * Provides a generic means to evaluate continued fractions. Subclasses simply + * provided the a and b coefficients to evaluate the continued fraction. + *

+ * + *

+ * Reference:
+ * + * Continued Fraction + *

+ * + * @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); + } + + /** + *

+ * Evaluates the continued fraction at the value x. + *

+ * + *

+ * The implementation of this method is based on: + *

+ *

+ * + * @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; + } +} diff --git a/src/java/org/apache/commons/math/distribution/DistributionFactory.java b/src/java/org/apache/commons/math/distribution/DistributionFactory.java index 0bd38cd22..048bb0d35 100644 --- a/src/java/org/apache/commons/math/distribution/DistributionFactory.java +++ b/src/java/org/apache/commons/math/distribution/DistributionFactory.java @@ -59,7 +59,9 @@ package org.apache.commons.math.stat.distribution; * The following distributions are supported: * *

* @@ -99,8 +101,16 @@ public abstract class DistributionFactory { * @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. diff --git a/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java b/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java index 648b6d06e..b172dca77 100644 --- a/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java +++ b/src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java @@ -99,4 +99,18 @@ public class DistributionFactoryImpl extends DistributionFactory { 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); + } + } diff --git a/src/java/org/apache/commons/math/distribution/FDistribution.java b/src/java/org/apache/commons/math/distribution/FDistribution.java new file mode 100644 index 000000000..610c5ff77 --- /dev/null +++ b/src/java/org/apache/commons/math/distribution/FDistribution.java @@ -0,0 +1,98 @@ +/* ==================================================================== + * 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 + * . + */ +package org.apache.commons.math.stat.distribution; + +/** + *

+ * F-Distribution. + *

+ * + *

+ * Instances of FDistribution objects should be created using + * {@link DistributionFactory#createFDistribution(double)} + *

+ * + *

+ * Reference:
+ * + * F-Distribution + *

+ * + * @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(); +} diff --git a/src/java/org/apache/commons/math/distribution/FDistributionImpl.java b/src/java/org/apache/commons/math/distribution/FDistributionImpl.java new file mode 100644 index 000000000..76db34fc8 --- /dev/null +++ b/src/java/org/apache/commons/math/distribution/FDistributionImpl.java @@ -0,0 +1,193 @@ +/* ==================================================================== + * 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 + * . + */ +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); + } + + /** + *

+ * For this disbution, X, this method returns P(X < x). + *

+ * + *

+ * The implementation of this method is based on: + *

    + *
  • + * + * F-Distribution, equation (4).
    • + *

    + * + * @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 p, 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 < lower bound) < p + */ + protected double getDomainLowerBound(double p){ + return 0.0; + } + + /** + * Access the domain value upper bound, based on p, 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 < upper bound) > p + */ + protected double getDomainUpperBound(double p){ + return Double.MAX_VALUE; + } + + /** + * Access the initial domain value, based on p, 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; + } +} diff --git a/src/java/org/apache/commons/math/special/Beta.java b/src/java/org/apache/commons/math/special/Beta.java index a6f25fc22..39f91848e 100644 --- a/src/java/org/apache/commons/math/special/Beta.java +++ b/src/java/org/apache/commons/math/special/Beta.java @@ -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 @@ public class Beta { * * Regularized Beta Function. *
  • - * - * Incomplete Beta Function.
    • + * + * Regularized Beta Function. *
*

* @@ -149,20 +151,63 @@ public class Beta { } 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; } + /** + *

+ * Returns the natural logarithm of the beta function B(a, b). + *

+ * + * @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); + } + /** *

* Returns the natural logarithm of the beta function B(a, b). @@ -176,10 +221,11 @@ public class Beta { * *

* - * @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 @@ public class Beta { } 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; diff --git a/src/java/org/apache/commons/math/special/Gamma.java b/src/java/org/apache/commons/math/special/Gamma.java index 84821c55c..dbd98f388 100644 --- a/src/java/org/apache/commons/math/special/Gamma.java +++ b/src/java/org/apache/commons/math/special/Gamma.java @@ -62,14 +62,9 @@ import org.apache.commons.math.ConvergenceException; * @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. */ @@ -77,6 +72,19 @@ public class Gamma { super(); } + /** + *

+ * Returns the regularized gamma function P(a, x). + *

+ * + * @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); + } + /** *

* Returns the regularized gamma function P(a, x). @@ -102,7 +110,7 @@ public class Gamma { * @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 @@ public class Gamma { 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 @@ public class Gamma { // 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 @@ public class Gamma { * @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) { diff --git a/src/test/org/apache/commons/math/ContinuedFractionTest.java b/src/test/org/apache/commons/math/ContinuedFractionTest.java new file mode 100644 index 000000000..fa55e26f1 --- /dev/null +++ b/src/test/org/apache/commons/math/ContinuedFractionTest.java @@ -0,0 +1,83 @@ +/* ==================================================================== + * 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 + * . + */ +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); + } +} diff --git a/src/test/org/apache/commons/math/distribution/DistributionFactoryImplTest.java b/src/test/org/apache/commons/math/distribution/DistributionFactoryImplTest.java index 4ee86fbd4..e66659663 100644 --- a/src/test/org/apache/commons/math/distribution/DistributionFactoryImplTest.java +++ b/src/test/org/apache/commons/math/distribution/DistributionFactoryImplTest.java @@ -58,6 +58,50 @@ public class DistributionFactoryImplTest extends TestCase { } } + 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); diff --git a/src/test/org/apache/commons/math/distribution/FDistributionTest.java b/src/test/org/apache/commons/math/distribution/FDistributionTest.java new file mode 100644 index 000000000..bad6b2d22 --- /dev/null +++ b/src/test/org/apache/commons/math/distribution/FDistributionTest.java @@ -0,0 +1,125 @@ +/* ==================================================================== + * 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 + * . + */ +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); + } +}