diff --git a/src/main/java/org/apache/commons/math/MessagesResources_fr.java b/src/main/java/org/apache/commons/math/MessagesResources_fr.java
index dceb1af95..3c00fc6f2 100644
--- a/src/main/java/org/apache/commons/math/MessagesResources_fr.java
+++ b/src/main/java/org/apache/commons/math/MessagesResources_fr.java
@@ -215,6 +215,8 @@ public class MessagesResources_fr
"Divergence de fraction continue \u00e0 l''infini pour la valeur {0}" },
{ "Continued fraction convergents failed to converge for value {0}",
"\u00c9chec de convergence de fraction continue pour la valeur {0}" },
+ { "Continued fraction diverged to NaN for value {0}",
+ "Divergence de fraction continue \u00e0 NaN pour la valeur {0}"},
// org.apache.commons.math.util.DefaultTransformer
{ "Conversion Exception in Transformation, Object is null",
@@ -735,6 +737,15 @@ public class MessagesResources_fr
"la borne inf\u00e9rieure ({0}) devrait \u00eatre inf\u00e9rieure " +
"ou \u00e9gale \u00e0 la borne sup\u00e9rieure ({1})" },
+ // org.apache.commons.math.distribution.AbstractContinuousDistribution
+ { "Cumulative probability function returned NaN for argument {0} p = {1}",
+ "Fonction de probabilit\u00e9 cumulative retourn\u00e9 NaN \u00e0 l''argument de {0} p = {1}" },
+
+ // org.apache.commons.math.distribution.AbstractIntegerDistribution
+ { "Discrete cumulative probability function returned NaN for argument {0}",
+ "Discr\u00e8tes fonction de probabilit\u00e9 cumulative retourn\u00e9 NaN \u00e0 l''argument de {0}" },
+
+
// org.apache.commons.math.distribution.BinomialDistributionImpl
{ "number of trials must be non-negative ({0})",
"le nombre d''essais ne doit pas \u00eatre n\u00e9gatif ({0})" },
diff --git a/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java b/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
index e8c27e89c..a2a726e6e 100644
--- a/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
+++ b/src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
@@ -32,6 +32,12 @@ import org.apache.commons.math.analysis.UnivariateRealFunction;
*/
public class BrentSolver extends UnivariateRealSolverImpl {
+ /** Default absolute accuracy */
+ public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
+
+ /** Default maximum number of iterations */
+ public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
+
/** Error message for non-bracketing interval. */
private static final String NON_BRACKETING_MESSAGE =
"function values at endpoints do not have different signs. " +
@@ -51,14 +57,33 @@ public class BrentSolver extends UnivariateRealSolverImpl {
*/
@Deprecated
public BrentSolver(UnivariateRealFunction f) {
- super(f, 100, 1E-6);
+ super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
}
/**
- * Construct a solver.
+ * Construct a solver with default properties.
*/
public BrentSolver() {
- super(100, 1E-6);
+ super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
+ }
+
+ /**
+ * Construct a solver with the given absolute accuracy.
+ *
+ * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
+ */
+ public BrentSolver(double absoluteAccuracy) {
+ super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
+ }
+
+ /**
+ * Contstruct a solver with the given maximum iterations and absolute accuracy.
+ *
+ * @param maximumIterations maximum number of iterations
+ * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
+ */
+ public BrentSolver(int maximumIterations, double absoluteAccuracy) {
+ super(maximumIterations, absoluteAccuracy);
}
/** {@inheritDoc} */
diff --git a/src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java b/src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java
index 815b8b1fd..1e4fccabb 100644
--- a/src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java
+++ b/src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java
@@ -23,6 +23,7 @@ import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.analysis.solvers.BrentSolver;
import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
/**
@@ -39,6 +40,9 @@ public abstract class AbstractContinuousDistribution
/** Serializable version identifier */
private static final long serialVersionUID = -38038050983108802L;
+ /** Solver absolute accuracy for inverse cum computation */
+ private double solverAbsoluteAccuracy = BrentSolver.DEFAULT_ABSOLUTE_ACCURACY;
+
/**
* Default constructor.
*/
@@ -69,11 +73,17 @@ public abstract class AbstractContinuousDistribution
UnivariateRealFunction rootFindingFunction =
new UnivariateRealFunction() {
public double value(double x) throws FunctionEvaluationException {
+ double ret = Double.NaN;
try {
- return cumulativeProbability(x) - p;
+ ret = cumulativeProbability(x) - p;
} catch (MathException ex) {
throw new FunctionEvaluationException(ex, x, ex.getPattern(), ex.getArguments());
}
+ if (Double.isNaN(ret)) {
+ throw new FunctionEvaluationException(x,
+ "Cumulative probability function returned NaN for argument {0} p = {1}", x, p);
+ }
+ return ret;
}
};
@@ -90,9 +100,6 @@ public abstract class AbstractContinuousDistribution
* Check domain endpoints to see if one gives value that is within
* the default solver's defaultAbsoluteAccuracy of 0 (will be the
* case if density has bounded support and p is 0 or 1).
- *
- * TODO: expose the default solver, defaultAbsoluteAccuracy as
- * a constant.
*/
if (Math.abs(rootFindingFunction.value(lowerBound)) < 1E-6) {
return lowerBound;
@@ -106,7 +113,9 @@ public abstract class AbstractContinuousDistribution
// find root
double root = UnivariateRealSolverUtils.solve(rootFindingFunction,
- bracket[0],bracket[1]);
+ // override getSolverAbsoluteAccuracy() to use a Brent solver with
+ // absolute accuracy different from BrentSolver default
+ bracket[0],bracket[1], getSolverAbsoluteAccuracy());
return root;
}
@@ -141,4 +150,13 @@ public abstract class AbstractContinuousDistribution
* P(X < upper bound) > p
*/
protected abstract double getDomainUpperBound(double p);
+
+ /**
+ * Returns the solver absolute accuracy for inverse cum computation.
+ *
+ * @return the maximum absolute error in inverse cumulative probability estimates
+ */
+ protected double getSolverAbsoluteAccuracy() {
+ return solverAbsoluteAccuracy;
+ }
}
diff --git a/src/main/java/org/apache/commons/math/distribution/AbstractIntegerDistribution.java b/src/main/java/org/apache/commons/math/distribution/AbstractIntegerDistribution.java
index 661f1e5b7..284af3864 100644
--- a/src/main/java/org/apache/commons/math/distribution/AbstractIntegerDistribution.java
+++ b/src/main/java/org/apache/commons/math/distribution/AbstractIntegerDistribution.java
@@ -18,6 +18,7 @@ package org.apache.commons.math.distribution;
import java.io.Serializable;
+import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
@@ -173,7 +174,7 @@ public abstract class AbstractIntegerDistribution extends AbstractDistribution
double pm;
while (x0 < x1) {
int xm = x0 + (x1 - x0) / 2;
- pm = cumulativeProbability(xm);
+ pm = checkedCumulativeProbability(xm);
if (pm > p) {
// update x1
if (xm == x1) {
@@ -198,15 +199,39 @@ public abstract class AbstractIntegerDistribution extends AbstractDistribution
}
// insure x0 is the correct critical point
- pm = cumulativeProbability(x0);
+ pm = checkedCumulativeProbability(x0);
while (pm > p) {
--x0;
- pm = cumulativeProbability(x0);
+ pm = checkedCumulativeProbability(x0);
}
return x0;
}
+ /**
+ * Computes the cumulative probablity function and checks for NaN values returned.
+ * Throws MathException if the value is NaN. Wraps and rethrows any MathException encountered
+ * evaluating the cumulative probability function in a FunctionEvaluationException. Throws
+ * FunctionEvaluationException of the cumulative probability function returns NaN.
+ *
+ * @param argument input value
+ * @return cumulative probability
+ * @throws FunctionEvaluationException if a MathException occurs computing the cumulative probability
+ */
+ private double checkedCumulativeProbability(int argument) throws FunctionEvaluationException {
+ double result = Double.NaN;
+ try {
+ result = cumulativeProbability(argument);
+ } catch (MathException ex) {
+ throw new FunctionEvaluationException(ex, argument, ex.getPattern(), ex.getArguments());
+ }
+ if (Double.isNaN(result)) {
+ throw new FunctionEvaluationException(argument,
+ "Discrete cumulative probability function returned NaN for argument {0}", argument);
+ }
+ return result;
+ }
+
/**
* Access the domain value lower bound, based on p
, used to
* bracket a PDF root. This method is used by
diff --git a/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java
index 010060c48..e5bdce781 100644
--- a/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java
+++ b/src/main/java/org/apache/commons/math/distribution/NormalDistributionImpl.java
@@ -33,6 +33,9 @@ import org.apache.commons.math.special.Erf;
public class NormalDistributionImpl extends AbstractContinuousDistribution
implements NormalDistribution, Serializable {
+ /** Default inverse cumulative probability accuracy */
+ public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
+
/** Serializable version identifier */
private static final long serialVersionUID = 8589540077390120676L;
@@ -45,15 +48,31 @@ public class NormalDistributionImpl extends AbstractContinuousDistribution
/** The standard deviation of this distribution. */
private double standardDeviation = 1;
+ /** Inverse cumulative probability accuracy */
+ private final double solverAbsoluteAccuracy;
+
/**
* Create a normal distribution using the given mean and standard deviation.
* @param mean mean for this distribution
* @param sd standard deviation for this distribution
*/
public NormalDistributionImpl(double mean, double sd){
+ this(mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
+ }
+
+ /**
+ * Create a normal distribution using the given mean, standard deviation and
+ * inverse cumulative distribution accuracy.
+ *
+ * @param mean mean for this distribution
+ * @param sd standard deviation for this distribution
+ * @param inverseCumAccuracy inverse cumulative probability accuracy
+ */
+ public NormalDistributionImpl(double mean, double sd, double inverseCumAccuracy) {
super();
- setMean(mean);
- setStandardDeviation(sd);
+ this.mean = mean;
+ this.standardDeviation = sd;
+ solverAbsoluteAccuracy = inverseCumAccuracy;
}
/**
@@ -136,6 +155,17 @@ public class NormalDistributionImpl extends AbstractContinuousDistribution
}
}
+ /**
+ * Return the absolute accuracy setting of the solver used to estimate
+ * inverse cumulative probabilities.
+ *
+ * @return the solver absolute accuracy
+ */
+ @Override
+ protected double getSolverAbsoluteAccuracy() {
+ return solverAbsoluteAccuracy;
+ }
+
/**
* For this distribution, X, this method returns the critical point x, such
* that P(X < x) = p
.
diff --git a/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java b/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
index 5e09dba43..9e2509978 100644
--- a/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
+++ b/src/main/java/org/apache/commons/math/distribution/PoissonDistributionImpl.java
@@ -31,6 +31,16 @@ import org.apache.commons.math.util.MathUtils;
public class PoissonDistributionImpl extends AbstractIntegerDistribution
implements PoissonDistribution, Serializable {
+ /**
+ * Default maximum number of iterations for cumulative probability calculations.
+ */
+ public static final int DEFAULT_MAX_ITERATIONS = 10000000;
+
+ /**
+ * Default convergence criterion
+ */
+ public static final double DEFAULT_EPSILON = 1E-12;
+
/** Serializable version identifier */
private static final long serialVersionUID = -3349935121172596109L;
@@ -42,6 +52,19 @@ public class PoissonDistributionImpl extends AbstractIntegerDistribution
*/
private double mean;
+ /**
+ * Maximum number of iterations for cumulative probability.
+ *
+ * Cumulative probabilities are estimated using either Lanczos series approximation of
+ * Gamma#regularizedGammaP or continued fraction approximation of Gamma#regularizedGammaQ.
+ */
+ private int maxIterations = DEFAULT_MAX_ITERATIONS;
+
+ /**
+ * Convergence criterion for cumulative probability.
+ */
+ private double epsilon = DEFAULT_EPSILON;
+
/**
* Create a new Poisson distribution with the given the mean. The mean value
* must be positive; otherwise an IllegalArgument
is thrown.
@@ -53,6 +76,43 @@ public class PoissonDistributionImpl extends AbstractIntegerDistribution
this(p, new NormalDistributionImpl());
}
+ /**
+ * Create a new Poisson distribution with the given mean, convergence criterion
+ * and maximum number of iterations.
+ *
+ * @param p the Poisson mean
+ * @param epsilon the convergence criteria for cumulative probabilites
+ * @param maxIterations the maximum number of iterations for cumulative probabilites
+ */
+ public PoissonDistributionImpl(double p, double epsilon, int maxIterations) {
+ setMean(p);
+ this.epsilon = epsilon;
+ this.maxIterations = maxIterations;
+ }
+
+ /**
+ * Create a new Poisson distribution with the given mean and convergence criterion.
+ *
+ * @param p the Poisson mean
+ * @param epsilon the convergence criteria for cumulative probabilites
+ */
+ public PoissonDistributionImpl(double p, double epsilon) {
+ setMean(p);
+ this.epsilon = epsilon;
+ }
+
+ /**
+ * Create a new Poisson distribution with the given mean and maximum number of iterations.
+ *
+ * @param p the Poisson mean
+ * @param maxIterations the maximum number of iterations for cumulative probabilites
+ */
+ public PoissonDistributionImpl(double p, int maxIterations) {
+ setMean(p);
+ this.maxIterations = maxIterations;
+ }
+
+
/**
* Create a new Poisson distribution with the given the mean. The mean value
* must be positive; otherwise an IllegalArgument
is thrown.
@@ -132,8 +192,7 @@ public class PoissonDistributionImpl extends AbstractIntegerDistribution
if (x == Integer.MAX_VALUE) {
return 1;
}
- return Gamma.regularizedGammaQ((double) x + 1, mean, 1E-12,
- Integer.MAX_VALUE);
+ return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon, maxIterations);
}
/**
diff --git a/src/main/java/org/apache/commons/math/special/Gamma.java b/src/main/java/org/apache/commons/math/special/Gamma.java
index a80c01948..af3086b1e 100644
--- a/src/main/java/org/apache/commons/math/special/Gamma.java
+++ b/src/main/java/org/apache/commons/math/special/Gamma.java
@@ -166,7 +166,7 @@ public class Gamma {
ret = Double.NaN;
} else if (x == 0.0) {
ret = 0.0;
- } else if (a >= 1.0 && x > a) {
+ } else if (x >= a + 1) {
// use regularizedGammaQ because it should converge faster in this
// case.
ret = 1.0 - regularizedGammaQ(a, x, epsilon, maxIterations);
@@ -175,7 +175,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 < maxIterations) {
+ while (Math.abs(an/sum) > epsilon && n < maxIterations && sum < Double.POSITIVE_INFINITY) {
// compute next element in the series
n = n + 1.0;
an = an * (x / (a + n));
@@ -185,6 +185,8 @@ public class Gamma {
}
if (n >= maxIterations) {
throw new MaxIterationsExceededException(maxIterations);
+ } else if (Double.isInfinite(sum)) {
+ ret = 1.0;
} else {
ret = Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
}
@@ -216,7 +218,7 @@ public class Gamma {
*
* Regularized Gamma Function, equation (1).
*