Skip to content

Commit

Permalink
MATH-1416: Depend on "Commons Numbers" (module "commons-numbers-gamma").
Browse files Browse the repository at this point in the history 
Transitional state (until issue NUMBERS-39 is done).
  • Loading branch information
Gilles committed May 13, 2017
1 parent cec35ba commit e3eff1e
Show file tree
Hide file tree
Showing 16 changed files with 86 additions and 76 deletions.
6 changes: 6 additions & 0 deletions pom.xml
View file
Expand Up @@ -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>
Expand Down
4 changes: 2 additions & 2 deletions src/main/java/org/apache/commons/math4/distribution/BetaDistribution.java
View file
Expand Up @@ -19,7 +19,7 @@
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;
Expand Down Expand Up @@ -74,7 +74,7 @@ public BetaDistribution(double alpha,
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;
}

Expand Down
26 changes: 13 additions & 13 deletions src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java
View file
Expand Up @@ -18,7 +18,8 @@

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;
Expand All @@ -31,6 +32,7 @@
* @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
Expand All @@ -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;
/**
Expand Down Expand Up @@ -136,18 +138,18 @@ public GammaDistribution(double shape,
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);
}

Expand Down Expand Up @@ -222,8 +224,7 @@ public double density(double x) {
*/
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);
}
/*
Expand All @@ -248,8 +249,7 @@ public double logDensity(double x) {
*/
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;
}
/*
Expand Down Expand Up @@ -279,7 +279,7 @@ public double cumulativeProbability(double x) {
if (x <= 0) {
ret = 0;
} else {
ret = Gamma.regularizedGammaP(shape, x / scale);
ret = RegularizedGamma.P.value(shape, x / scale);
}

return ret;
Expand Down
11 changes: 6 additions & 5 deletions src/main/java/org/apache/commons/math4/distribution/NakagamiDistribution.java
View file
Expand Up @@ -19,7 +19,8 @@
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;

/**
Expand Down Expand Up @@ -112,26 +113,26 @@ public double density(double x) {
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);
}

Expand Down
10 changes: 5 additions & 5 deletions src/main/java/org/apache/commons/math4/distribution/PoissonDistribution.java
View file
Expand Up @@ -18,7 +18,7 @@

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;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -166,8 +166,8 @@ public double cumulativeProbability(int x) {
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);
}

/**
Expand Down
4 changes: 2 additions & 2 deletions src/main/java/org/apache/commons/math4/distribution/SaddlePointExpansion.java
View file
Expand Up @@ -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;

Expand Down Expand Up @@ -110,7 +110,7 @@ static double getStirlingError(double z) {
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 {
Expand Down
6 changes: 3 additions & 3 deletions src/main/java/org/apache/commons/math4/distribution/TDistribution.java
View file
Expand Up @@ -19,7 +19,7 @@
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;

/**
Expand Down Expand Up @@ -75,9 +75,9 @@ public TDistribution(double degreesOfFreedom,

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

/**
Expand Down
6 changes: 3 additions & 3 deletions src/main/java/org/apache/commons/math4/distribution/WeibullDistribution.java
View file
Expand Up @@ -20,7 +20,7 @@
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;

/**
Expand Down Expand Up @@ -224,7 +224,7 @@ protected double calculateNumericalMean() {
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)));
}

/**
Expand Down Expand Up @@ -252,7 +252,7 @@ protected double calculateNumericalVariance() {
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);
}

Expand Down
5 changes: 3 additions & 2 deletions src/main/java/org/apache/commons/math4/special/BesselJ.java
View file
Expand Up @@ -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;
Expand Down Expand Up @@ -283,7 +284,7 @@ public static BesselJResult rjBesl(double x, double alpha, int nb) {
}
if (alpha != 0) {
tempa = FastMath.pow(halfx, alpha) /
(alpha * Gamma.gamma(alpha));
(alpha * Gamma.value(alpha));
}
tempb = 0;
if (x + 1 > 1) {
Expand Down Expand Up @@ -620,7 +621,7 @@ public static BesselJResult rjBesl(double x, double alpha, int nb) {
// ---------------------------------------------------------------------

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) {
Expand Down
36 changes: 19 additions & 17 deletions src/main/java/org/apache/commons/math4/special/Beta.java
View file
Expand Up @@ -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;
Expand Down Expand Up @@ -249,11 +251,11 @@ private static double logGammaSum(final double a, final double b)

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

Expand Down Expand Up @@ -415,7 +417,7 @@ public static double logBeta(final double p, final double q) {
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;
Expand All @@ -434,12 +436,12 @@ public static double logBeta(final double p, final double q) {
}
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);
}
}
Expand All @@ -453,29 +455,29 @@ public static double logBeta(final double p, final double q) {
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));
}
}
}
Expand Down

0 comments on commit e3eff1e

Please sign in to comment.