You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/16 23:40:03 UTC
[33/82] [partial] [math] Update for next development iteration:
commons-math4
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);
- }
-}