MATH-1416: Depend on "Commons Numbers" (module "commons-numbers-gamma").

Transitional state (until issue NUMBERS-39 is done).
This commit is contained in:
Gilles 2017-05-13 21:21:24 +02:00
parent cec35baf0b
commit e3eff1e3a4
16 changed files with 86 additions and 76 deletions

View File

@ -366,6 +366,12 @@
<version>1.0-SNAPSHOT</version>
</dependency>
<dependency>
<groupId>org.apache.commons</groupId>
<artifactId>commons-numbers-gamma</artifactId>
<version>1.0-SNAPSHOT</version>
</dependency>
<dependency>
<groupId>org.apache.commons</groupId>
<artifactId>commons-rng-client-api</artifactId>

View File

@ -19,7 +19,7 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Beta;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
@ -74,7 +74,7 @@ public class BetaDistribution extends AbstractRealDistribution {
double inverseCumAccuracy) {
this.alpha = alpha;
this.beta = beta;
z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
z = LogGamma.value(alpha) + LogGamma.value(beta) - LogGamma.value(alpha + beta);
solverAbsoluteAccuracy = inverseCumAccuracy;
}

View File

@ -18,7 +18,8 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LanczosApproximation;
import org.apache.commons.numbers.gamma.RegularizedGamma;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
@ -31,6 +32,7 @@ import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGa
* @see <a href="http://mathworld.wolfram.com/GammaDistribution.html">Gamma distribution (MathWorld)</a>
*/
public class GammaDistribution extends AbstractRealDistribution {
private static final double LANCZOS_G = LanczosApproximation.g();
/**
* Default inverse cumulative probability accuracy.
* @since 2.1
@ -44,7 +46,7 @@ public class GammaDistribution extends AbstractRealDistribution {
private final double scale;
/**
* The constant value of {@code shape + g + 0.5}, where {@code g} is the
* Lanczos constant {@link Gamma#LANCZOS_G}.
* Lanczos constant {@link LanczosApproximation#g()}.
*/
private final double shiftedShape;
/**
@ -136,18 +138,18 @@ public class GammaDistribution extends AbstractRealDistribution {
this.shape = shape;
this.scale = scale;
this.solverAbsoluteAccuracy = inverseCumAccuracy;
this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
this.shiftedShape = shape + LANCZOS_G + 0.5;
final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
this.densityPrefactor2 = shape * FastMath.sqrt(aux) / Gamma.lanczos(shape);
this.densityPrefactor2 = shape * FastMath.sqrt(aux) / LanczosApproximation.value(shape);
this.logDensityPrefactor2 = FastMath.log(shape) + 0.5 * FastMath.log(aux) -
FastMath.log(Gamma.lanczos(shape));
FastMath.log(LanczosApproximation.value(shape));
this.densityPrefactor1 = this.densityPrefactor2 / scale *
FastMath.pow(shiftedShape, -shape) *
FastMath.exp(shape + Gamma.LANCZOS_G);
FastMath.exp(shape + LANCZOS_G);
this.logDensityPrefactor1 = this.logDensityPrefactor2 - FastMath.log(scale) -
FastMath.log(shiftedShape) * shape +
shape + Gamma.LANCZOS_G;
this.minY = shape + Gamma.LANCZOS_G - FastMath.log(Double.MAX_VALUE);
shape + LANCZOS_G;
this.minY = shape + LANCZOS_G - FastMath.log(Double.MAX_VALUE);
this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
}
@ -222,8 +224,7 @@ public class GammaDistribution extends AbstractRealDistribution {
*/
final double aux1 = (y - shiftedShape) / shiftedShape;
final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape +
Gamma.LANCZOS_G + aux2;
final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
return densityPrefactor2 / x * FastMath.exp(aux3);
}
/*
@ -248,8 +249,7 @@ public class GammaDistribution extends AbstractRealDistribution {
*/
final double aux1 = (y - shiftedShape) / shiftedShape;
final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape +
Gamma.LANCZOS_G + aux2;
final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
return logDensityPrefactor2 - FastMath.log(x) + aux3;
}
/*
@ -279,7 +279,7 @@ public class GammaDistribution extends AbstractRealDistribution {
if (x <= 0) {
ret = 0;
} else {
ret = Gamma.regularizedGammaP(shape, x / scale);
ret = RegularizedGamma.P.value(shape, x / scale);
}
return ret;

View File

@ -19,7 +19,8 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.Gamma;
import org.apache.commons.numbers.gamma.RegularizedGamma;
import org.apache.commons.math4.util.FastMath;
/**
@ -112,26 +113,26 @@ public class NakagamiDistribution extends AbstractRealDistribution {
if (x <= 0) {
return 0.0;
}
return 2.0 * FastMath.pow(mu, mu) / (Gamma.gamma(mu) * FastMath.pow(omega, mu)) *
return 2.0 * FastMath.pow(mu, mu) / (Gamma.value(mu) * FastMath.pow(omega, mu)) *
FastMath.pow(x, 2 * mu - 1) * FastMath.exp(-mu * x * x / omega);
}
/** {@inheritDoc} */
@Override
public double cumulativeProbability(double x) {
return Gamma.regularizedGammaP(mu, mu * x * x / omega);
return RegularizedGamma.P.value(mu, mu * x * x / omega);
}
/** {@inheritDoc} */
@Override
public double getNumericalMean() {
return Gamma.gamma(mu + 0.5) / Gamma.gamma(mu) * FastMath.sqrt(omega / mu);
return Gamma.value(mu + 0.5) / Gamma.value(mu) * FastMath.sqrt(omega / mu);
}
/** {@inheritDoc} */
@Override
public double getNumericalVariance() {
double v = Gamma.gamma(mu + 0.5) / Gamma.gamma(mu);
double v = Gamma.value(mu + 0.5) / Gamma.value(mu);
return omega * (1 - 1 / mu * v * v);
}

View File

@ -18,7 +18,7 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.RegularizedGamma;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathUtils;
import org.apache.commons.rng.UniformRandomProvider;
@ -52,9 +52,9 @@ public class PoissonDistribution extends AbstractIntegerDistribution {
/**
* Maximum number of iterations for cumulative probability. Cumulative
* probabilities are estimated using either Lanczos series approximation
* of {@link Gamma#regularizedGammaP(double, double, double, int)}
* of {@link RegularizedGamma.P#value(double, double, double, int)}
* or continued fraction approximation of
* {@link Gamma#regularizedGammaQ(double, double, double, int)}.
* {@link RegularizedGamma.Q#value(double, double, double, int)}.
*/
private final int maxIterations;
@ -166,8 +166,8 @@ public class PoissonDistribution extends AbstractIntegerDistribution {
if (x == Integer.MAX_VALUE) {
return 1;
}
return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
maxIterations);
return RegularizedGamma.Q.value((double) x + 1, mean, epsilon,
maxIterations);
}
/**

View File

@ -16,7 +16,7 @@
*/
package org.apache.commons.math4.distribution;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.math4.util.MathUtils;
@ -110,7 +110,7 @@ final class SaddlePointExpansion {
if (FastMath.floor(z2) == z2) {
ret = EXACT_STIRLING_ERRORS[(int) z2];
} else {
ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) +
ret = LogGamma.value(z + 1.0) - (z + 0.5) * FastMath.log(z) +
z - HALF_LOG_2_PI;
}
} else {

View File

@ -19,7 +19,7 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Beta;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.util.FastMath;
/**
@ -75,9 +75,9 @@ public class TDistribution extends AbstractRealDistribution {
final double n = degreesOfFreedom;
final double nPlus1Over2 = (n + 1) / 2;
factor = Gamma.logGamma(nPlus1Over2) -
factor = LogGamma.value(nPlus1Over2) -
0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
Gamma.logGamma(n / 2);
LogGamma.value(n / 2);
}
/**

View File

@ -20,7 +20,7 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.util.FastMath;
/**
@ -224,7 +224,7 @@ public class WeibullDistribution extends AbstractRealDistribution {
final double sh = getShape();
final double sc = getScale();
return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
return sc * FastMath.exp(LogGamma.value(1 + (1 / sh)));
}
/**
@ -252,7 +252,7 @@ public class WeibullDistribution extends AbstractRealDistribution {
final double sc = getScale();
final double mn = getNumericalMean();
return (sc * sc) * FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
return (sc * sc) * FastMath.exp(LogGamma.value(1 + (2 / sh))) -
(mn * mn);
}

View File

@ -17,6 +17,7 @@
package org.apache.commons.math4.special;
import org.apache.commons.numbers.gamma.Gamma;
import org.apache.commons.math4.analysis.UnivariateFunction;
import org.apache.commons.math4.exception.ConvergenceException;
import org.apache.commons.math4.exception.MathIllegalArgumentException;
@ -283,7 +284,7 @@ public class BesselJ
}
if (alpha != 0) {
tempa = FastMath.pow(halfx, alpha) /
(alpha * Gamma.gamma(alpha));
(alpha * Gamma.value(alpha));
}
tempb = 0;
if (x + 1 > 1) {
@ -620,7 +621,7 @@ public class BesselJ
// ---------------------------------------------------------------------
if (FastMath.abs(alpha) > 1e-16) {
sum *= Gamma.gamma(alpha) * FastMath.pow(x * 0.5, -alpha);
sum *= Gamma.value(alpha) * FastMath.pow(x * 0.5, -alpha);
}
tempa = ENMTEN;
if (sum > 1) {

View File

@ -16,6 +16,8 @@
*/
package org.apache.commons.math4.special;
import org.apache.commons.numbers.gamma.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.exception.NumberIsTooSmallException;
import org.apache.commons.math4.exception.OutOfRangeException;
import org.apache.commons.math4.util.ContinuedFraction;
@ -249,11 +251,11 @@ public class Beta {
final double x = (a - 1.0) + (b - 1.0);
if (x <= 0.5) {
return Gamma.logGamma1p(1.0 + x);
return org.apache.commons.math4.special.Gamma.logGamma1p(1.0 + x);
} else if (x <= 1.5) {
return Gamma.logGamma1p(x) + FastMath.log1p(x);
return org.apache.commons.math4.special.Gamma.logGamma1p(x) + FastMath.log1p(x);
} else {
return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
return org.apache.commons.math4.special.Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
}
}
@ -415,7 +417,7 @@ public class Beta {
prod *= ared / (1.0 + ared / b);
}
return (FastMath.log(prod) - n * FastMath.log(b)) +
(Gamma.logGamma(ared) +
(LogGamma.value(ared) +
logGammaMinusLogGammaSum(ared, b));
} else {
double prod1 = 1.0;
@ -434,12 +436,12 @@ public class Beta {
}
return FastMath.log(prod1) +
FastMath.log(prod2) +
(Gamma.logGamma(ared) +
(Gamma.logGamma(bred) -
(LogGamma.value(ared) +
(LogGamma.value(bred) -
logGammaSum(ared, bred)));
} else {
return FastMath.log(prod1) +
Gamma.logGamma(ared) +
LogGamma.value(ared) +
logGammaMinusLogGammaSum(ared, b);
}
}
@ -453,29 +455,29 @@ public class Beta {
prod *= bred / (a + bred);
}
return FastMath.log(prod) +
(Gamma.logGamma(a) +
(Gamma.logGamma(bred) -
(LogGamma.value(a) +
(LogGamma.value(bred) -
logGammaSum(a, bred)));
} else {
return Gamma.logGamma(a) +
return LogGamma.value(a) +
logGammaMinusLogGammaSum(a, b);
}
} else {
return Gamma.logGamma(a) +
Gamma.logGamma(b) -
return LogGamma.value(a) +
LogGamma.value(b) -
logGammaSum(a, b);
}
} else {
if (b >= 10.0) {
return Gamma.logGamma(a) +
return LogGamma.value(a) +
logGammaMinusLogGammaSum(a, b);
} else {
// The following command is the original NSWC implementation.
// return Gamma.logGamma(a) +
// (Gamma.logGamma(b) - Gamma.logGamma(a + b));
// return LogGamma.value(a) +
// (LogGamma.value(b) - LogGamma.value(a + b));
// The following command turns out to be more accurate.
return FastMath.log(Gamma.gamma(a) * Gamma.gamma(b) /
Gamma.gamma(a + b));
return FastMath.log(Gamma.value(a) * Gamma.value(b) /
Gamma.value(a + b));
}
}
}

View File

@ -17,11 +17,11 @@
package org.apache.commons.math4.special;
import org.apache.commons.math4.util.FastMath;
import org.apache.commons.numbers.gamma.RegularizedGamma;
/**
* This is a utility class that provides computation methods related to the
* error functions.
*
*/
public class Erf {
@ -48,7 +48,7 @@ public class Erf {
* <p>erf(x) = 2/&radic;&pi; <sub>0</sub>&int;<sup>x</sup> e<sup>-t<span style="position: relative; top: -.5em">2</span></sup>dt </p>
*
* <p>This implementation computes erf(x) using the
* {@link Gamma#regularizedGammaP(double, double, double, int) regularized gamma function},
* {@link RegularizedGamma.P.value(double, double, double, int) regularized gamma function},
* following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3)</p>
*
* <p>The value returned is always between -1 and 1 (inclusive).
@ -60,13 +60,13 @@ public class Erf {
* @return the error function erf(x)
* @throws org.apache.commons.math4.exception.MaxCountExceededException
* if the algorithm fails to converge.
* @see Gamma#regularizedGammaP(double, double, double, int)
* @see RegularizedGamma.P#value(double, double, double, int)
*/
public static double erf(double x) {
if (FastMath.abs(x) > 40) {
return x > 0 ? 1 : -1;
}
final double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
final double ret = RegularizedGamma.P.value(0.5, x * x, 1.0e-15, 10000);
return x < 0 ? -ret : ret;
}
@ -78,7 +78,7 @@ public class Erf {
* = 1 - {@link #erf(double) erf(x)} </p>
*
* <p>This implementation computes erfc(x) using the
* {@link Gamma#regularizedGammaQ(double, double, double, int) regularized gamma function},
* {@link RegularizedGamma.Q#value(double, double, double, int) regularized gamma function},
* following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3).</p>
*
* <p>The value returned is always between 0 and 2 (inclusive).
@ -90,14 +90,14 @@ public class Erf {
* @return the complementary error function erfc(x)
* @throws org.apache.commons.math4.exception.MaxCountExceededException
* if the algorithm fails to converge.
* @see Gamma#regularizedGammaQ(double, double, double, int)
* @see RegularizedGamma.Q#value(double, double, double, int)
* @since 2.2
*/
public static double erfc(double x) {
if (FastMath.abs(x) > 40) {
return x > 0 ? 0 : 2;
}
final double ret = Gamma.regularizedGammaQ(0.5, x * x, 1.0e-15, 10000);
final double ret = RegularizedGamma.Q.value(0.5, x * x, 1.0e-15, 10000);
return x < 0 ? 2 - ret : ret;
}

View File

@ -24,7 +24,7 @@ import org.apache.commons.math4.exception.MathArithmeticException;
import org.apache.commons.math4.exception.NotPositiveException;
import org.apache.commons.math4.exception.NumberIsTooLargeException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
/**
* Combinatorial utilities.
@ -462,8 +462,8 @@ public final class CombinatoricsUtils {
/**
* Class for computing the natural logarithm of the factorial of {@code n}.
* It allows to allocate a cache of precomputed values.
* In case of cache miss, computation is preformed by a call to
* {@link Gamma#logGamma(double)}.
* In case of cache miss, computation is performed by a call to
* {@link LogGamma#value(double)}.
*/
public static final class FactorialLog {
/**
@ -547,7 +547,7 @@ public final class CombinatoricsUtils {
}
// Delegate.
return Gamma.logGamma(n + 1);
return LogGamma.value(n + 1);
}
}
}

View File

@ -17,7 +17,7 @@
package org.apache.commons.math4.analysis.integration.gauss;
import org.apache.commons.math4.analysis.UnivariateFunction;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.Gamma;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@ -44,7 +44,7 @@ public class LaguerreTest {
final GaussIntegrator integrator = factory.laguerre(7);
final double s = integrator.integrate(f);
Assert.assertEquals(1d, Gamma.gamma(t) / s, tol);
Assert.assertEquals(1d, Gamma.value(t) / s, tol);
}
}
}

View File

@ -24,7 +24,7 @@ import java.io.InputStreamReader;
import org.apache.commons.math4.distribution.GammaDistribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LanczosApproximation;
import org.apache.commons.math4.stat.descriptive.SummaryStatistics;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
@ -189,8 +189,8 @@ public class GammaDistributionTest extends RealDistributionAbstractTest {
if (Double.isNaN(x) || (x <= 0.0)) {
ret = Double.NaN;
} else {
double sum = Gamma.lanczos(x);
double tmp = x + Gamma.LANCZOS_G + .5;
double sum = LanczosApproximation.value(x);
double tmp = x + LanczosApproximation.g() + .5;
ret = ((x + .5) * FastMath.log(tmp)) - tmp +
HALF_LOG_2_PI + FastMath.log(sum / x);
}

View File

@ -19,7 +19,7 @@ package org.apache.commons.math4.distribution;
import org.apache.commons.math4.distribution.WeibullDistribution;
import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.apache.commons.math4.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@ -112,15 +112,15 @@ public class WeibullDistributionTest extends RealDistributionAbstractTest {
dist = new WeibullDistribution(2.5, 3.5);
// In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(LogGamma.value(1 + (1 / 2.5))), tol);
Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) *
FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
FastMath.exp(LogGamma.value(1 + (2 / 2.5))) -
(dist.getNumericalMean() * dist.getNumericalMean()), tol);
dist = new WeibullDistribution(10.4, 2.222);
Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(LogGamma.value(1 + (1 / 10.4))), tol);
Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) *
FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
FastMath.exp(LogGamma.value(1 + (2 / 10.4))) -
(dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
}

View File

@ -17,7 +17,7 @@
package org.apache.commons.math4.util;
import org.apache.commons.math4.exception.NotPositiveException;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.numbers.gamma.LogGamma;
import org.junit.Assert;
import org.junit.Test;
@ -43,9 +43,9 @@ public class FactorialLogTest {
final CombinatoricsUtils.FactorialLog f = CombinatoricsUtils.FactorialLog.create();
// Starting at 21 because for smaller arguments, there is no delegation to the
// "Gamma" class.
// "LogGamma" class.
for (int i = 21; i < 10000; i++) {
final double expected = Gamma.logGamma(i + 1);
final double expected = LogGamma.value(i + 1);
Assert.assertEquals(i + "! ",
expected, f.value(i), 0d);
}