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