mirror of
https://github.com/apache/commons-math.git
synced 2025-03-04 15:39:17 +00:00
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 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1533974 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Phil Steitz
parent
97bd8fe5de
commit
422fc2a426
3
pom.xml
3
pom.xml
@ -180,6 +180,9 @@
|
||||
<contributor>
|
||||
<name>Larry Diamond</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Aleksei Dievskii</name>
|
||||
</contributor>
|
||||
<contributor>
|
||||
<name>Rodrigo di Lorenzo Lopes</name>
|
||||
</contributor>
|
||||
|
||||
@ -51,6 +51,10 @@ If the output is not quite correct, check for invisible trailing spaces!
|
||||
</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.
|
||||
|
||||
@ -232,4 +232,23 @@ public abstract class AbstractIntegerDistribution implements IntegerDistribution
|
||||
}
|
||||
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));
|
||||
}
|
||||
}
|
||||
|
||||
@ -286,5 +286,23 @@ implements RealDistribution, Serializable {
|
||||
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));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -156,6 +156,29 @@ public class BetaDistribution extends AbstractRealDistribution {
|
||||
}
|
||||
}
|
||||
|
||||
/** {@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) {
|
||||
|
||||
@ -110,6 +110,20 @@ public class BinomialDistribution extends AbstractIntegerDistribution {
|
||||
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;
|
||||
|
||||
@ -108,6 +108,12 @@ public class ChiSquaredDistribution extends AbstractRealDistribution {
|
||||
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);
|
||||
|
||||
@ -56,6 +56,8 @@ public class ExponentialDistribution extends AbstractRealDistribution {
|
||||
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 extends AbstractRealDistribution {
|
||||
throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
|
||||
}
|
||||
this.mean = mean;
|
||||
logMean = FastMath.log(mean);
|
||||
solverAbsoluteAccuracy = inverseCumAccuracy;
|
||||
}
|
||||
|
||||
@ -164,6 +167,12 @@ public class ExponentialDistribution extends AbstractRealDistribution {
|
||||
return FastMath.exp(-x / mean) / mean;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} **/
|
||||
@Override
|
||||
public double logDensity(double x) {
|
||||
return -x / mean - logMean;
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
|
||||
@ -154,6 +154,21 @@ public class FDistribution extends AbstractRealDistribution {
|
||||
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}
|
||||
*
|
||||
|
||||
@ -56,6 +56,15 @@ public class GammaDistribution extends AbstractRealDistribution {
|
||||
* calculation.
|
||||
*/
|
||||
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)},
|
||||
@ -65,6 +74,15 @@ public class GammaDistribution extends AbstractRealDistribution {
|
||||
* calculation.
|
||||
*/
|
||||
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
|
||||
@ -159,9 +177,14 @@ public class GammaDistribution extends AbstractRealDistribution {
|
||||
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 AbstractRealDistribution {
|
||||
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}
|
||||
*
|
||||
|
||||
@ -85,6 +85,19 @@ public class GeometricDistribution extends AbstractIntegerDistribution {
|
||||
return ret;
|
||||
}
|
||||
|
||||
/** {@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;
|
||||
|
||||
@ -214,6 +214,30 @@ public class HypergeometricDistribution extends AbstractIntegerDistribution {
|
||||
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)}.
|
||||
*
|
||||
|
||||
@ -79,6 +79,21 @@ public class LevyDistribution extends AbstractRealDistribution {
|
||||
return FastMath.sqrt(f / FastMath.PI) * FastMath.exp(-f) /delta;
|
||||
}
|
||||
|
||||
/** {@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
|
||||
|
||||
@ -71,6 +71,8 @@ public class LogNormalDistribution extends AbstractRealDistribution {
|
||||
|
||||
/** 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 extends AbstractRealDistribution {
|
||||
|
||||
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 extends AbstractRealDistribution {
|
||||
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}
|
||||
*
|
||||
|
||||
@ -49,6 +49,8 @@ public class NormalDistribution extends AbstractRealDistribution {
|
||||
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 AbstractRealDistribution {
|
||||
|
||||
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 AbstractRealDistribution {
|
||||
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}
|
||||
*
|
||||
|
||||
@ -174,6 +174,18 @@ public class ParetoDistribution extends AbstractRealDistribution {
|
||||
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>
|
||||
|
||||
@ -68,6 +68,12 @@ public class PascalDistribution extends AbstractIntegerDistribution {
|
||||
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 AbstractIntegerDistribution {
|
||||
|
||||
numberOfSuccesses = r;
|
||||
probabilityOfSuccess = p;
|
||||
logProbabilityOfSuccess = FastMath.log(p);
|
||||
log1mProbabilityOfSuccess = FastMath.log1p(-p);
|
||||
}
|
||||
|
||||
/**
|
||||
@ -146,6 +154,21 @@ public class PascalDistribution extends AbstractIntegerDistribution {
|
||||
return ret;
|
||||
}
|
||||
|
||||
/** {@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;
|
||||
|
||||
@ -174,6 +174,22 @@ public class PoissonDistribution extends AbstractIntegerDistribution {
|
||||
return ret;
|
||||
}
|
||||
|
||||
/** {@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) {
|
||||
|
||||
@ -129,6 +129,18 @@ public class TDistribution extends AbstractRealDistribution {
|
||||
nPlus1Over2 * FastMath.log(1 + x * x / n));
|
||||
}
|
||||
|
||||
/** {@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;
|
||||
|
||||
@ -174,6 +174,26 @@ public class WeibullDistribution extends AbstractRealDistribution {
|
||||
return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape);
|
||||
}
|
||||
|
||||
/** {@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;
|
||||
|
||||
@ -114,6 +114,16 @@ public class ZipfDistribution extends AbstractIntegerDistribution {
|
||||
return (1.0 / FastMath.pow(x, exponent)) / generalizedHarmonic(numberOfElements, exponent);
|
||||
}
|
||||
|
||||
/** {@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) {
|
||||
|
||||
@ -18,6 +18,7 @@ package org.apache.commons.math3.distribution;
|
||||
|
||||
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 IntegerDistributionAbstractTest {
|
||||
/** 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 IntegerDistributionAbstractTest {
|
||||
/** 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 IntegerDistributionAbstractTest {
|
||||
distribution = makeDistribution();
|
||||
densityTestPoints = makeDensityTestPoints();
|
||||
densityTestValues = makeDensityTestValues();
|
||||
logDensityTestValues = makeLogDensityTestValues();
|
||||
cumulativeTestPoints = makeCumulativeTestPoints();
|
||||
cumulativeTestValues = makeCumulativeTestValues();
|
||||
inverseCumulativeTestPoints = makeInverseCumulativeTestPoints();
|
||||
@ -119,6 +140,7 @@ public abstract class IntegerDistributionAbstractTest {
|
||||
distribution = null;
|
||||
densityTestPoints = null;
|
||||
densityTestValues = null;
|
||||
logDensityTestValues = null;
|
||||
cumulativeTestPoints = null;
|
||||
cumulativeTestValues = null;
|
||||
inverseCumulativeTestPoints = null;
|
||||
@ -139,6 +161,19 @@ public abstract class IntegerDistributionAbstractTest {
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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
|
||||
@ -175,6 +210,15 @@ public abstract class IntegerDistributionAbstractTest {
|
||||
verifyDensities();
|
||||
}
|
||||
|
||||
/**
|
||||
* 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
|
||||
|
||||
@ -68,4 +68,15 @@ public class LevyDistributionTest extends RealDistributionAbstractTest {
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* 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};
|
||||
}
|
||||
}
|
||||
|
||||
@ -67,6 +67,18 @@ public class PoissonDistributionTest extends IntegerDistributionAbstractTest {
|
||||
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.
|
||||
*/
|
||||
|
||||
@ -95,6 +95,9 @@ public abstract class RealDistributionAbstractTest {
|
||||
/** 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 RealDistributionAbstractTest {
|
||||
/** 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 RealDistributionAbstractTest {
|
||||
inverseCumulativeTestPoints = makeInverseCumulativeTestPoints();
|
||||
inverseCumulativeTestValues = makeInverseCumulativeTestValues();
|
||||
densityTestValues = makeDensityTestValues();
|
||||
logDensityTestValues = makeLogDensityTestValues();
|
||||
}
|
||||
|
||||
/**
|
||||
@ -147,6 +163,7 @@ public abstract class RealDistributionAbstractTest {
|
||||
inverseCumulativeTestPoints = null;
|
||||
inverseCumulativeTestValues = null;
|
||||
densityTestValues = null;
|
||||
logDensityTestValues = null;
|
||||
}
|
||||
|
||||
//-------------------- Verification methods -------------------------------
|
||||
@ -208,6 +225,19 @@ public abstract class RealDistributionAbstractTest {
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* 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 -----------------------------
|
||||
|
||||
/**
|
||||
@ -237,6 +267,15 @@ public abstract class RealDistributionAbstractTest {
|
||||
verifyDensities();
|
||||
}
|
||||
|
||||
/**
|
||||
* Verifies that logarithmic density calculations return expected values
|
||||
* for default test instance data
|
||||
*/
|
||||
@Test
|
||||
public void testLogDensities() {
|
||||
verifyLogDensities();
|
||||
}
|
||||
|
||||
/**
|
||||
* Verifies that probability computations are consistent
|
||||
*/
|
||||
|
||||
@ -73,6 +73,17 @@ public class ZipfDistributionTest extends IntegerDistributionAbstractTest {
|
||||
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() {
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user