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:
*
* - Chi-Squared
+ * - F
* - Gamma
+ * - Student's t
*
*
*
@@ -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);
+ }
+}