MATH-227. fixed F distribution inverse CDF computation for small denominator degrees of freedom.
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@699157 13f79535-47bb-0310-9956-ffa450edef68
diff --git a/src/java/org/apache/commons/math/distribution/FDistributionImpl.java b/src/java/org/apache/commons/math/distribution/FDistributionImpl.java
index 3959403..4e85614 100644
--- a/src/java/org/apache/commons/math/distribution/FDistributionImpl.java
+++ b/src/java/org/apache/commons/math/distribution/FDistributionImpl.java
@@ -24,168 +24,188 @@
/**
* Default implementation of
* {@link org.apache.commons.math.distribution.FDistribution}.
- *
- * @version $Revision$ $Date$
+ *
+ * @version $Revision$ $Date: 2008-02-08 09:44:11 -0600 (Fri, 08 Feb
+ * 2008) $
*/
-public class FDistributionImpl
- extends AbstractContinuousDistribution
- implements FDistribution, Serializable {
+public class FDistributionImpl extends AbstractContinuousDistribution implements
+ FDistribution, Serializable {
- /** Serializable version identifier */
- private static final long serialVersionUID = -8516354193418641566L;
+ /** Serializable version identifier */
+ private static final long serialVersionUID = -8516354193418641566L;
- /** The numerator degrees of freedom*/
- private double numeratorDegreesOfFreedom;
+ /** 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 numeratorDegreesOfFreedom the numerator degrees of freedom.
- * @param denominatorDegreesOfFreedom the denominator degrees of freedom.
- */
- public FDistributionImpl(double numeratorDegreesOfFreedom,
- double denominatorDegreesOfFreedom) {
- super();
- setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
- setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
- }
-
- /**
- * For this distribution, X, this method returns P(X < x).
- *
- * 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>
- * </ul>
- *
- * @param x the value at which the CDF is evaluated.
- * @return CDF for this distribution.
- * @throws MathException if the cumulative probability can not be
- * computed due to convergence or other numerical errors.
- */
- public double cumulativeProbability(double x) throws MathException {
- 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;
- }
-
- /**
- * For this distribution, X, this method returns the critical point x, such
- * that P(X < x) = <code>p</code>.
- * <p>
- * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
- *
- * @param p the desired probability
- * @return x, such that P(X < x) = <code>p</code>
- * @throws MathException if the inverse cumulative probability can not be
- * computed due to convergence or other numerical errors.
- * @throws IllegalArgumentException if <code>p</code> is not a valid
- * probability.
- */
- public double inverseCumulativeProbability(final double p)
- throws MathException {
- if (p == 0) {
- return 0d;
- }
- if (p == 1) {
- return Double.POSITIVE_INFINITY;
- }
- return super.inverseCumulativeProbability(p);
- }
-
- /**
- * Access the domain value lower bound, based on <code>p</code>, used to
- * bracket a CDF root. This method is used by
- * {@link #inverseCumulativeProbability(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;
- }
+ /** The numerator degrees of freedom */
+ private double denominatorDegreesOfFreedom;
- /**
- * Access the domain value upper bound, based on <code>p</code>, used to
- * bracket a CDF root. This method is used by
- * {@link #inverseCumulativeProbability(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;
- }
+ /**
+ * Create a F distribution using the given degrees of freedom.
+ *
+ * @param numeratorDegreesOfFreedom
+ * the numerator degrees of freedom.
+ * @param denominatorDegreesOfFreedom
+ * the denominator degrees of freedom.
+ */
+ public FDistributionImpl(double numeratorDegreesOfFreedom,
+ double denominatorDegreesOfFreedom) {
+ super();
+ setNumeratorDegreesOfFreedom(numeratorDegreesOfFreedom);
+ setDenominatorDegreesOfFreedom(denominatorDegreesOfFreedom);
+ }
- /**
- * Access the initial domain value, based on <code>p</code>, used to
- * bracket a CDF root. This method is used by
- * {@link #inverseCumulativeProbability(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.
- * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
- * positive.
- */
- 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.
- * @throws IllegalArgumentException if <code>degreesOfFreedom</code> is not
- * positive.
- */
- 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;
- }
+ /**
+ * For this distribution, X, this method returns P(X < x).
+ *
+ * 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>
+ * </ul>
+ *
+ * @param x
+ * the value at which the CDF is evaluated.
+ * @return CDF for this distribution.
+ * @throws MathException
+ * if the cumulative probability can not be computed due to
+ * convergence or other numerical errors.
+ */
+ public double cumulativeProbability(double x) throws MathException {
+ 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;
+ }
+
+ /**
+ * For this distribution, X, this method returns the critical point x, such
+ * that P(X < x) = <code>p</code>.
+ * <p>
+ * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.
+ * </p>
+ *
+ * @param p
+ * the desired probability
+ * @return x, such that P(X < x) = <code>p</code>
+ * @throws MathException
+ * if the inverse cumulative probability can not be computed due
+ * to convergence or other numerical errors.
+ * @throws IllegalArgumentException
+ * if <code>p</code> is not a valid probability.
+ */
+ public double inverseCumulativeProbability(final double p)
+ throws MathException {
+ if (p == 0) {
+ return 0d;
+ }
+ if (p == 1) {
+ return Double.POSITIVE_INFINITY;
+ }
+ return super.inverseCumulativeProbability(p);
+ }
+
+ /**
+ * Access the domain value lower bound, based on <code>p</code>, used to
+ * bracket a CDF root. This method is used by
+ * {@link #inverseCumulativeProbability(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 #inverseCumulativeProbability(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 #inverseCumulativeProbability(double)} to find critical values.
+ *
+ * @param p
+ * the desired probability for the critical value
+ * @return initial domain value
+ */
+ protected double getInitialDomain(double p) {
+ double ret = 1.0;
+ double d = getDenominatorDegreesOfFreedom();
+ if (d > 2.0) {
+ // use mean
+ ret = d / (d - 2.0);
+ }
+ return ret;
+ }
+
+ /**
+ * Modify the numerator degrees of freedom.
+ *
+ * @param degreesOfFreedom
+ * the new numerator degrees of freedom.
+ * @throws IllegalArgumentException
+ * if <code>degreesOfFreedom</code> is not positive.
+ */
+ 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.
+ * @throws IllegalArgumentException
+ * if <code>degreesOfFreedom</code> is not positive.
+ */
+ 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/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index cc87af4..009f63f 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -65,6 +65,9 @@
<action dev="brentworden" type="fix" issue="MATH-204" due-to="Mick">
Added root checks for the endpoints.
</action>
+ <action dev="brentworden" type="fix" issue="MATH-227" due-to="Joerg Henning">
+ Fixed F distribution inverse CDF computation for small denominator degrees of freedom.
+ </action>
</release>
<release version="1.2" date="2008-02-24"
description="This release combines bug fixes and new features. Most notable
diff --git a/src/test/org/apache/commons/math/distribution/FDistributionTest.java b/src/test/org/apache/commons/math/distribution/FDistributionTest.java
index d345d68..bd7c143 100644
--- a/src/test/org/apache/commons/math/distribution/FDistributionTest.java
+++ b/src/test/org/apache/commons/math/distribution/FDistributionTest.java
@@ -105,4 +105,18 @@
double x = fd.inverseCumulativeProbability(p);
assertEquals(.999, x, 1.0e-5);
}
+
+ public void testSmallDegreesOfFreedom() throws Exception {
+ org.apache.commons.math.distribution.FDistributionImpl fd =
+ new org.apache.commons.math.distribution.FDistributionImpl(
+ 1.0, 1.0);
+ double p = fd.cumulativeProbability(0.975);
+ double x = fd.inverseCumulativeProbability(p);
+ assertEquals(0.975, x, 1.0e-5);
+
+ fd.setDenominatorDegreesOfFreedom(2.0);
+ p = fd.cumulativeProbability(0.975);
+ x = fd.inverseCumulativeProbability(p);
+ assertEquals(0.975, x, 1.0e-5);
+ }
}
