http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java b/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java deleted file mode 100644 index 411f1a2..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/ExponentialDistribution.java +++ /dev/null @@ -1,351 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.util.CombinatoricsUtils; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.ResizableDoubleArray; - -/** - * Implementation of the exponential distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/Exponential_distribution";>Exponential distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/ExponentialDistribution.html";>Exponential distribution (MathWorld)</a> - */ -public class ExponentialDistribution extends AbstractRealDistribution { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier */ - private static final long serialVersionUID = 2401296428283614780L; - /** - * Used when generating Exponential samples. - * Table containing the constants - * q_i = sum_{j=1}^i (ln 2)^j/j! = ln 2 + (ln 2)^2/2 + ... + (ln 2)^i/i! - * until the largest representable fraction below 1 is exceeded. - * - * Note that - * 1 = 2 - 1 = exp(ln 2) - 1 = sum_{n=1}^infty (ln 2)^n / n! - * thus q_i -> 1 as i -> +inf, - * so the higher i, the closer to one we get (the series is not alternating). - * - * By trying, n = 16 in Java is enough to reach 1.0. - */ - 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; - - /** - * Initialize tables. - */ - static { - /** - * Filling EXPONENTIAL_SA_QI table. - * Note that we don't want qi = 0 in the table. - */ - final double LN2 = FastMath.log(2); - double qi = 0; - int i = 1; - - /** - * ArithmeticUtils provides factorials up to 20, so let's use that - * limit together with Precision.EPSILON to generate the following - * code (a priori, we know that there will be 16 elements, but it is - * better to not hardcode it). - */ - final ResizableDoubleArray ra = new ResizableDoubleArray(20); - - while (qi < 1) { - qi += FastMath.pow(LN2, i) / CombinatoricsUtils.factorial(i); - ra.addElement(qi); - ++i; - } - - EXPONENTIAL_SA_QI = ra.getElements(); - } - - /** - * Create an exponential distribution with the given mean. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param mean mean of this distribution. - */ - public ExponentialDistribution(double mean) { - this(mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Create an exponential distribution with the given mean. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param mean Mean of this distribution. - * @param inverseCumAccuracy Maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code mean <= 0}. - * @since 2.1 - */ - public ExponentialDistribution(double mean, double inverseCumAccuracy) { - this(new Well19937c(), mean, inverseCumAccuracy); - } - - /** - * Creates an exponential distribution. - * - * @param rng Random number generator. - * @param mean Mean of this distribution. - * @throws NotStrictlyPositiveException if {@code mean <= 0}. - * @since 3.3 - */ - public ExponentialDistribution(RandomGenerator rng, double mean) - throws NotStrictlyPositiveException { - this(rng, mean, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Creates an exponential distribution. - * - * @param rng Random number generator. - * @param mean Mean of this distribution. - * @param inverseCumAccuracy Maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code mean <= 0}. - * @since 3.1 - */ - public ExponentialDistribution(RandomGenerator rng, - double mean, - double inverseCumAccuracy) - throws NotStrictlyPositiveException { - super(rng); - - if (mean <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean); - } - this.mean = mean; - logMean = FastMath.log(mean); - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * Access the mean. - * - * @return the mean. - */ - public double getMean() { - return mean; - } - - /** {@inheritDoc} */ - public double density(double x) { - final double logDensity = logDensity(x); - return logDensity == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logDensity); - } - - /** {@inheritDoc} **/ - @Override - public double logDensity(double x) { - if (x < 0) { - return Double.NEGATIVE_INFINITY; - } - return -x / mean - logMean; - } - - /** - * {@inheritDoc} - * - * The implementation of this method is based on: - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/ExponentialDistribution.html";> - * Exponential Distribution</a>, equation (1).</li> - * </ul> - */ - public double cumulativeProbability(double x) { - double ret; - if (x <= 0.0) { - ret = 0.0; - } else { - ret = 1.0 - FastMath.exp(-x / mean); - } - return ret; - } - - /** - * {@inheritDoc} - * - * Returns {@code 0} when {@code p= = 0} and - * {@code Double.POSITIVE_INFINITY} when {@code p == 1}. - */ - @Override - public double inverseCumulativeProbability(double p) throws OutOfRangeException { - double ret; - - if (p < 0.0 || p > 1.0) { - throw new OutOfRangeException(p, 0.0, 1.0); - } else if (p == 1.0) { - ret = Double.POSITIVE_INFINITY; - } else { - ret = -mean * FastMath.log(1.0 - p); - } - - return ret; - } - - /** - * {@inheritDoc} - * - * <p><strong>Algorithm Description</strong>: this implementation uses the - * <a href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html";> - * Inversion Method</a> to generate exponentially distributed random values - * from uniform deviates.</p> - * - * @return a random value. - * @since 2.2 - */ - @Override - public double sample() { - // Step 1: - double a = 0; - double u = random.nextDouble(); - - // Step 2 and 3: - while (u < 0.5) { - a += EXPONENTIAL_SA_QI[0]; - u *= 2; - } - - // Step 4 (now u >= 0.5): - u += u - 1; - - // Step 5: - if (u <= EXPONENTIAL_SA_QI[0]) { - return mean * (a + u); - } - - // Step 6: - int i = 0; // Should be 1, be we iterate before it in while using 0 - double u2 = random.nextDouble(); - double umin = u2; - - // Step 7 and 8: - do { - ++i; - u2 = random.nextDouble(); - - if (u2 < umin) { - umin = u2; - } - - // Step 8: - } while (u > EXPONENTIAL_SA_QI[i]); // Ensured to exit since EXPONENTIAL_SA_QI[MAX] = 1 - - return mean * (a + umin * EXPONENTIAL_SA_QI[0]); - } - - /** {@inheritDoc} */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@inheritDoc} - * - * For mean parameter {@code k}, the mean is {@code k}. - */ - public double getNumericalMean() { - return getMean(); - } - - /** - * {@inheritDoc} - * - * For mean parameter {@code k}, the variance is {@code k^2}. - */ - public double getNumericalVariance() { - final double m = getMean(); - return m * m; - } - - /** - * {@inheritDoc} - * - * The lower bound of the support is always 0 no matter the mean parameter. - * - * @return lower bound of the support (always 0) - */ - public double getSupportLowerBound() { - return 0; - } - - /** - * {@inheritDoc} - * - * The upper bound of the support is always positive infinity - * no matter the mean parameter. - * - * @return upper bound of the support (always Double.POSITIVE_INFINITY) - */ - public double getSupportUpperBound() { - return Double.POSITIVE_INFINITY; - } - - /** {@inheritDoc} */ - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** {@inheritDoc} */ - public boolean isSupportUpperBoundInclusive() { - return false; - } - - /** - * {@inheritDoc} - * - * The support of this distribution is connected. - * - * @return {@code true} - */ - public boolean isSupportConnected() { - return true; - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/FDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/FDistribution.java b/src/main/java/org/apache/commons/math3/distribution/FDistribution.java deleted file mode 100644 index bd98c37..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/FDistribution.java +++ /dev/null @@ -1,328 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.special.Beta; -import org.apache.commons.math3.util.FastMath; - -/** - * Implementation of the F-distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/F-distribution";>F-distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/F-Distribution.html";>F-distribution (MathWorld)</a> - */ -public class FDistribution extends AbstractRealDistribution { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier. */ - private static final long serialVersionUID = -8516354193418641566L; - /** The numerator degrees of freedom. */ - private final double numeratorDegreesOfFreedom; - /** The numerator degrees of freedom. */ - private final double denominatorDegreesOfFreedom; - /** Inverse cumulative probability accuracy. */ - private final double solverAbsoluteAccuracy; - /** Cached numerical variance */ - private double numericalVariance = Double.NaN; - /** Whether or not the numerical variance has been calculated */ - private boolean numericalVarianceIsCalculated = false; - - /** - * Creates an F distribution using the given degrees of freedom. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @throws NotStrictlyPositiveException if - * {@code numeratorDegreesOfFreedom <= 0} or - * {@code denominatorDegreesOfFreedom <= 0}. - */ - public FDistribution(double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom) - throws NotStrictlyPositiveException { - this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, - DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Creates an F distribution using the given degrees of freedom - * and inverse cumulative probability accuracy. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates. - * @throws NotStrictlyPositiveException if - * {@code numeratorDegreesOfFreedom <= 0} or - * {@code denominatorDegreesOfFreedom <= 0}. - * @since 2.1 - */ - public FDistribution(double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom, - double inverseCumAccuracy) - throws NotStrictlyPositiveException { - this(new Well19937c(), numeratorDegreesOfFreedom, - denominatorDegreesOfFreedom, inverseCumAccuracy); - } - - /** - * Creates an F distribution. - * - * @param rng Random number generator. - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or - * {@code denominatorDegreesOfFreedom <= 0}. - * @since 3.3 - */ - public FDistribution(RandomGenerator rng, - double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom) - throws NotStrictlyPositiveException { - this(rng, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Creates an F distribution. - * - * @param rng Random number generator. - * @param numeratorDegreesOfFreedom Numerator degrees of freedom. - * @param denominatorDegreesOfFreedom Denominator degrees of freedom. - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates. - * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or - * {@code denominatorDegreesOfFreedom <= 0}. - * @since 3.1 - */ - public FDistribution(RandomGenerator rng, - double numeratorDegreesOfFreedom, - double denominatorDegreesOfFreedom, - double inverseCumAccuracy) - throws NotStrictlyPositiveException { - super(rng); - - if (numeratorDegreesOfFreedom <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, - numeratorDegreesOfFreedom); - } - if (denominatorDegreesOfFreedom <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM, - denominatorDegreesOfFreedom); - } - this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom; - this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom; - solverAbsoluteAccuracy = inverseCumAccuracy; - } - - /** - * {@inheritDoc} - * - * @since 2.1 - */ - public double density(double x) { - return FastMath.exp(logDensity(x)); - } - - /** {@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} - * - * The implementation of this method is based on - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/F-Distribution.html";> - * F-Distribution</a>, equation (4). - * </li> - * </ul> - */ - public double cumulativeProbability(double x) { - double ret; - if (x <= 0) { - ret = 0; - } else { - double n = numeratorDegreesOfFreedom; - double m = denominatorDegreesOfFreedom; - - ret = Beta.regularizedBeta((n * x) / (m + n * x), - 0.5 * n, - 0.5 * m); - } - return ret; - } - - /** - * Access the numerator degrees of freedom. - * - * @return the numerator degrees of freedom. - */ - public double getNumeratorDegreesOfFreedom() { - return numeratorDegreesOfFreedom; - } - - /** - * Access the denominator degrees of freedom. - * - * @return the denominator degrees of freedom. - */ - public double getDenominatorDegreesOfFreedom() { - return denominatorDegreesOfFreedom; - } - - /** {@inheritDoc} */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@inheritDoc} - * - * For denominator degrees of freedom parameter {@code b}, the mean is - * <ul> - * <li>if {@code b > 2} then {@code b / (b - 2)},</li> - * <li>else undefined ({@code Double.NaN}). - * </ul> - */ - public double getNumericalMean() { - final double denominatorDF = getDenominatorDegreesOfFreedom(); - - if (denominatorDF > 2) { - return denominatorDF / (denominatorDF - 2); - } - - return Double.NaN; - } - - /** - * {@inheritDoc} - * - * For numerator degrees of freedom parameter {@code a} and denominator - * degrees of freedom parameter {@code b}, the variance is - * <ul> - * <li> - * if {@code b > 4} then - * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]}, - * </li> - * <li>else undefined ({@code Double.NaN}). - * </ul> - */ - public double getNumericalVariance() { - if (!numericalVarianceIsCalculated) { - numericalVariance = calculateNumericalVariance(); - numericalVarianceIsCalculated = true; - } - return numericalVariance; - } - - /** - * used by {@link #getNumericalVariance()} - * - * @return the variance of this distribution - */ - protected double calculateNumericalVariance() { - final double denominatorDF = getDenominatorDegreesOfFreedom(); - - if (denominatorDF > 4) { - final double numeratorDF = getNumeratorDegreesOfFreedom(); - final double denomDFMinusTwo = denominatorDF - 2; - - return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) / - ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) ); - } - - return Double.NaN; - } - - /** - * {@inheritDoc} - * - * The lower bound of the support is always 0 no matter the parameters. - * - * @return lower bound of the support (always 0) - */ - public double getSupportLowerBound() { - return 0; - } - - /** - * {@inheritDoc} - * - * The upper bound of the support is always positive infinity - * no matter the parameters. - * - * @return upper bound of the support (always Double.POSITIVE_INFINITY) - */ - public double getSupportUpperBound() { - return Double.POSITIVE_INFINITY; - } - - /** {@inheritDoc} */ - public boolean isSupportLowerBoundInclusive() { - return false; - } - - /** {@inheritDoc} */ - public boolean isSupportUpperBoundInclusive() { - return false; - } - - /** - * {@inheritDoc} - * - * The support of this distribution is connected. - * - * @return {@code true} - */ - public boolean isSupportConnected() { - return true; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java b/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java deleted file mode 100644 index 4f60fa9..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/GammaDistribution.java +++ /dev/null @@ -1,513 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.special.Gamma; -import org.apache.commons.math3.util.FastMath; - -/** - * Implementation of the Gamma distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/Gamma_distribution";>Gamma distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/GammaDistribution.html";>Gamma distribution (MathWorld)</a> - */ -public class GammaDistribution extends AbstractRealDistribution { - /** - * Default inverse cumulative probability accuracy. - * @since 2.1 - */ - public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; - /** Serializable version identifier. */ - private static final long serialVersionUID = 20120524L; - /** The shape parameter. */ - private final double shape; - /** The scale parameter. */ - private final double scale; - /** - * The constant value of {@code shape + g + 0.5}, where {@code g} is the - * Lanczos constant {@link Gamma#LANCZOS_G}. - */ - private final double shiftedShape; - /** - * The constant value of - * {@code 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 #density(double)}, when no overflow occurs with the natural - * 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)}, - * where {@code L(shape)} is the Lanczos approximation returned by - * {@link Gamma#lanczos(double)}. This prefactor is used in - * {@link #density(double)}, when overflow occurs with the natural - * 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 - * calculation overflows. - */ - private final double minY; - /** - * Upper bound on {@code log(y)} ({@code y = x / scale}) for the selection - * of the computation method in {@link #density(double)}. For - * {@code log(y) >= maxLogY}, the natural calculation overflows. - */ - private final double maxLogY; - /** Inverse cumulative probability accuracy. */ - private final double solverAbsoluteAccuracy; - - /** - * Creates a new gamma distribution with specified values of the shape and - * scale parameters. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param shape the shape parameter - * @param scale the scale parameter - * @throws NotStrictlyPositiveException if {@code shape <= 0} or - * {@code scale <= 0}. - */ - public GammaDistribution(double shape, double scale) throws NotStrictlyPositiveException { - this(shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Creates a new gamma distribution with specified values of the shape and - * scale parameters. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param shape the shape parameter - * @param scale the scale parameter - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code shape <= 0} or - * {@code scale <= 0}. - * @since 2.1 - */ - public GammaDistribution(double shape, double scale, double inverseCumAccuracy) - throws NotStrictlyPositiveException { - this(new Well19937c(), shape, scale, inverseCumAccuracy); - } - - /** - * Creates a Gamma distribution. - * - * @param rng Random number generator. - * @param shape the shape parameter - * @param scale the scale parameter - * @throws NotStrictlyPositiveException if {@code shape <= 0} or - * {@code scale <= 0}. - * @since 3.3 - */ - public GammaDistribution(RandomGenerator rng, double shape, double scale) - throws NotStrictlyPositiveException { - this(rng, shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); - } - - /** - * Creates a Gamma distribution. - * - * @param rng Random number generator. - * @param shape the shape parameter - * @param scale the scale parameter - * @param inverseCumAccuracy the maximum absolute error in inverse - * cumulative probability estimates (defaults to - * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). - * @throws NotStrictlyPositiveException if {@code shape <= 0} or - * {@code scale <= 0}. - * @since 3.1 - */ - public GammaDistribution(RandomGenerator rng, - double shape, - double scale, - double inverseCumAccuracy) - throws NotStrictlyPositiveException { - super(rng); - - if (shape <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape); - } - if (scale <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale); - } - - this.shape = shape; - this.scale = scale; - this.solverAbsoluteAccuracy = inverseCumAccuracy; - 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); - } - - /** - * Returns the shape parameter of {@code this} distribution. - * - * @return the shape parameter - * @deprecated as of version 3.1, {@link #getShape()} should be preferred. - * This method will be removed in version 4.0. - */ - @Deprecated - public double getAlpha() { - return shape; - } - - /** - * Returns the shape parameter of {@code this} distribution. - * - * @return the shape parameter - * @since 3.1 - */ - public double getShape() { - return shape; - } - - /** - * Returns the scale parameter of {@code this} distribution. - * - * @return the scale parameter - * @deprecated as of version 3.1, {@link #getScale()} should be preferred. - * This method will be removed in version 4.0. - */ - @Deprecated - public double getBeta() { - return scale; - } - - /** - * Returns the scale parameter of {@code this} distribution. - * - * @return the scale parameter - * @since 3.1 - */ - public double getScale() { - return scale; - } - - /** {@inheritDoc} */ - public double density(double x) { - /* The present method must return the value of - * - * 1 x a - x - * ---------- (-) exp(---) - * x Gamma(a) b b - * - * where a is the shape parameter, and b the scale parameter. - * Substituting the Lanczos approximation of Gamma(a) leads to the - * following expression of the density - * - * a e 1 y a - * - sqrt(------------------) ---- (-----------) exp(a - y + g), - * x 2 pi (a + g + 0.5) L(a) a + g + 0.5 - * - * where y = x / b. The above formula is the "natural" computation, which - * is implemented when no overflow is likely to occur. If overflow occurs - * with the natural computation, the following identity is used. It is - * based on the BOOST library - * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html - * Formula (15) needs adaptations, which are detailed below. - * - * y a - * (-----------) exp(a - y + g) - * a + g + 0.5 - * y - a - g - 0.5 y (g + 0.5) - * = exp(a log1pm(---------------) - ----------- + g), - * a + g + 0.5 a + g + 0.5 - * - * where log1pm(z) = log(1 + z) - z. Therefore, the value to be - * returned is - * - * a e 1 - * - sqrt(------------------) ---- - * x 2 pi (a + g + 0.5) L(a) - * y - a - g - 0.5 y (g + 0.5) - * * exp(a log1pm(---------------) - ----------- + g). - * a + g + 0.5 a + g + 0.5 - */ - if (x < 0) { - return 0; - } - 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 densityPrefactor2 / x * FastMath.exp(aux3); - } - /* - * Natural calculation. - */ - return densityPrefactor1 * FastMath.exp(-y) * 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} - * - * The implementation of this method is based on: - * <ul> - * <li> - * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html";> - * Chi-Squared Distribution</a>, equation (9). - * </li> - * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. - * Belmont, CA: Duxbury Press. - * </li> - * </ul> - */ - public double cumulativeProbability(double x) { - double ret; - - if (x <= 0) { - ret = 0; - } else { - ret = Gamma.regularizedGammaP(shape, x / scale); - } - - return ret; - } - - /** {@inheritDoc} */ - @Override - protected double getSolverAbsoluteAccuracy() { - return solverAbsoluteAccuracy; - } - - /** - * {@inheritDoc} - * - * For shape parameter {@code alpha} and scale parameter {@code beta}, the - * mean is {@code alpha * beta}. - */ - public double getNumericalMean() { - return shape * scale; - } - - /** - * {@inheritDoc} - * - * For shape parameter {@code alpha} and scale parameter {@code beta}, the - * variance is {@code alpha * beta^2}. - * - * @return {@inheritDoc} - */ - public double getNumericalVariance() { - return shape * scale * scale; - } - - /** - * {@inheritDoc} - * - * The lower bound of the support is always 0 no matter the parameters. - * - * @return lower bound of the support (always 0) - */ - public double getSupportLowerBound() { - return 0; - } - - /** - * {@inheritDoc} - * - * The upper bound of the support is always positive infinity - * no matter the parameters. - * - * @return upper bound of the support (always Double.POSITIVE_INFINITY) - */ - public double getSupportUpperBound() { - return Double.POSITIVE_INFINITY; - } - - /** {@inheritDoc} */ - public boolean isSupportLowerBoundInclusive() { - return true; - } - - /** {@inheritDoc} */ - public boolean isSupportUpperBoundInclusive() { - return false; - } - - /** - * {@inheritDoc} - * - * The support of this distribution is connected. - * - * @return {@code true} - */ - public boolean isSupportConnected() { - return true; - } - - /** - * <p>This implementation uses the following algorithms: </p> - * - * <p>For 0 < shape < 1: <br/> - * Ahrens, J. H. and Dieter, U., <i>Computer methods for - * sampling from gamma, beta, Poisson and binomial distributions.</i> - * Computing, 12, 223-246, 1974.</p> - * - * <p>For shape >= 1: <br/> - * Marsaglia and Tsang, <i>A Simple Method for Generating - * Gamma Variables.</i> ACM Transactions on Mathematical Software, - * Volume 26 Issue 3, September, 2000.</p> - * - * @return random value sampled from the Gamma(shape, scale) distribution - */ - @Override - public double sample() { - if (shape < 1) { - // [1]: p. 228, Algorithm GS - - while (true) { - // Step 1: - final double u = random.nextDouble(); - final double bGS = 1 + shape / FastMath.E; - final double p = bGS * u; - - if (p <= 1) { - // Step 2: - - final double x = FastMath.pow(p, 1 / shape); - final double u2 = random.nextDouble(); - - if (u2 > FastMath.exp(-x)) { - // Reject - continue; - } else { - return scale * x; - } - } else { - // Step 3: - - final double x = -1 * FastMath.log((bGS - p) / shape); - final double u2 = random.nextDouble(); - - if (u2 > FastMath.pow(x, shape - 1)) { - // Reject - continue; - } else { - return scale * x; - } - } - } - } - - // Now shape >= 1 - - final double d = shape - 0.333333333333333333; - final double c = 1 / (3 * FastMath.sqrt(d)); - - while (true) { - final double x = random.nextGaussian(); - final double v = (1 + c * x) * (1 + c * x) * (1 + c * x); - - if (v <= 0) { - continue; - } - - final double x2 = x * x; - final double u = random.nextDouble(); - - // Squeeze - if (u < 1 - 0.0331 * x2 * x2) { - return scale * d * v; - } - - if (FastMath.log(u) < 0.5 * x2 + d * (1 - v + FastMath.log(v))) { - return scale * d * v; - } - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java b/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java deleted file mode 100644 index f82a3ec..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/GeometricDistribution.java +++ /dev/null @@ -1,173 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.util.FastMath; - -/** - * Implementation of the geometric distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/Geometric_distribution";>Geometric distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/GeometricDistribution.html";>Geometric Distribution (MathWorld)</a> - * @since 3.3 - */ -public class GeometricDistribution extends AbstractIntegerDistribution { - - /** Serializable version identifier. */ - private static final long serialVersionUID = 20130507L; - /** The probability of success. */ - private final double probabilityOfSuccess; - - /** - * Create a geometric distribution with the given probability of success. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param p probability of success. - * @throws OutOfRangeException if {@code p <= 0} or {@code p > 1}. - */ - public GeometricDistribution(double p) { - this(new Well19937c(), p); - } - - /** - * Creates a geometric distribution. - * - * @param rng Random number generator. - * @param p Probability of success. - * @throws OutOfRangeException if {@code p <= 0} or {@code p > 1}. - */ - public GeometricDistribution(RandomGenerator rng, double p) { - super(rng); - - if (p <= 0 || p > 1) { - throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT, p, 0, 1); - } - - probabilityOfSuccess = p; - } - - /** - * Access the probability of success for this distribution. - * - * @return the probability of success. - */ - public double getProbabilityOfSuccess() { - return probabilityOfSuccess; - } - - /** {@inheritDoc} */ - public double probability(int x) { - double ret; - if (x < 0) { - ret = 0.0; - } else { - final double p = probabilityOfSuccess; - ret = FastMath.pow(1 - p, x) * p; - } - 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; - if (x < 0) { - ret = 0.0; - } else { - final double p = probabilityOfSuccess; - ret = 1.0 - FastMath.pow(1 - p, x + 1); - } - return ret; - } - - /** - * {@inheritDoc} - * - * For probability parameter {@code p}, the mean is {@code (1 - p) / p}. - */ - public double getNumericalMean() { - final double p = probabilityOfSuccess; - return (1 - p) / p; - } - - /** - * {@inheritDoc} - * - * For probability parameter {@code p}, the variance is - * {@code (1 - p) / (p * p)}. - */ - public double getNumericalVariance() { - final double p = probabilityOfSuccess; - return (1 - p) / (p * p); - } - - /** - * {@inheritDoc} - * - * The lower bound of the support is always 0. - * - * @return lower bound of the support (always 0) - */ - public int getSupportLowerBound() { - return 0; - } - - /** - * {@inheritDoc} - * - * The upper bound of the support is infinite (which we approximate as - * {@code Integer.MAX_VALUE}). - * - * @return upper bound of the support (always Integer.MAX_VALUE) - */ - public int getSupportUpperBound() { - return Integer.MAX_VALUE; - } - - /** - * {@inheritDoc} - * - * The support of this distribution is connected. - * - * @return {@code true} - */ - public boolean isSupportConnected() { - return true; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/GumbelDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/GumbelDistribution.java b/src/main/java/org/apache/commons/math3/distribution/GumbelDistribution.java deleted file mode 100644 index 85dbedd..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/GumbelDistribution.java +++ /dev/null @@ -1,166 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathUtils; - -/** - * This class implements the Gumbel distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/Gumbel_distribution";>Gumbel Distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/GumbelDistribution.html";>Gumbel Distribution (Mathworld)</a> - * - * @since 3.4 - */ -public class GumbelDistribution extends AbstractRealDistribution { - - /** Serializable version identifier. */ - private static final long serialVersionUID = 20141003; - - /** - * Approximation of Euler's constant - * see http://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html - */ - private static final double EULER = FastMath.PI / (2 * FastMath.E); - - /** The location parameter. */ - private final double mu; - /** The scale parameter. */ - private final double beta; - - /** - * Build a new instance. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param mu location parameter - * @param beta scale parameter (must be positive) - * @throws NotStrictlyPositiveException if {@code beta <= 0} - */ - public GumbelDistribution(double mu, double beta) { - this(new Well19937c(), mu, beta); - } - - /** - * Build a new instance. - * - * @param rng Random number generator - * @param mu location parameter - * @param beta scale parameter (must be positive) - * @throws NotStrictlyPositiveException if {@code beta <= 0} - */ - public GumbelDistribution(RandomGenerator rng, double mu, double beta) { - super(rng); - - if (beta <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, beta); - } - - this.beta = beta; - this.mu = mu; - } - - /** - * Access the location parameter, {@code mu}. - * - * @return the location parameter. - */ - public double getLocation() { - return mu; - } - - /** - * Access the scale parameter, {@code beta}. - * - * @return the scale parameter. - */ - public double getScale() { - return beta; - } - - /** {@inheritDoc} */ - public double density(double x) { - final double z = (x - mu) / beta; - final double t = FastMath.exp(-z); - return FastMath.exp(-z - t) / beta; - } - - /** {@inheritDoc} */ - public double cumulativeProbability(double x) { - final double z = (x - mu) / beta; - return FastMath.exp(-FastMath.exp(-z)); - } - - @Override - public double inverseCumulativeProbability(double p) throws OutOfRangeException { - if (p < 0.0 || p > 1.0) { - throw new OutOfRangeException(p, 0.0, 1.0); - } else if (p == 0) { - return Double.NEGATIVE_INFINITY; - } else if (p == 1) { - return Double.POSITIVE_INFINITY; - } - return mu - FastMath.log(-FastMath.log(p)) * beta; - } - - /** {@inheritDoc} */ - public double getNumericalMean() { - return mu + EULER * beta; - } - - /** {@inheritDoc} */ - public double getNumericalVariance() { - return (MathUtils.PI_SQUARED) / 6.0 * (beta * beta); - } - - /** {@inheritDoc} */ - public double getSupportLowerBound() { - return Double.NEGATIVE_INFINITY; - } - - /** {@inheritDoc} */ - public double getSupportUpperBound() { - return Double.POSITIVE_INFINITY; - } - - /** {@inheritDoc} */ - public boolean isSupportLowerBoundInclusive() { - return false; - } - - /** {@inheritDoc} */ - public boolean isSupportUpperBoundInclusive() { - return false; - } - - /** {@inheritDoc} */ - public boolean isSupportConnected() { - return true; - } - -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java b/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java deleted file mode 100644 index 7a1436a..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/HypergeometricDistribution.java +++ /dev/null @@ -1,347 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NotPositiveException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.random.RandomGenerator; -import org.apache.commons.math3.random.Well19937c; -import org.apache.commons.math3.util.FastMath; - -/** - * Implementation of the hypergeometric distribution. - * - * @see <a href="http://en.wikipedia.org/wiki/Hypergeometric_distribution";>Hypergeometric distribution (Wikipedia)</a> - * @see <a href="http://mathworld.wolfram.com/HypergeometricDistribution.html";>Hypergeometric distribution (MathWorld)</a> - */ -public class HypergeometricDistribution extends AbstractIntegerDistribution { - /** Serializable version identifier. */ - private static final long serialVersionUID = -436928820673516179L; - /** The number of successes in the population. */ - private final int numberOfSuccesses; - /** The population size. */ - private final int populationSize; - /** The sample size. */ - private final int sampleSize; - /** Cached numerical variance */ - private double numericalVariance = Double.NaN; - /** Whether or not the numerical variance has been calculated */ - private boolean numericalVarianceIsCalculated = false; - - /** - * Construct a new hypergeometric distribution with the specified population - * size, number of successes in the population, and sample size. - * <p> - * <b>Note:</b> this constructor will implicitly create an instance of - * {@link Well19937c} as random generator to be used for sampling only (see - * {@link #sample()} and {@link #sample(int)}). In case no sampling is - * needed for the created distribution, it is advised to pass {@code null} - * as random generator via the appropriate constructors to avoid the - * additional initialisation overhead. - * - * @param populationSize Population size. - * @param numberOfSuccesses Number of successes in the population. - * @param sampleSize Sample size. - * @throws NotPositiveException if {@code numberOfSuccesses < 0}. - * @throws NotStrictlyPositiveException if {@code populationSize <= 0}. - * @throws NumberIsTooLargeException if {@code numberOfSuccesses > populationSize}, - * or {@code sampleSize > populationSize}. - */ - public HypergeometricDistribution(int populationSize, int numberOfSuccesses, int sampleSize) - throws NotPositiveException, NotStrictlyPositiveException, NumberIsTooLargeException { - this(new Well19937c(), populationSize, numberOfSuccesses, sampleSize); - } - - /** - * Creates a new hypergeometric distribution. - * - * @param rng Random number generator. - * @param populationSize Population size. - * @param numberOfSuccesses Number of successes in the population. - * @param sampleSize Sample size. - * @throws NotPositiveException if {@code numberOfSuccesses < 0}. - * @throws NotStrictlyPositiveException if {@code populationSize <= 0}. - * @throws NumberIsTooLargeException if {@code numberOfSuccesses > populationSize}, - * or {@code sampleSize > populationSize}. - * @since 3.1 - */ - public HypergeometricDistribution(RandomGenerator rng, - int populationSize, - int numberOfSuccesses, - int sampleSize) - throws NotPositiveException, NotStrictlyPositiveException, NumberIsTooLargeException { - super(rng); - - if (populationSize <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.POPULATION_SIZE, - populationSize); - } - if (numberOfSuccesses < 0) { - throw new NotPositiveException(LocalizedFormats.NUMBER_OF_SUCCESSES, - numberOfSuccesses); - } - if (sampleSize < 0) { - throw new NotPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, - sampleSize); - } - - if (numberOfSuccesses > populationSize) { - throw new NumberIsTooLargeException(LocalizedFormats.NUMBER_OF_SUCCESS_LARGER_THAN_POPULATION_SIZE, - numberOfSuccesses, populationSize, true); - } - if (sampleSize > populationSize) { - throw new NumberIsTooLargeException(LocalizedFormats.SAMPLE_SIZE_LARGER_THAN_POPULATION_SIZE, - sampleSize, populationSize, true); - } - - this.numberOfSuccesses = numberOfSuccesses; - this.populationSize = populationSize; - this.sampleSize = sampleSize; - } - - /** {@inheritDoc} */ - public double cumulativeProbability(int x) { - double ret; - - int[] domain = getDomain(populationSize, numberOfSuccesses, sampleSize); - if (x < domain[0]) { - ret = 0.0; - } else if (x >= domain[1]) { - ret = 1.0; - } else { - ret = innerCumulativeProbability(domain[0], x, 1); - } - - return ret; - } - - /** - * Return the domain for the given hypergeometric distribution parameters. - * - * @param n Population size. - * @param m Number of successes in the population. - * @param k Sample size. - * @return a two element array containing the lower and upper bounds of the - * hypergeometric distribution. - */ - private int[] getDomain(int n, int m, int k) { - return new int[] { getLowerDomain(n, m, k), getUpperDomain(m, k) }; - } - - /** - * Return the lowest domain value for the given hypergeometric distribution - * parameters. - * - * @param n Population size. - * @param m Number of successes in the population. - * @param k Sample size. - * @return the lowest domain value of the hypergeometric distribution. - */ - private int getLowerDomain(int n, int m, int k) { - return FastMath.max(0, m - (n - k)); - } - - /** - * Access the number of successes. - * - * @return the number of successes. - */ - public int getNumberOfSuccesses() { - return numberOfSuccesses; - } - - /** - * Access the population size. - * - * @return the population size. - */ - public int getPopulationSize() { - return populationSize; - } - - /** - * Access the sample size. - * - * @return the sample size. - */ - public int getSampleSize() { - return sampleSize; - } - - /** - * Return the highest domain value for the given hypergeometric distribution - * parameters. - * - * @param m Number of successes in the population. - * @param k Sample size. - * @return the highest domain value of the hypergeometric distribution. - */ - private int getUpperDomain(int m, int k) { - return FastMath.min(k, m); - } - - /** {@inheritDoc} */ - public double probability(int x) { - final double logProbability = logProbability(x); - return logProbability == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logProbability); - } - - /** {@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)}. - * - * @param x Value at which the CDF is evaluated. - * @return the upper tail CDF for this distribution. - * @since 1.1 - */ - public double upperCumulativeProbability(int x) { - double ret; - - final int[] domain = getDomain(populationSize, numberOfSuccesses, sampleSize); - if (x <= domain[0]) { - ret = 1.0; - } else if (x > domain[1]) { - ret = 0.0; - } else { - ret = innerCumulativeProbability(domain[1], x, -1); - } - - return ret; - } - - /** - * For this distribution, {@code X}, this method returns - * {@code P(x0 <= X <= x1)}. - * This probability is computed by summing the point probabilities for the - * values {@code x0, x0 + 1, x0 + 2, ..., x1}, in the order directed by - * {@code dx}. - * - * @param x0 Inclusive lower bound. - * @param x1 Inclusive upper bound. - * @param dx Direction of summation (1 indicates summing from x0 to x1, and - * 0 indicates summing from x1 to x0). - * @return {@code P(x0 <= X <= x1)}. - */ - private double innerCumulativeProbability(int x0, int x1, int dx) { - double ret = probability(x0); - while (x0 != x1) { - x0 += dx; - ret += probability(x0); - } - return ret; - } - - /** - * {@inheritDoc} - * - * For population size {@code N}, number of successes {@code m}, and sample - * size {@code n}, the mean is {@code n * m / N}. - */ - public double getNumericalMean() { - return getSampleSize() * (getNumberOfSuccesses() / (double) getPopulationSize()); - } - - /** - * {@inheritDoc} - * - * For population size {@code N}, number of successes {@code m}, and sample - * size {@code n}, the variance is - * {@code [n * m * (N - n) * (N - m)] / [N^2 * (N - 1)]}. - */ - public double getNumericalVariance() { - if (!numericalVarianceIsCalculated) { - numericalVariance = calculateNumericalVariance(); - numericalVarianceIsCalculated = true; - } - return numericalVariance; - } - - /** - * Used by {@link #getNumericalVariance()}. - * - * @return the variance of this distribution - */ - protected double calculateNumericalVariance() { - final double N = getPopulationSize(); - final double m = getNumberOfSuccesses(); - final double n = getSampleSize(); - return (n * m * (N - n) * (N - m)) / (N * N * (N - 1)); - } - - /** - * {@inheritDoc} - * - * For population size {@code N}, number of successes {@code m}, and sample - * size {@code n}, the lower bound of the support is - * {@code max(0, n + m - N)}. - * - * @return lower bound of the support - */ - public int getSupportLowerBound() { - return FastMath.max(0, - getSampleSize() + getNumberOfSuccesses() - getPopulationSize()); - } - - /** - * {@inheritDoc} - * - * For number of successes {@code m} and sample size {@code n}, the upper - * bound of the support is {@code min(m, n)}. - * - * @return upper bound of the support - */ - public int getSupportUpperBound() { - return FastMath.min(getNumberOfSuccesses(), getSampleSize()); - } - - /** - * {@inheritDoc} - * - * The support of this distribution is connected. - * - * @return {@code true} - */ - public boolean isSupportConnected() { - return true; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/IntegerDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/IntegerDistribution.java b/src/main/java/org/apache/commons/math3/distribution/IntegerDistribution.java deleted file mode 100644 index 9ab4a04..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/IntegerDistribution.java +++ /dev/null @@ -1,155 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.distribution; - -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.OutOfRangeException; - -/** - * Interface for distributions on the integers. - * - */ -public interface IntegerDistribution { - /** - * For a random variable {@code X} whose values are distributed according - * to this distribution, this method returns {@code P(X = x)}. In other - * words, this method represents the probability mass function (PMF) - * for the distribution. - * - * @param x the point at which the PMF is evaluated - * @return the value of the probability mass function at {@code x} - */ - double probability(int x); - - /** - * For a random variable {@code X} whose values are distributed according - * to this distribution, this method returns {@code P(X <= x)}. In other - * words, this method represents the (cumulative) distribution function - * (CDF) for this distribution. - * - * @param x the point at which the CDF is evaluated - * @return the probability that a random variable with this - * distribution takes a value less than or equal to {@code x} - */ - double cumulativeProbability(int x); - - /** - * For a random variable {@code X} whose values are distributed according - * to this distribution, this method returns {@code P(x0 < X <= x1)}. - * - * @param x0 the exclusive lower bound - * @param x1 the inclusive upper bound - * @return the probability that a random variable with this distribution - * will take a value between {@code x0} and {@code x1}, - * excluding the lower and including the upper endpoint - * @throws NumberIsTooLargeException if {@code x0 > x1} - */ - double cumulativeProbability(int x0, int x1) throws NumberIsTooLargeException; - - /** - * Computes the quantile function of this distribution. - * For a random variable {@code X} distributed according to this distribution, - * the returned value is - * <ul> - * <li><code>inf{x in Z | P(X<=x) >= p}</code> for {@code 0 < p <= 1},</li> - * <li><code>inf{x in Z | P(X<=x) > 0}</code> for {@code p = 0}.</li> - * </ul> - * If the result exceeds the range of the data type {@code int}, - * then {@code Integer.MIN_VALUE} or {@code Integer.MAX_VALUE} is returned. - * - * @param p the cumulative probability - * @return the smallest {@code p}-quantile of this distribution - * (largest 0-quantile for {@code p = 0}) - * @throws OutOfRangeException if {@code p < 0} or {@code p > 1} - */ - int inverseCumulativeProbability(double p) throws OutOfRangeException; - - /** - * Use this method to get the numerical value of the mean of this - * distribution. - * - * @return the mean or {@code Double.NaN} if it is not defined - */ - double getNumericalMean(); - - /** - * Use this method to get the numerical value of the variance of this - * distribution. - * - * @return the variance (possibly {@code Double.POSITIVE_INFINITY} or - * {@code Double.NaN} if it is not defined) - */ - double getNumericalVariance(); - - /** - * Access the lower bound of the support. This method must return the same - * value as {@code inverseCumulativeProbability(0)}. In other words, this - * method must return - * <p><code>inf {x in Z | P(X <= x) > 0}</code>.</p> - * - * @return lower bound of the support ({@code Integer.MIN_VALUE} - * for negative infinity) - */ - int getSupportLowerBound(); - - /** - * Access the upper bound of the support. This method must return the same - * value as {@code inverseCumulativeProbability(1)}. In other words, this - * method must return - * <p><code>inf {x in R | P(X <= x) = 1}</code>.</p> - * - * @return upper bound of the support ({@code Integer.MAX_VALUE} - * for positive infinity) - */ - int getSupportUpperBound(); - - /** - * Use this method to get information about whether the support is - * connected, i.e. whether all integers between the lower and upper bound of - * the support are included in the support. - * - * @return whether the support is connected or not - */ - boolean isSupportConnected(); - - /** - * Reseed the random generator used to generate samples. - * - * @param seed the new seed - * @since 3.0 - */ - void reseedRandomGenerator(long seed); - - /** - * Generate a random value sampled from this distribution. - * - * @return a random value - * @since 3.0 - */ - int sample(); - - /** - * Generate a random sample from the distribution. - * - * @param sampleSize the number of random values to generate - * @return an array representing the random sample - * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException - * if {@code sampleSize} is not positive - * @since 3.0 - */ - int[] sample(int sampleSize); -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java b/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java deleted file mode 100644 index 7af514d..0000000 --- a/src/main/java/org/apache/commons/math3/distribution/KolmogorovSmirnovDistribution.java +++ /dev/null @@ -1,357 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.distribution; - -import java.io.Serializable; -import java.math.BigDecimal; - -import org.apache.commons.math3.exception.MathArithmeticException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.fraction.BigFraction; -import org.apache.commons.math3.fraction.BigFractionField; -import org.apache.commons.math3.fraction.FractionConversionException; -import org.apache.commons.math3.linear.Array2DRowFieldMatrix; -import org.apache.commons.math3.linear.Array2DRowRealMatrix; -import org.apache.commons.math3.linear.FieldMatrix; -import org.apache.commons.math3.linear.RealMatrix; -import org.apache.commons.math3.util.FastMath; - -/** - * Implementation of the Kolmogorov-Smirnov distribution. - * - * <p> - * Treats the distribution of the two-sided {@code P(D_n < d)} where - * {@code D_n = sup_x |G(x) - G_n (x)|} for the theoretical cdf {@code G} and - * the empirical cdf {@code G_n}. - * </p> - * <p> - * This implementation is based on [1] with certain quick decisions for extreme - * values given in [2]. - * </p> - * <p> - * In short, when wanting to evaluate {@code P(D_n < d)}, the method in [1] is - * to write {@code d = (k - h) / n} for positive integer {@code k} and - * {@code 0 <= h < 1}. Then {@code P(D_n < d) = (n! / n^n) * t_kk}, where - * {@code t_kk} is the {@code (k, k)}'th entry in the special matrix - * {@code H^n}, i.e. {@code H} to the {@code n}'th power. - * </p> - * <p> - * References: - * <ul> - * <li>[1] <a href="http://www.jstatsoft.org/v08/i18/";> - * Evaluating Kolmogorov's Distribution</a> by George Marsaglia, Wai - * Wan Tsang, and Jingbo Wang</li> - * <li>[2] <a href="http://www.jstatsoft.org/v39/i11/";> - * Computing the Two-Sided Kolmogorov-Smirnov Distribution</a> by Richard Simard - * and Pierre L'Ecuyer</li> - * </ul> - * Note that [1] contains an error in computing h, refer to - * <a href="https://issues.apache.org/jira/browse/MATH-437";>MATH-437</a> for details. - * </p> - * - * @see <a href="http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test";> - * Kolmogorov-Smirnov test (Wikipedia)</a> - * @deprecated to be removed in version 4.0 - - * use {@link org.apache.commons.math3.stat.inference.KolmogorovSmirnovTest} - */ -public class KolmogorovSmirnovDistribution implements Serializable { - - /** Serializable version identifier. */ - private static final long serialVersionUID = -4670676796862967187L; - - /** Number of observations. */ - private int n; - - /** - * @param n Number of observations - * @throws NotStrictlyPositiveException if {@code n <= 0} - */ - public KolmogorovSmirnovDistribution(int n) - throws NotStrictlyPositiveException { - if (n <= 0) { - throw new NotStrictlyPositiveException(LocalizedFormats.NOT_POSITIVE_NUMBER_OF_SAMPLES, n); - } - - this.n = n; - } - - /** - * Calculates {@code P(D_n < d)} using method described in [1] with quick - * decisions for extreme values given in [2] (see above). The result is not - * exact as with - * {@link KolmogorovSmirnovDistribution#cdfExact(double)} because - * calculations are based on {@code double} rather than - * {@link org.apache.commons.math3.fraction.BigFraction}. - * - * @param d statistic - * @return the two-sided probability of {@code P(D_n < d)} - * @throws MathArithmeticException if algorithm fails to convert {@code h} - * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing - * {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - public double cdf(double d) throws MathArithmeticException { - return this.cdf(d, false); - } - - /** - * Calculates {@code P(D_n < d)} using method described in [1] with quick - * decisions for extreme values given in [2] (see above). The result is - * exact in the sense that BigFraction/BigReal is used everywhere at the - * expense of very slow execution time. Almost never choose this in real - * applications unless you are very sure; this is almost solely for - * verification purposes. Normally, you would choose - * {@link KolmogorovSmirnovDistribution#cdf(double)} - * - * @param d statistic - * @return the two-sided probability of {@code P(D_n < d)} - * @throws MathArithmeticException if algorithm fails to convert {@code h} - * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing - * {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - public double cdfExact(double d) throws MathArithmeticException { - return this.cdf(d, true); - } - - /** - * Calculates {@code P(D_n < d)} using method described in [1] with quick - * decisions for extreme values given in [2] (see above). - * - * @param d statistic - * @param exact whether the probability should be calculated exact using - * {@link org.apache.commons.math3.fraction.BigFraction} everywhere at the - * expense of very slow execution time, or if {@code double} should be used - * convenient places to gain speed. Almost never choose {@code true} in real - * applications unless you are very sure; {@code true} is almost solely for - * verification purposes. - * @return the two-sided probability of {@code P(D_n < d)} - * @throws MathArithmeticException if algorithm fails to convert {@code h} - * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing - * {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - public double cdf(double d, boolean exact) throws MathArithmeticException { - - final double ninv = 1 / ((double) n); - final double ninvhalf = 0.5 * ninv; - - if (d <= ninvhalf) { - - return 0; - - } else if (ninvhalf < d && d <= ninv) { - - double res = 1; - double f = 2 * d - ninv; - - // n! f^n = n*f * (n-1)*f * ... * 1*x - for (int i = 1; i <= n; ++i) { - res *= i * f; - } - - return res; - - } else if (1 - ninv <= d && d < 1) { - - return 1 - 2 * FastMath.pow(1 - d, n); - - } else if (1 <= d) { - - return 1; - } - - return exact ? exactK(d) : roundedK(d); - } - - /** - * Calculates the exact value of {@code P(D_n < d)} using method described - * in [1] and {@link org.apache.commons.math3.fraction.BigFraction} (see - * above). - * - * @param d statistic - * @return the two-sided probability of {@code P(D_n < d)} - * @throws MathArithmeticException if algorithm fails to convert {@code h} - * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing - * {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - private double exactK(double d) throws MathArithmeticException { - - final int k = (int) FastMath.ceil(n * d); - - final FieldMatrix<BigFraction> H = this.createH(d); - final FieldMatrix<BigFraction> Hpower = H.power(n); - - BigFraction pFrac = Hpower.getEntry(k - 1, k - 1); - - for (int i = 1; i <= n; ++i) { - pFrac = pFrac.multiply(i).divide(n); - } - - /* - * BigFraction.doubleValue converts numerator to double and the - * denominator to double and divides afterwards. That gives NaN quite - * easy. This does not (scale is the number of digits): - */ - return pFrac.bigDecimalValue(20, BigDecimal.ROUND_HALF_UP).doubleValue(); - } - - /** - * Calculates {@code P(D_n < d)} using method described in [1] and doubles - * (see above). - * - * @param d statistic - * @return the two-sided probability of {@code P(D_n < d)} - * @throws MathArithmeticException if algorithm fails to convert {@code h} - * to a {@link org.apache.commons.math3.fraction.BigFraction} in expressing - * {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - private double roundedK(double d) throws MathArithmeticException { - - final int k = (int) FastMath.ceil(n * d); - final FieldMatrix<BigFraction> HBigFraction = this.createH(d); - final int m = HBigFraction.getRowDimension(); - - /* - * Here the rounding part comes into play: use - * RealMatrix instead of FieldMatrix<BigFraction> - */ - final RealMatrix H = new Array2DRowRealMatrix(m, m); - - for (int i = 0; i < m; ++i) { - for (int j = 0; j < m; ++j) { - H.setEntry(i, j, HBigFraction.getEntry(i, j).doubleValue()); - } - } - - final RealMatrix Hpower = H.power(n); - - double pFrac = Hpower.getEntry(k - 1, k - 1); - - for (int i = 1; i <= n; ++i) { - pFrac *= (double) i / (double) n; - } - - return pFrac; - } - - /*** - * Creates {@code H} of size {@code m x m} as described in [1] (see above). - * - * @param d statistic - * @return H matrix - * @throws NumberIsTooLargeException if fractional part is greater than 1 - * @throws FractionConversionException if algorithm fails to convert - * {@code h} to a {@link org.apache.commons.math3.fraction.BigFraction} in - * expressing {@code d} as {@code (k - h) / m} for integer {@code k, m} and - * {@code 0 <= h < 1}. - */ - private FieldMatrix<BigFraction> createH(double d) - throws NumberIsTooLargeException, FractionConversionException { - - int k = (int) FastMath.ceil(n * d); - - int m = 2 * k - 1; - double hDouble = k - n * d; - - if (hDouble >= 1) { - throw new NumberIsTooLargeException(hDouble, 1.0, false); - } - - BigFraction h = null; - - try { - h = new BigFraction(hDouble, 1.0e-20, 10000); - } catch (FractionConversionException e1) { - try { - h = new BigFraction(hDouble, 1.0e-10, 10000); - } catch (FractionConversionException e2) { - h = new BigFraction(hDouble, 1.0e-5, 10000); - } - } - - final BigFraction[][] Hdata = new BigFraction[m][m]; - - /* - * Start by filling everything with either 0 or 1. - */ - for (int i = 0; i < m; ++i) { - for (int j = 0; j < m; ++j) { - if (i - j + 1 < 0) { - Hdata[i][j] = BigFraction.ZERO; - } else { - Hdata[i][j] = BigFraction.ONE; - } - } - } - - /* - * Setting up power-array to avoid calculating the same value twice: - * hPowers[0] = h^1 ... hPowers[m-1] = h^m - */ - final BigFraction[] hPowers = new BigFraction[m]; - hPowers[0] = h; - for (int i = 1; i < m; ++i) { - hPowers[i] = h.multiply(hPowers[i - 1]); - } - - /* - * First column and last row has special values (each other reversed). - */ - for (int i = 0; i < m; ++i) { - Hdata[i][0] = Hdata[i][0].subtract(hPowers[i]); - Hdata[m - 1][i] = Hdata[m - 1][i].subtract(hPowers[m - i - 1]); - } - - /* - * [1] states: "For 1/2 < h < 1 the bottom left element of the matrix - * should be (1 - 2*h^m + (2h - 1)^m )/m!" Since 0 <= h < 1, then if h > - * 1/2 is sufficient to check: - */ - if (h.compareTo(BigFraction.ONE_HALF) == 1) { - Hdata[m - 1][0] = Hdata[m - 1][0].add(h.multiply(2).subtract(1).pow(m)); - } - - /* - * Aside from the first column and last row, the (i, j)-th element is - * 1/(i - j + 1)! if i - j + 1 >= 0, else 0. 1's and 0's are already - * put, so only division with (i - j + 1)! is needed in the elements - * that have 1's. There is no need to calculate (i - j + 1)! and then - * divide - small steps avoid overflows. - * - * Note that i - j + 1 > 0 <=> i + 1 > j instead of j'ing all the way to - * m. Also note that it is started at g = 2 because dividing by 1 isn't - * really necessary. - */ - for (int i = 0; i < m; ++i) { - for (int j = 0; j < i + 1; ++j) { - if (i - j + 1 > 0) { - for (int g = 2; g <= i - j + 1; ++g) { - Hdata[i][j] = Hdata[i][j].divide(g); - } - } - } - } - - return new Array2DRowFieldMatrix<BigFraction>(BigFractionField.getInstance(), Hdata); - } -}
