Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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
This commit is contained in:
Mark R. Diggory 2003-06-11 01:19:19 +00:00
parent 8d1482c461
commit b58585fb8d
10 changed files with 851 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

200
src/java/org/apache/commons/math/ContinuedFraction.java Normal file
View File

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

14
src/java/org/apache/commons/math/distribution/DistributionFactory.java
View File

@ -59,7 +59,9 @@ package org.apache.commons.math.stat.distribution;
* The following distributions are supported: * The following distributions are supported:
* <ul> * <ul>
* <li>Chi-Squared</li> * <li>Chi-Squared</li>
* <li>F</li>
* <li>Gamma</li> * <li>Gamma</li>
* <li>Student's t</li>
* </ul> * </ul>
* </p> * </p>
* *
@ -99,8 +101,16 @@ public abstract class DistributionFactory {
* @return a new chi-square distribution. * @return a new chi-square distribution.
*/ */
public abstract ChiSquaredDistribution createChiSquareDistribution( 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. * Create a new gamma distribution with the given alpha and beta values.

14
src/java/org/apache/commons/math/distribution/DistributionFactoryImpl.java
View File

@ -99,4 +99,18 @@ public class DistributionFactoryImpl extends DistributionFactory {
public TDistribution createTDistribution(double degreesOfFreedom) { public TDistribution createTDistribution(double degreesOfFreedom) {
return new TDistributionImpl(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);
}
} }

98
src/java/org/apache/commons/math/distribution/FDistribution.java Normal file
View File

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

193
src/java/org/apache/commons/math/distribution/FDistributionImpl.java Normal file
View File

@ -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
* <http://www.apache.org/>.
*/
package org.apache.commons.math.stat.distribution;
import org.apache.commons.math.special.Beta;
/**
* Default implementation of
* {@link org.apache.commons.math.stat.distribution.TDistribution}.
*
* @author Brent Worden
*/
public class FDistributionImpl
extends AbstractContinuousDistribution
implements FDistribution {
/** The numerator degrees of freedom*/
private double numeratorDegreesOfFreedom;
/** The numerator degrees of freedom*/
private double denominatorDegreesOfFreedom;
/**
* Create a F distribution using the given degrees of freedom.
* @param degreesOfFreedom the degrees of freedom.
*/
public FDistributionImpl(double numeratorDegreesOfFreedom,
double denominatorDegreesOfFreedom){
super();
setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
}
/**
* <p>
* For this disbution, X, this method returns P(X &lt; x).
* </p>
*
* <p>
* The implementation of this method is based on:
* <ul>
* <li>
* <a href="http://mathworld.wolfram.com/F-Distribution.html">
* F-Distribution</a>, equation (4).</li>
* </p>
*
* @param x the value at which the CDF is evaluated.
* @return CDF for this distribution.
*/
public double cummulativeProbability(double x) {
double ret;
if(x <= 0.0){
ret = 0.0;
} else {
double n = getNumeratorDegreesOfFreedom();
double m = getDenominatorDegreesOfFreedom();
ret = Beta.regularizedBeta((n * x) / (m + n * x),
0.5 * n,
0.5 * m);
}
return ret;
}
/**
* Access the domain value lower bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* {@link #inverseCummulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
* @return domain value lower bound, i.e.
* P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
*/
protected double getDomainLowerBound(double p){
return 0.0;
}
/**
* Access the domain value upper bound, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* {@link #inverseCummulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
* @return domain value upper bound, i.e.
* P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
*/
protected double getDomainUpperBound(double p){
return Double.MAX_VALUE;
}
/**
* Access the initial domain value, based on <code>p</code>, used to
* bracket a CDF root. This method is used by
* {@link #inverseCummulativeProbability(double)} to find critical values.
*
* @param p the desired probability for the critical value
* @return initial domain value
*/
protected double getInitialDomain(double p){
return getDenominatorDegreesOfFreedom() / (getDenominatorDegreesOfFreedom() - 2.0);
}
/**
* Modify the numerator degrees of freedom.
* @param degreesOfFreedom the new numerator degrees of freedom.
*/
public void setNumeratorDegreesOfFreedom(double degreesOfFreedom){
if(degreesOfFreedom <= 0.0){
throw new IllegalArgumentException(
"degrees of freedom must be positive.");
}
this.numeratorDegreesOfFreedom = degreesOfFreedom;
}
/**
* Access the numerator degrees of freedom.
* @return the numerator degrees of freedom.
*/
public double getNumeratorDegreesOfFreedom(){
return numeratorDegreesOfFreedom;
}
/**
* Modify the denominator degrees of freedom.
* @param degreesOfFreedom the new denominator degrees of freedom.
*/
public void setDenominatorDegreesOfFreedom(double degreesOfFreedom){
if(degreesOfFreedom <= 0.0){
throw new IllegalArgumentException(
"degrees of freedom must be positive.");
}
this.denominatorDegreesOfFreedom = degreesOfFreedom;
}
/**
* Access the denominator degrees of freedom.
* @return the denominator degrees of freedom.
*/
public double getDenominatorDegreesOfFreedom(){
return denominatorDegreesOfFreedom;
}
}

76
src/java/org/apache/commons/math/special/Beta.java
View File

@ -53,9 +53,11 @@
*/ */
package org.apache.commons.math.special; package org.apache.commons.math.special;
import org.apache.commons.math.ContinuedFraction;
/** /**
* This is a utility class that provides computation methods related to the * This is a utility class that provides computation methods related to the
* Gamma family of functions. * Beta family of functions.
* *
* @author Brent Worden * @author Brent Worden
*/ */
@ -124,8 +126,8 @@ public class Beta {
* <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html"> * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
* Regularized Beta Function</a>.</li> * Regularized Beta Function</a>.</li>
* <li> * <li>
* <a href="http://mathworld.wolfram.com/IncompleteBetaFunction.html"> * <a href="http://functions.wolfram.com/06.21.10.0001.01">
* Incomplete Beta Function</a>.</li> * Regularized Beta Function</a>.</li>
* </ul> * </ul>
* </p> * </p>
* *
@ -149,20 +151,63 @@ public class Beta {
} else if(x == 1.0){ } else if(x == 1.0){
ret = 1.0; ret = 1.0;
} else { } else {
double n = 0.0; ContinuedFraction fraction = new ContinuedFraction() {
double an = 1.0 / a; protected double getB(int n, double x) {
double s = an; double ret;
while(Math.abs(an) > epsilon && n < maxIterations){ double m;
n = n + 1.0; switch (n) {
an = an * (n - b) / n * x / (a + n) * (a + n - 1); case 1 :
s = s + an; 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)));
} }
ret = Math.exp(a * Math.log(x) - logBeta(a, b)) * s; 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; return ret;
} }
/**
* <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> * <p>
* Returns the natural logarithm of the beta function B(a, b). * Returns the natural logarithm of the beta function B(a, b).
@ -176,10 +221,11 @@ public class Beta {
* </ul> * </ul>
* </p> * </p>
* *
* @param x ??? * @param a ???
* @param b ???
* @return log(B(a, 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; double ret;
if (a <= 0.0) { if (a <= 0.0) {
@ -187,8 +233,8 @@ public class Beta {
} else if (b <= 0.0) { } else if (b <= 0.0) {
throw new IllegalArgumentException("b must be positive"); throw new IllegalArgumentException("b must be positive");
} else { } else {
ret = Gamma.logGamma(a) + Gamma.logGamma(b) ret = Gamma.logGamma(a, epsilon, maxIterations) + Gamma.logGamma(b, epsilon, maxIterations)
- Gamma.logGamma(a + b); - Gamma.logGamma(a + b, epsilon, maxIterations);
} }
return ret; return ret;

30
src/java/org/apache/commons/math/special/Gamma.java
View File

@ -62,13 +62,8 @@ import org.apache.commons.math.ConvergenceException;
* @author Brent Worden * @author Brent Worden
*/ */
public class Gamma { 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. */ /** Maximum allowed numerical error. */
private static final double EPSILON = 10e-9; private static final double DEFAULT_EPSILON = 10e-9;
/** /**
* Default constructor. Prohibit instantiation. * Default constructor. Prohibit instantiation.
@ -77,6 +72,19 @@ public class Gamma {
super(); super();
} }
/**
* <p>
* 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> * <p>
* Returns the regularized gamma function P(a, x). * Returns the regularized gamma function P(a, x).
@ -102,7 +110,7 @@ public class Gamma {
* @param x ??? * @param x ???
* @return the regularized gamma function P(a, 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; double ret;
if (a <= 0.0) { if (a <= 0.0) {
@ -114,7 +122,7 @@ public class Gamma {
double n = 0.0; // current element index double n = 0.0; // current element index
double an = 1.0 / a; // n-th element in the series double an = 1.0 / a; // n-th element in the series
double sum = an; // partial sum 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 // compute next element in the series
n = n + 1.0; n = n + 1.0;
an = an * (x / (a + n)); an = an * (x / (a + n));
@ -122,11 +130,11 @@ public class Gamma {
// update partial sum // update partial sum
sum = sum + an; sum = sum + an;
} }
if (n >= MAXIMUM_ITERATIONS) { if (n >= maxIterations) {
throw new ConvergenceException( throw new ConvergenceException(
"maximum number of iterations reached"); "maximum number of iterations reached");
} else { } 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 ??? * @param x ???
* @return log(&#915;(x)) * @return log(&#915;(x))
*/ */
public static double logGamma(double x) { public static double logGamma(double x, double epsilon, int maxIterations) {
double ret; double ret;
if (x <= 0.0) { if (x <= 0.0) {

83
src/test/org/apache/commons/math/ContinuedFractionTest.java Normal file
View File

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

44
src/test/org/apache/commons/math/distribution/DistributionFactoryImplTest.java
View File

@ -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(){ public void testCreateGammaDistributionNegativePositive(){
try { try {
factory.createGammaDistribution(-1.0, 1.0); factory.createGammaDistribution(-1.0, 1.0);

125
src/test/org/apache/commons/math/distribution/FDistributionTest.java Normal file
View File

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