Author: psteitz
Date: Sun Oct 20 20:42:41 2013
New Revision: 1533974
URL: http://svn.apache.org/r1533974
Log:
Added logDensity methods to AbstractReal/IntegerDistribution with naive default
implementations and improved implementations for some current distributions.
JIRA: MATH-1039
Patch provided by Aleksei Dievskii
Modified:
commons/proper/math/trunk/pom.xml
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BinomialDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ChiSquaredDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/FDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LogNormalDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/NormalDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ParetoDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PoissonDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/TDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/WeibullDistribution.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ZipfDistribution.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/IntegerDistributionAbstractTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/PoissonDistributionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/RealDistributionAbstractTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/ZipfDistributionTest.java
Modified: commons/proper/math/trunk/pom.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
--- commons/proper/math/trunk/pom.xml (original)
+++ commons/proper/math/trunk/pom.xml Sun Oct 20 20:42:41 2013
@@ -181,6 +181,9 @@
<name>Larry Diamond</name>
</contributor>
<contributor>
+ <name>Aleksei Dievskii</name>
+ </contributor>
+ <contributor>
<name>Rodrigo di Lorenzo Lopes</name>
</contributor>
<contributor>
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Sun Oct 20 20:42:41 2013
@@ -51,6 +51,10 @@ If the output is not quite correct, chec
</properties>
<body>
<release version="x.y" date="TBD" description="TBD">
+ <action dev="psteitz" type="update" issue="MATH-1039" due-to="Aleksei
Dievskii">
+ Added logDensity methods to AbstractReal/IntegerDistribution with
naive default
+ implementations and improved implementations for some current
distributions.
+ </action>
<action dev="psteitz" type="add" issue="MATH-1038" due-to="Thorsten
Schaefer">
Added ConfidenceInterval class and BinomialConfidenceInterval
providing several
estimators for confidence intervals for binomial probabilities.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractIntegerDistribution.java
Sun Oct 20 20:42:41 2013
@@ -232,4 +232,23 @@ public abstract class AbstractIntegerDis
}
return result;
}
+
+ /**
+ * For a random variable {@code X} whose values are distributed according
to
+ * this distribution, this method returns {@code log(P(X = x))}, where
+ * {@code log} is the natural logarithm. In other words, this method
+ * represents the logarithm of the probability mass function (PMF) for the
+ * distribution. Note that due to the floating point precision and
+ * under/overflow issues, this method will for some distributions be more
+ * precise and faster than computing the logarithm of
+ * {@link #probability(int)}.
+ * <p>
+ * The default implementation simply computes the logarithm of {@code
probability(x)}.</p>
+ *
+ * @param x the point at which the PMF is evaluated
+ * @return the logarithm of the value of the probability mass function at
{@code x}
+ */
+ public double logProbability(int x) {
+ return FastMath.log(probability(x));
+ }
}
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/AbstractRealDistribution.java
Sun Oct 20 20:42:41 2013
@@ -286,5 +286,23 @@ implements RealDistribution, Serializabl
public double probability(double x) {
return 0d;
}
+
+ /**
+ * Returns the natural logarithm of the probability density function (PDF)
of this distribution
+ * evaluated at the specified point {@code x}. In general, the PDF is the
derivative of the
+ * {@link #cumulativeProbability(double) CDF}. If the derivative does not
exist at {@code x},
+ * then an appropriate replacement should be returned, e.g. {@code
Double.POSITIVE_INFINITY},
+ * {@code Double.NaN}, or the limit inferior or limit superior of the
difference quotient. Note
+ * that due to the floating point precision and under/overflow issues,
this method will for some
+ * distributions be more precise and faster than computing the logarithm of
+ * {@link #density(double)}. The default implementation simply computes
the logarithm of
+ * {@code density(x)}.
+ *
+ * @param x the point at which the PDF is evaluated
+ * @return the logarithm of the value of the probability density function
at point {@code x}
+ */
+ public double logDensity(double x) {
+ return FastMath.log(density(x));
+ }
}
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java
Sun Oct 20 20:42:41 2013
@@ -156,6 +156,29 @@ public class BetaDistribution extends Ab
}
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ recomputeZ();
+ if (x < 0 || x > 1) {
+ return Double.NEGATIVE_INFINITY;
+ } else if (x == 0) {
+ if (alpha < 1) {
+ throw new
NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA,
alpha, 1, false);
+ }
+ return Double.NEGATIVE_INFINITY;
+ } else if (x == 1) {
+ if (beta < 1) {
+ throw new
NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA,
beta, 1, false);
+ }
+ return Double.NEGATIVE_INFINITY;
+ } else {
+ double logX = FastMath.log(x);
+ double log1mX = FastMath.log1p(-x);
+ return (alpha - 1) * logX + (beta - 1) * log1mX - z;
+ }
+ }
+
/** {@inheritDoc} */
public double cumulativeProbability(double x) {
if (x <= 0) {
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BinomialDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BinomialDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BinomialDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/BinomialDistribution.java
Sun Oct 20 20:42:41 2013
@@ -110,6 +110,20 @@ public class BinomialDistribution extend
return ret;
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logProbability(int x) {
+ double ret;
+ if (x < 0 || x > numberOfTrials) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else {
+ ret = SaddlePointExpansion.logBinomialProbability(x,
+ numberOfTrials, probabilityOfSuccess,
+ 1.0 - probabilityOfSuccess);
+ }
+ return ret;
+ }
+
/** {@inheritDoc} */
public double cumulativeProbability(int x) {
double ret;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ChiSquaredDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ChiSquaredDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ChiSquaredDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ChiSquaredDistribution.java
Sun Oct 20 20:42:41 2013
@@ -108,6 +108,12 @@ public class ChiSquaredDistribution exte
return gamma.density(x);
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ return gamma.logDensity(x);
+ }
+
/** {@inheritDoc} */
public double cumulativeProbability(double x) {
return gamma.cumulativeProbability(x);
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java
Sun Oct 20 20:42:41 2013
@@ -56,6 +56,8 @@ public class ExponentialDistribution ext
private static final double[] EXPONENTIAL_SA_QI;
/** The mean of this distribution. */
private final double mean;
+ /** The logarithm of the mean, stored to reduce computing time. **/
+ private final double logMean;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
@@ -144,6 +146,7 @@ public class ExponentialDistribution ext
throw new NotStrictlyPositiveException(LocalizedFormats.MEAN,
mean);
}
this.mean = mean;
+ logMean = FastMath.log(mean);
solverAbsoluteAccuracy = inverseCumAccuracy;
}
@@ -164,6 +167,12 @@ public class ExponentialDistribution ext
return FastMath.exp(-x / mean) / mean;
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ return -x / mean - logMean;
+ }
+
/**
* {@inheritDoc}
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/FDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/FDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/FDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/FDistribution.java
Sun Oct 20 20:42:41 2013
@@ -154,6 +154,21 @@ public class FDistribution extends Abstr
Beta.logBeta(nhalf, mhalf));
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ final double nhalf = numeratorDegreesOfFreedom / 2;
+ final double mhalf = denominatorDegreesOfFreedom / 2;
+ final double logx = FastMath.log(x);
+ final double logn = FastMath.log(numeratorDegreesOfFreedom);
+ final double logm = FastMath.log(denominatorDegreesOfFreedom);
+ final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
+ denominatorDegreesOfFreedom);
+ return nhalf * logn + nhalf * logx - logx +
+ mhalf * logm - nhalf * lognxm - mhalf * lognxm -
+ Beta.logBeta(nhalf, mhalf);
+ }
+
/**
* {@inheritDoc}
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java
Sun Oct 20 20:42:41 2013
@@ -58,6 +58,15 @@ public class GammaDistribution extends A
private final double densityPrefactor1;
/**
* The constant value of
+ * {@code log(shape / scale * sqrt(e / (2 * pi * (shape + g + 0.5))) /
L(shape))},
+ * where {@code L(shape)} is the Lanczos approximation returned by
+ * {@link Gamma#lanczos(double)}. This prefactor is used in
+ * {@link #logDensity(double)}, when no overflow occurs with the natural
+ * calculation.
+ */
+ private final double logDensityPrefactor1;
+ /**
+ * The constant value of
* {@code shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape)},
* where {@code L(shape)} is the Lanczos approximation returned by
* {@link Gamma#lanczos(double)}. This prefactor is used in
@@ -66,6 +75,15 @@ public class GammaDistribution extends A
*/
private final double densityPrefactor2;
/**
+ * The constant value of
+ * {@code log(shape * sqrt(e / (2 * pi * (shape + g + 0.5))) / L(shape))},
+ * where {@code L(shape)} is the Lanczos approximation returned by
+ * {@link Gamma#lanczos(double)}. This prefactor is used in
+ * {@link #logDensity(double)}, when overflow occurs with the natural
+ * calculation.
+ */
+ private final double logDensityPrefactor2;
+ /**
* Lower bound on {@code y = x / scale} for the selection of the
computation
* method in {@link #density(double)}. For {@code y <= minY}, the natural
* calculation overflows.
@@ -159,9 +177,14 @@ public class GammaDistribution extends A
this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
this.densityPrefactor2 = shape * FastMath.sqrt(aux) /
Gamma.lanczos(shape);
+ this.logDensityPrefactor2 = FastMath.log(shape) + 0.5 *
FastMath.log(aux) -
+ FastMath.log(Gamma.lanczos(shape));
this.densityPrefactor1 = this.densityPrefactor2 / scale *
FastMath.pow(shiftedShape, -shape) *
FastMath.exp(shape + Gamma.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);
this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
}
@@ -271,6 +294,32 @@ public class GammaDistribution extends A
FastMath.pow(y, shape - 1);
}
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ /* see the comment in {@link #density(double)} for computation details
+ */
+ if (x < 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ final double y = x / scale;
+ if ((y <= minY) || (FastMath.log(y) >= maxLogY)) {
+ /*
+ * Overflow.
+ */
+ 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;
+ return logDensityPrefactor2 - FastMath.log(x) + aux3;
+ }
+ /*
+ * Natural calculation.
+ */
+ return logDensityPrefactor1 - y +
+ FastMath.log(y) * (shape - 1);
+ }
+
/**
* {@inheritDoc}
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java
Sun Oct 20 20:42:41 2013
@@ -86,6 +86,19 @@ public class GeometricDistribution exten
}
/** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ double ret;
+ if (x < 0) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else {
+ final double p = probabilityOfSuccess;
+ ret = x * FastMath.log1p(-p) + FastMath.log(p);
+ }
+ return ret;
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(int x) {
double ret;
if (x < 0) {
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java
Sun Oct 20 20:42:41 2013
@@ -214,6 +214,30 @@ public class HypergeometricDistribution
return ret;
}
+ /** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ double ret;
+
+ int[] domain = getDomain(populationSize, numberOfSuccesses,
sampleSize);
+ if (x < domain[0] || x > domain[1]) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else {
+ double p = (double) sampleSize / (double) populationSize;
+ double q = (double) (populationSize - sampleSize) / (double)
populationSize;
+ double p1 = SaddlePointExpansion.logBinomialProbability(x,
+ numberOfSuccesses, p, q);
+ double p2 =
+ SaddlePointExpansion.logBinomialProbability(sampleSize - x,
+ populationSize - numberOfSuccesses, p, q);
+ double p3 =
+ SaddlePointExpansion.logBinomialProbability(sampleSize,
populationSize, p, q);
+ ret = p1 + p2 - p3;
+ }
+
+ return ret;
+ }
+
/**
* For this distribution, {@code X}, this method returns {@code P(X >= x)}.
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
Sun Oct 20 20:42:41 2013
@@ -80,6 +80,21 @@ public class LevyDistribution extends Ab
}
/** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+
+ final double delta = x - mu;
+ final double f = halfC / delta;
+ return 0.5 * FastMath.log(f / FastMath.PI) - f - FastMath.log(delta);
+ }
+
+ /** {@inheritDoc}
* <p>
* From Wikipedia: the cumulative distribution function is
* </p>
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LogNormalDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LogNormalDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LogNormalDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/LogNormalDistribution.java
Sun Oct 20 20:42:41 2013
@@ -71,6 +71,8 @@ public class LogNormalDistribution exten
/** The shape parameter of this distribution. */
private final double shape;
+ /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster
computation. */
+ private final double logShapePlusHalfLog2Pi;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
@@ -149,6 +151,7 @@ public class LogNormalDistribution exten
this.scale = scale;
this.shape = shape;
+ this.logShapePlusHalfLog2Pi = FastMath.log(shape) + 0.5 *
FastMath.log(2 * FastMath.PI);
this.solverAbsoluteAccuracy = inverseCumAccuracy;
}
@@ -190,6 +193,21 @@ public class LogNormalDistribution exten
return FastMath.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
}
+ /** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x <= 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ final double logX = FastMath.log(x);
+ final double x0 = logX - scale;
+ final double x1 = x0 / shape;
+ return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX);
+ }
+
/**
* {@inheritDoc}
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/NormalDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/NormalDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/NormalDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/NormalDistribution.java
Sun Oct 20 20:42:41 2013
@@ -49,6 +49,8 @@ public class NormalDistribution extends
private final double mean;
/** Standard deviation of this distribution. */
private final double standardDeviation;
+ /** The value of {@code log(sd) + 0.5*log(2*pi)} stored for faster
computation. */
+ private final double logStandardDeviationPlusHalfLog2Pi;
/** Inverse cumulative probability accuracy. */
private final double solverAbsoluteAccuracy;
@@ -124,6 +126,7 @@ public class NormalDistribution extends
this.mean = mean;
standardDeviation = sd;
+ logStandardDeviationPlusHalfLog2Pi = FastMath.log(sd) + 0.5 *
FastMath.log(2 * FastMath.PI);
solverAbsoluteAccuracy = inverseCumAccuracy;
}
@@ -152,6 +155,13 @@ public class NormalDistribution extends
return FastMath.exp(-0.5 * x1 * x1) / (standardDeviation * SQRT2PI);
}
+ /** {@inheritDoc} */
+ public double logDensity(double x) {
+ final double x0 = x - mean;
+ final double x1 = x0 / standardDeviation;
+ return -0.5 * x1 * x1 - logStandardDeviationPlusHalfLog2Pi;
+ }
+
/**
* {@inheritDoc}
*
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ParetoDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ParetoDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ParetoDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ParetoDistribution.java
Sun Oct 20 20:42:41 2013
@@ -174,6 +174,18 @@ public class ParetoDistribution extends
return FastMath.pow(scale, shape) / FastMath.pow(x, shape + 1) * shape;
}
+ /** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x < scale) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ return FastMath.log(scale) * shape - FastMath.log(x) * (shape + 1) +
FastMath.log(shape);
+ }
+
/**
* {@inheritDoc}
* <p>
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PascalDistribution.java
Sun Oct 20 20:42:41 2013
@@ -68,6 +68,12 @@ public class PascalDistribution extends
private final int numberOfSuccesses;
/** The probability of success. */
private final double probabilityOfSuccess;
+ /** The value of {@code log(p)}, where {@code p} is the probability of
success,
+ * stored for faster computation. */
+ private final double logProbabilityOfSuccess;
+ /** The value of {@code log(1-p)}, where {@code p} is the probability of
success,
+ * stored for faster computation. */
+ private final double log1mProbabilityOfSuccess;
/**
* Create a Pascal distribution with the given number of successes and
@@ -112,6 +118,8 @@ public class PascalDistribution extends
numberOfSuccesses = r;
probabilityOfSuccess = p;
+ logProbabilityOfSuccess = FastMath.log(p);
+ log1mProbabilityOfSuccess = FastMath.log1p(-p);
}
/**
@@ -147,6 +155,21 @@ public class PascalDistribution extends
}
/** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ double ret;
+ if (x < 0) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else {
+ ret = CombinatoricsUtils.binomialCoefficientLog(x +
+ numberOfSuccesses - 1, numberOfSuccesses - 1) +
+ logProbabilityOfSuccess * numberOfSuccesses +
+ log1mProbabilityOfSuccess * x;
+ }
+ return ret;
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(int x) {
double ret;
if (x < 0) {
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PoissonDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PoissonDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PoissonDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/PoissonDistribution.java
Sun Oct 20 20:42:41 2013
@@ -175,6 +175,22 @@ public class PoissonDistribution extends
}
/** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ double ret;
+ if (x < 0 || x == Integer.MAX_VALUE) {
+ ret = Double.NEGATIVE_INFINITY;
+ } else if (x == 0) {
+ ret = -mean;
+ } else {
+ ret = -SaddlePointExpansion.getStirlingError(x) -
+ SaddlePointExpansion.getDeviancePart(x, mean) -
+ 0.5 * FastMath.log(MathUtils.TWO_PI) - 0.5 * FastMath.log(x);
+ }
+ return ret;
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(int x) {
if (x < 0) {
return 0;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/TDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/TDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/TDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/TDistribution.java
Sun Oct 20 20:42:41 2013
@@ -130,6 +130,18 @@ public class TDistribution extends Abstr
}
/** {@inheritDoc} */
+ @Override
+ public double logDensity(double x) {
+ final double n = degreesOfFreedom;
+ final double nPlus1Over2 = (n + 1) / 2;
+ return Gamma.logGamma(nPlus1Over2) -
+ 0.5 * (FastMath.log(FastMath.PI) +
+ FastMath.log(n)) -
+ Gamma.logGamma(n / 2) -
+ nPlus1Over2 * FastMath.log(1 + x * x / n);
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(double x) {
double ret;
if (x == 0) {
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/WeibullDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/WeibullDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/WeibullDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/WeibullDistribution.java
Sun Oct 20 20:42:41 2013
@@ -175,6 +175,26 @@ public class WeibullDistribution extends
}
/** {@inheritDoc} */
+ @Override
+ public double logDensity(double x) {
+ if (x < 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ final double xscale = x / scale;
+ final double logxscalepow = FastMath.log(xscale) * (shape - 1);
+
+ /*
+ * FastMath.pow(x / scale, shape) =
+ * FastMath.pow(xscale, shape) =
+ * FastMath.pow(xscale, shape - 1) * xscale
+ */
+ final double xscalepowshape = FastMath.exp(logxscalepow) * xscale;
+
+ return FastMath.log(shape / scale) + logxscalepow - xscalepowshape;
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(double x) {
double ret;
if (x <= 0.0) {
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ZipfDistribution.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ZipfDistribution.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ZipfDistribution.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/distribution/ZipfDistribution.java
Sun Oct 20 20:42:41 2013
@@ -115,6 +115,16 @@ public class ZipfDistribution extends Ab
}
/** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ if (x <= 0 || x > numberOfElements) {
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ return -FastMath.log(x) * exponent -
FastMath.log(generalizedHarmonic(numberOfElements, exponent));
+ }
+
+ /** {@inheritDoc} */
public double cumulativeProbability(final int x) {
if (x <= 0) {
return 0.0;
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/IntegerDistributionAbstractTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/IntegerDistributionAbstractTest.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/IntegerDistributionAbstractTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/IntegerDistributionAbstractTest.java
Sun Oct 20 20:42:41 2013
@@ -18,6 +18,7 @@ package org.apache.commons.math3.distrib
import org.apache.commons.math3.TestUtils;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
+import org.apache.commons.math3.util.FastMath;
import org.junit.After;
import org.junit.Assert;
import org.junit.Before;
@@ -60,6 +61,9 @@ public abstract class IntegerDistributio
/** Values used to test probability density calculations */
private double[] densityTestValues;
+ /** Values used to test logarithmic probability density calculations */
+ private double[] logDensityTestValues;
+
/** Arguments used to test cumulative probability density calculations */
private int[] cumulativeTestPoints;
@@ -83,6 +87,22 @@ public abstract class IntegerDistributio
/** Creates the default probability density test expected values */
public abstract double[] makeDensityTestValues();
+ /** Creates the default logarithmic probability density test expected
values.
+ *
+ * The default implementation simply computes the logarithm of all the
values in
+ * {@link #makeDensityTestValues()}.
+ *
+ * @return double[] the default logarithmic probability density test
expected values.
+ */
+ public double[] makeLogDensityTestValues() {
+ final double[] densityTestValues = makeDensityTestValues();
+ final double[] logDensityTestValues = new
double[densityTestValues.length];
+ for (int i = 0; i < densityTestValues.length; i++) {
+ logDensityTestValues[i] = FastMath.log(densityTestValues[i]);
+ }
+ return logDensityTestValues;
+ }
+
/** Creates the default cumulative probability density test input values */
public abstract int[] makeCumulativeTestPoints();
@@ -105,6 +125,7 @@ public abstract class IntegerDistributio
distribution = makeDistribution();
densityTestPoints = makeDensityTestPoints();
densityTestValues = makeDensityTestValues();
+ logDensityTestValues = makeLogDensityTestValues();
cumulativeTestPoints = makeCumulativeTestPoints();
cumulativeTestValues = makeCumulativeTestValues();
inverseCumulativeTestPoints = makeInverseCumulativeTestPoints();
@@ -119,6 +140,7 @@ public abstract class IntegerDistributio
distribution = null;
densityTestPoints = null;
densityTestValues = null;
+ logDensityTestValues = null;
cumulativeTestPoints = null;
cumulativeTestValues = null;
inverseCumulativeTestPoints = null;
@@ -140,6 +162,19 @@ public abstract class IntegerDistributio
}
/**
+ * Verifies that logarithmic probability density calculations match
expected values
+ * using current test instance data.
+ */
+ protected void verifyLogDensities() {
+ for (int i = 0; i < densityTestPoints.length; i++) {
+ // FIXME: when logProbability methods are added to
IntegerDistribution in 4.0, remove cast below
+ Assert.assertEquals("Incorrect log density value returned for " +
densityTestPoints[i],
+ logDensityTestValues[i],
+ ((AbstractIntegerDistribution)
distribution).logProbability(densityTestPoints[i]), tolerance);
+ }
+ }
+
+ /**
* Verifies that cumulative probability density calculations match
expected values
* using current test instance data
*/
@@ -176,6 +211,15 @@ public abstract class IntegerDistributio
}
/**
+ * Verifies that logarithmic probability density calculations match
expected values
+ * using default test instance data
+ */
+ @Test
+ public void testLogDensities() {
+ verifyLogDensities();
+ }
+
+ /**
* Verifies that cumulative probability density calculations match
expected values
* using default test instance data
*/
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
Sun Oct 20 20:42:41 2013
@@ -68,4 +68,15 @@ public class LevyDistributionTest extend
};
}
+ /**
+ * Creates the default logarithmic probability density test expected
values.
+ * Reference values are from R, version 2.14.1.
+ */
+ @Override
+ public double[] makeLogDensityTestValues() {
+ return new double[] {
+ -1987.561573341398d, -14.469328620160d, -3.843764717971d,
+ -0.883485488811d, 0.076793740349d, -1.127785768948d,
+ -2.650679030597d, -3.644945255983d};
+ }
}
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/PoissonDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/PoissonDistributionTest.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/PoissonDistributionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/PoissonDistributionTest.java
Sun Oct 20 20:42:41 2013
@@ -67,6 +67,18 @@ public class PoissonDistributionTest ext
0.156293451851d, 0.00529247667642d, 8.27746364655e-09};
}
+ /**
+ * Creates the default logarithmic probability density test expected
values.
+ * Reference values are from R, version 2.14.1.
+ */
+ @Override
+ public double[] makeLogDensityTestValues() {
+ return new double[] { Double.NEGATIVE_INFINITY, -4.000000000000d,
+ -2.613705638880d, -1.920558458320d, -1.632876385868d,
+ -1.632876385868d, -1.856019937183d, -5.241468961877d,
+ -18.609729238356d};
+ }
+
/**
* Creates the default cumulative probability density test input values.
*/
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/RealDistributionAbstractTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/RealDistributionAbstractTest.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/RealDistributionAbstractTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/RealDistributionAbstractTest.java
Sun Oct 20 20:42:41 2013
@@ -95,6 +95,9 @@ public abstract class RealDistributionAb
/** Values used to test density calculations */
private double[] densityTestValues;
+ /** Values used to test logarithmic density calculations */
+ private double[] logDensityTestValues;
+
//-------------------- Abstract methods -----------------------------------
/** Creates the default continuous distribution instance to use in tests.
*/
@@ -109,6 +112,18 @@ public abstract class RealDistributionAb
/** Creates the default density test expected values */
public abstract double[] makeDensityTestValues();
+ /** Creates the default logarithmic density test expected values.
+ * The default implementation simply computes the logarithm
+ * of each value returned by {@link #makeDensityTestValues()}.*/
+ public double[] makeLogDensityTestValues() {
+ final double[] densityTestValues = makeDensityTestValues();
+ final double[] logDensityTestValues = new
double[densityTestValues.length];
+ for (int i = 0; i < densityTestValues.length; i++) {
+ logDensityTestValues[i] = FastMath.log(densityTestValues[i]);
+ }
+ return logDensityTestValues;
+ }
+
//---- Default implementations of inverse test data generation methods ----
/** Creates the default inverse cumulative probability test input values */
@@ -134,6 +149,7 @@ public abstract class RealDistributionAb
inverseCumulativeTestPoints = makeInverseCumulativeTestPoints();
inverseCumulativeTestValues = makeInverseCumulativeTestValues();
densityTestValues = makeDensityTestValues();
+ logDensityTestValues = makeLogDensityTestValues();
}
/**
@@ -147,6 +163,7 @@ public abstract class RealDistributionAb
inverseCumulativeTestPoints = null;
inverseCumulativeTestValues = null;
densityTestValues = null;
+ logDensityTestValues = null;
}
//-------------------- Verification methods -------------------------------
@@ -208,6 +225,19 @@ public abstract class RealDistributionAb
}
}
+ /**
+ * Verifies that logarithmic density calculations match expected values
+ */
+ protected void verifyLogDensities() {
+ // FIXME: when logProbability methods are added to RealDistribution in
4.0, remove cast below
+ for (int i = 0; i < cumulativeTestPoints.length; i++) {
+ TestUtils.assertEquals("Incorrect probability density value
returned for "
+ + cumulativeTestPoints[i], logDensityTestValues[i],
+ ((AbstractRealDistribution)
distribution).logDensity(cumulativeTestPoints[i]),
+ getTolerance());
+ }
+ }
+
//------------------------ Default test cases -----------------------------
/**
@@ -238,6 +268,15 @@ public abstract class RealDistributionAb
}
/**
+ * Verifies that logarithmic density calculations return expected values
+ * for default test instance data
+ */
+ @Test
+ public void testLogDensities() {
+ verifyLogDensities();
+ }
+
+ /**
* Verifies that probability computations are consistent
*/
@Test
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/ZipfDistributionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/ZipfDistributionTest.java?rev=1533974&r1=1533973&r2=1533974&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/ZipfDistributionTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/distribution/ZipfDistributionTest.java
Sun Oct 20 20:42:41 2013
@@ -73,6 +73,17 @@ public class ZipfDistributionTest extend
0.0569028586912, 0.0487738788782, 0.0426771440184,
0.0379352391275, 0.0341417152147, 0};
}
+ /**
+ * Creates the default logarithmic probability density test expected
values.
+ * Reference values are from R, version 2.14.1.
+ */
+ @Override
+ public double[] makeLogDensityTestValues() {
+ return new double[] {Double.NEGATIVE_INFINITY,
Double.NEGATIVE_INFINITY,
+ -1.0746d, -1.7678d, -2.1733d, -2.4609d, -2.6841d, -2.8664d, -3.0206d,
-3.1541d,
+ -3.2719d, -3.3772d, Double.NEGATIVE_INFINITY};
+ }
+
/** Creates the default cumulative probability density test input values */
@Override
public int[] makeCumulativeTestPoints() {
