You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2018/01/22 01:15:33 UTC
[4/6] commons-statistics git commit: Fixed file paths.
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
deleted file mode 100644
index 225b8f1..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/PoissonDistribution.java
+++ /dev/null
@@ -1,238 +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.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.RegularizedGamma;
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-import org.apache.commons.rng.sampling.distribution.PoissonSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Poisson_distribution">Poisson distribution</a>.
- */
-public class PoissonDistribution extends AbstractDiscreteDistribution {
- /** ln(2 π). */
- private static final double LOG_TWO_PI = Math.log(2 * Math.PI);
- /** Default maximum number of iterations. */
- private static final int DEFAULT_MAX_ITERATIONS = 10000000;
- /** Default convergence criterion. */
- private static final double DEFAULT_EPSILON = 1e-12;
- /** Distribution used to compute normal approximation. */
- private final NormalDistribution normal;
- /** Mean of the distribution. */
- private final double mean;
- /** Maximum number of iterations for cumulative probability. */
- private final int maxIterations;
- /** Convergence criterion for cumulative probability. */
- private final double epsilon;
-
- /**
- * Creates a new Poisson distribution with specified mean.
- *
- * @param p the Poisson mean
- * @throws IllegalArgumentException if {@code p <= 0}.
- */
- public PoissonDistribution(double p) {
- this(p, DEFAULT_EPSILON, DEFAULT_MAX_ITERATIONS);
- }
-
- /**
- * Creates a new Poisson distribution with specified mean, convergence
- * criterion and maximum number of iterations.
- *
- * @param p Poisson mean.
- * @param epsilon Convergence criterion for cumulative probabilities.
- * @param maxIterations the maximum number of iterations for cumulative
- * probabilities.
- * @throws IllegalArgumentException if {@code p <= 0}.
- */
- public PoissonDistribution(double p,
- double epsilon,
- int maxIterations) {
- if (p <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE, p);
- }
- mean = p;
- this.epsilon = epsilon;
- this.maxIterations = maxIterations;
-
- normal = new NormalDistribution(p, Math.sqrt(p));
- }
-
- /**
- * Creates a new Poisson distribution with the specified mean and
- * convergence criterion.
- *
- * @param p Poisson mean.
- * @param epsilon Convergence criterion for cumulative probabilities.
- * @throws IllegalArgumentException if {@code p <= 0}.
- */
- public PoissonDistribution(double p,
- double epsilon) {
- this(p, epsilon, DEFAULT_MAX_ITERATIONS);
- }
-
- /**
- * Creates a new Poisson distribution with the specified mean and maximum
- * number of iterations.
- *
- * @param p Poisson mean.
- * @param maxIterations Maximum number of iterations for cumulative
- * probabilities.
- */
- public PoissonDistribution(double p,
- int maxIterations) {
- this(p, DEFAULT_EPSILON, maxIterations);
- }
-
- /**
- * Get the mean for the distribution.
- *
- * @return the mean for the distribution.
- */
- public double getMean() {
- return mean;
- }
-
- /** {@inheritDoc} */
- @Override
- public double probability(int x) {
- final double logProbability = logProbability(x);
- return logProbability == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logProbability);
- }
-
- /** {@inheritDoc} */
- @Override
- public double logProbability(int x) {
- double ret;
- if (x < 0 || x == Integer.MAX_VALUE) {
- ret = Double.NEGATIVE_INFINITY;
- } else if (x == 0) {
- ret = -mean;
- } else {
- ret = -SaddlePointExpansion.getStirlingError(x) -
- SaddlePointExpansion.getDeviancePart(x, mean) -
- 0.5 * LOG_TWO_PI - 0.5 * Math.log(x);
- }
- return ret;
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(int x) {
- if (x < 0) {
- return 0;
- }
- if (x == Integer.MAX_VALUE) {
- return 1;
- }
- return RegularizedGamma.Q.value((double) x + 1, mean, epsilon,
- maxIterations);
- }
-
- /**
- * Calculates the Poisson distribution function using a normal
- * approximation. The {@code N(mean, sqrt(mean))} distribution is used
- * to approximate the Poisson distribution. The computation uses
- * "half-correction" (evaluating the normal distribution function at
- * {@code x + 0.5}).
- *
- * @param x Upper bound, inclusive.
- * @return the distribution function value calculated using a normal
- * approximation.
- */
- public double normalApproximateProbability(int x) {
- // Calculate the probability using half-correction.
- return normal.cumulativeProbability(x + 0.5);
- }
-
- /**
- * {@inheritDoc}
- *
- * For mean parameter {@code p}, the mean is {@code p}.
- */
- @Override
- public double getNumericalMean() {
- return getMean();
- }
-
- /**
- * {@inheritDoc}
- *
- * For mean parameter {@code p}, the variance is {@code p}.
- */
- @Override
- public double getNumericalVariance() {
- return getMean();
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is always 0 no matter the mean parameter.
- *
- * @return lower bound of the support (always 0)
- */
- @Override
- public int getSupportLowerBound() {
- return 0;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is positive infinity,
- * regardless of the parameter values. There is no integer infinity,
- * so this method returns {@code Integer.MAX_VALUE}.
- *
- * @return upper bound of the support (always {@code Integer.MAX_VALUE} for
- * positive infinity)
- */
- @Override
- public int getSupportUpperBound() {
- return Integer.MAX_VALUE;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-
- /**{@inheritDoc} */
- @Override
- public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
- return new DiscreteDistribution.Sampler() {
- /**
- * Poisson distribution sampler.
- */
- private final DiscreteSampler sampler = new PoissonSampler(rng, mean);
-
- /**{@inheritDoc} */
- @Override
- public int sample() {
- return sampler.sample();
- }
- };
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
deleted file mode 100644
index 7bb847a..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/SaddlePointExpansion.java
+++ /dev/null
@@ -1,191 +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.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Utility class used by various distributions to accurately compute their
- * respective probability mass functions. The implementation for this class is
- * based on the Catherine Loader's
- * <a href="http://www.herine.net/stat/software/dbinom.html">dbinom</a> routines.
- *
- * This class is not intended to be called directly.
- *
- * @since 1.0
- */
-final class SaddlePointExpansion {
- /** 2 π */
- private static final double TWO_PI = 2 * Math.PI;
- /** 1/2 * log(2 π). */
- private static final double HALF_LOG_TWO_PI = 0.5 * Math.log(TWO_PI);
-
- /** exact Stirling expansion error for certain values. */
- private static final double[] EXACT_STIRLING_ERRORS = { 0.0, /* 0.0 */
- 0.1534264097200273452913848, /* 0.5 */
- 0.0810614667953272582196702, /* 1.0 */
- 0.0548141210519176538961390, /* 1.5 */
- 0.0413406959554092940938221, /* 2.0 */
- 0.03316287351993628748511048, /* 2.5 */
- 0.02767792568499833914878929, /* 3.0 */
- 0.02374616365629749597132920, /* 3.5 */
- 0.02079067210376509311152277, /* 4.0 */
- 0.01848845053267318523077934, /* 4.5 */
- 0.01664469118982119216319487, /* 5.0 */
- 0.01513497322191737887351255, /* 5.5 */
- 0.01387612882307074799874573, /* 6.0 */
- 0.01281046524292022692424986, /* 6.5 */
- 0.01189670994589177009505572, /* 7.0 */
- 0.01110455975820691732662991, /* 7.5 */
- 0.010411265261972096497478567, /* 8.0 */
- 0.009799416126158803298389475, /* 8.5 */
- 0.009255462182712732917728637, /* 9.0 */
- 0.008768700134139385462952823, /* 9.5 */
- 0.008330563433362871256469318, /* 10.0 */
- 0.007934114564314020547248100, /* 10.5 */
- 0.007573675487951840794972024, /* 11.0 */
- 0.007244554301320383179543912, /* 11.5 */
- 0.006942840107209529865664152, /* 12.0 */
- 0.006665247032707682442354394, /* 12.5 */
- 0.006408994188004207068439631, /* 13.0 */
- 0.006171712263039457647532867, /* 13.5 */
- 0.005951370112758847735624416, /* 14.0 */
- 0.005746216513010115682023589, /* 14.5 */
- 0.005554733551962801371038690 /* 15.0 */
- };
-
- /**
- * Forbid construction.
- */
- private SaddlePointExpansion() {}
-
- /**
- * Compute the error of Stirling's series at the given value.
- * <p>
- * References:
- * <ol>
- * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web
- * Resource. <a target="_blank"
- * href="http://mathworld.wolfram.com/StirlingsSeries.html">
- * http://mathworld.wolfram.com/StirlingsSeries.html</a></li>
- * </ol>
- * </p>
- *
- * @param z the value.
- * @return the Striling's series error.
- */
- static double getStirlingError(double z) {
- double ret;
- if (z < 15.0) {
- double z2 = 2.0 * z;
- if (Math.floor(z2) == z2) {
- ret = EXACT_STIRLING_ERRORS[(int) z2];
- } else {
- ret = LogGamma.value(z + 1.0) - (z + 0.5) * Math.log(z) +
- z - HALF_LOG_TWO_PI;
- }
- } else {
- double z2 = z * z;
- ret = (0.083333333333333333333 -
- (0.00277777777777777777778 -
- (0.00079365079365079365079365 -
- (0.000595238095238095238095238 -
- 0.0008417508417508417508417508 /
- z2) / z2) / z2) / z2) / z;
- }
- return ret;
- }
-
- /**
- * A part of the deviance portion of the saddle point approximation.
- * <p>
- * References:
- * <ol>
- * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
- * Probabilities.". <a target="_blank"
- * href="http://www.herine.net/stat/papers/dbinom.pdf">
- * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
- * </ol>
- * </p>
- *
- * @param x the x value.
- * @param mu the average.
- * @return a part of the deviance.
- */
- static double getDeviancePart(double x, double mu) {
- double ret;
- if (Math.abs(x - mu) < 0.1 * (x + mu)) {
- double d = x - mu;
- double v = d / (x + mu);
- double s1 = v * d;
- double s = Double.NaN;
- double ej = 2.0 * x * v;
- v *= v;
- int j = 1;
- while (s1 != s) {
- s = s1;
- ej *= v;
- s1 = s + ej / ((j * 2) + 1);
- ++j;
- }
- ret = s1;
- } else {
- if (x == 0) {
- return mu;
- }
- ret = x * Math.log(x / mu) + mu - x;
- }
- return ret;
- }
-
- /**
- * Compute the logarithm of the PMF for a binomial distribution
- * using the saddle point expansion.
- *
- * @param x the value at which the probability is evaluated.
- * @param n the number of trials.
- * @param p the probability of success.
- * @param q the probability of failure (1 - p).
- * @return log(p(x)).
- */
- static double logBinomialProbability(int x, int n, double p, double q) {
- double ret;
- if (x == 0) {
- if (p < 0.1) {
- ret = -getDeviancePart(n, n * q) - n * p;
- } else {
- if (n == 0) {
- return 0;
- }
- ret = n * Math.log(q);
- }
- } else if (x == n) {
- if (q < 0.1) {
- ret = -getDeviancePart(n, n * p) - n * q;
- } else {
- ret = n * Math.log(p);
- }
- } else {
- ret = getStirlingError(n) - getStirlingError(x) -
- getStirlingError(n - x) - getDeviancePart(x, n * p) -
- getDeviancePart(n - x, n * q);
- final double f = (TWO_PI * x * (n - x)) / n;
- ret = -0.5 * Math.log(f) + ret;
- }
- return ret;
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
deleted file mode 100644
index ef34c1f..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TDistribution.java
+++ /dev/null
@@ -1,180 +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.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.RegularizedBeta;
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Implementation of <a href='http://en.wikipedia.org/wiki/Student's_t-distribution'>Student's t-distribution</a>.
- */
-public class TDistribution extends AbstractContinuousDistribution {
- /** The degrees of freedom. */
- private final double degreesOfFreedom;
- /** degreesOfFreedom / 2 */
- private final double dofOver2;
- /** Cached value. */
- private final double factor;
-
- /**
- * Creates a distribution.
- *
- * @param degreesOfFreedom Degrees of freedom.
- * @throws IllegalArgumentException if {@code degreesOfFreedom <= 0}
- */
- public TDistribution(double degreesOfFreedom) {
- if (degreesOfFreedom <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE,
- degreesOfFreedom);
- }
- this.degreesOfFreedom = degreesOfFreedom;
-
- dofOver2 = 0.5 * degreesOfFreedom;
- factor = LogGamma.value(dofOver2 + 0.5) -
- 0.5 * (Math.log(Math.PI) + Math.log(degreesOfFreedom)) -
- LogGamma.value(dofOver2);
- }
-
- /**
- * Access the degrees of freedom.
- *
- * @return the degrees of freedom.
- */
- public double getDegreesOfFreedom() {
- return degreesOfFreedom;
- }
-
- /** {@inheritDoc} */
- @Override
- public double density(double x) {
- return Math.exp(logDensity(x));
- }
-
- /** {@inheritDoc} */
- @Override
- public double logDensity(double x) {
- final double nPlus1Over2 = dofOver2 + 0.5;
- return factor - nPlus1Over2 * Math.log(1 + x * x / degreesOfFreedom);
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(double x) {
- double ret;
- if (x == 0) {
- ret = 0.5;
- } else {
- final double t =
- RegularizedBeta.value(degreesOfFreedom / (degreesOfFreedom + (x * x)),
- dofOver2,
- 0.5);
- if (x < 0) {
- ret = 0.5 * t;
- } else {
- ret = 1 - 0.5 * t;
- }
- }
-
- return ret;
- }
-
- /**
- * {@inheritDoc}
- *
- * For degrees of freedom parameter {@code df}, the mean is
- * <ul>
- * <li>if {@code df > 1} then {@code 0},</li>
- * <li>else undefined ({@code Double.NaN}).</li>
- * </ul>
- */
- @Override
- public double getNumericalMean() {
- final double df = getDegreesOfFreedom();
-
- if (df > 1) {
- return 0;
- }
-
- return Double.NaN;
- }
-
- /**
- * {@inheritDoc}
- *
- * For degrees of freedom parameter {@code df}, the variance is
- * <ul>
- * <li>if {@code df > 2} then {@code df / (df - 2)},</li>
- * <li>if {@code 1 < df <= 2} then positive infinity
- * ({@code Double.POSITIVE_INFINITY}),</li>
- * <li>else undefined ({@code Double.NaN}).</li>
- * </ul>
- */
- @Override
- public double getNumericalVariance() {
- final double df = getDegreesOfFreedom();
-
- if (df > 2) {
- return df / (df - 2);
- }
-
- if (df > 1 && df <= 2) {
- return Double.POSITIVE_INFINITY;
- }
-
- return Double.NaN;
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is always negative infinity no matter the
- * parameters.
- *
- * @return lower bound of the support (always
- * {@code Double.NEGATIVE_INFINITY})
- */
- @Override
- public double getSupportLowerBound() {
- return Double.NEGATIVE_INFINITY;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is always positive infinity no matter the
- * parameters.
- *
- * @return upper bound of the support (always
- * {@code Double.POSITIVE_INFINITY})
- */
- @Override
- public double getSupportUpperBound() {
- return Double.POSITIVE_INFINITY;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
deleted file mode 100644
index 6a94b62..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/TriangularDistribution.java
+++ /dev/null
@@ -1,222 +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.statistics.distribution;
-
-/**
- * Implementation of the triangular real distribution.
- *
- * @see <a href="http://en.wikipedia.org/wiki/Triangular_distribution">
- * Triangular distribution (Wikipedia)</a>
- *
- * @since 3.0
- */
-public class TriangularDistribution extends AbstractContinuousDistribution {
- /** Serializable version identifier. */
- private static final long serialVersionUID = 20160311L;
- /** Lower limit of this distribution (inclusive). */
- private final double a;
- /** Upper limit of this distribution (inclusive). */
- private final double b;
- /** Mode of this distribution. */
- private final double c;
-
- /**
- * Creates a distribution.
- *
- * @param a Lower limit of this distribution (inclusive).
- * @param b Upper limit of this distribution (inclusive).
- * @param c Mode of this distribution.
- * @throws IllegalArgumentException if {@code a >= b}, if {@code c > b}
- * or if {@code c < a}.
- */
- public TriangularDistribution(double a,
- double c,
- double b) {
- if (a >= b) {
- throw new DistributionException(DistributionException.TOO_LARGE,
- a, b);
- }
- if (c < a) {
- throw new DistributionException(DistributionException.TOO_SMALL,
- c, a);
- }
- if (c > b) {
- throw new DistributionException(DistributionException.TOO_LARGE,
- c, b);
- }
-
- this.a = a;
- this.c = c;
- this.b = b;
- }
-
- /**
- * Gets the mode.
- *
- * @return the mode of the distribution.
- */
- public double getMode() {
- return c;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
- * PDF is given by
- * <ul>
- * <li>{@code 2 * (x - a) / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
- * <li>{@code 2 / (b - a)} if {@code x = c},</li>
- * <li>{@code 2 * (b - x) / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
- * <li>{@code 0} otherwise.
- * </ul>
- */
- @Override
- public double density(double x) {
- if (x < a) {
- return 0;
- }
- if (a <= x && x < c) {
- double divident = 2 * (x - a);
- double divisor = (b - a) * (c - a);
- return divident / divisor;
- }
- if (x == c) {
- return 2 / (b - a);
- }
- if (c < x && x <= b) {
- double divident = 2 * (b - x);
- double divisor = (b - a) * (b - c);
- return divident / divisor;
- }
- return 0;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower limit {@code a}, upper limit {@code b} and mode {@code c}, the
- * CDF is given by
- * <ul>
- * <li>{@code 0} if {@code x < a},</li>
- * <li>{@code (x - a)^2 / [(b - a) * (c - a)]} if {@code a <= x < c},</li>
- * <li>{@code (c - a) / (b - a)} if {@code x = c},</li>
- * <li>{@code 1 - (b - x)^2 / [(b - a) * (b - c)]} if {@code c < x <= b},</li>
- * <li>{@code 1} if {@code x > b}.</li>
- * </ul>
- */
- @Override
- public double cumulativeProbability(double x) {
- if (x < a) {
- return 0;
- }
- if (a <= x && x < c) {
- double divident = (x - a) * (x - a);
- double divisor = (b - a) * (c - a);
- return divident / divisor;
- }
- if (x == c) {
- return (c - a) / (b - a);
- }
- if (c < x && x <= b) {
- double divident = (b - x) * (b - x);
- double divisor = (b - a) * (b - c);
- return 1 - (divident / divisor);
- }
- return 1;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
- * the mean is {@code (a + b + c) / 3}.
- */
- @Override
- public double getNumericalMean() {
- return (a + b + c) / 3;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower limit {@code a}, upper limit {@code b}, and mode {@code c},
- * the variance is {@code (a^2 + b^2 + c^2 - a * b - a * c - b * c) / 18}.
- */
- @Override
- public double getNumericalVariance() {
- return (a * a + b * b + c * c - a * b - a * c - b * c) / 18;
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is equal to the lower limit parameter
- * {@code a} of the distribution.
- *
- * @return lower bound of the support
- */
- @Override
- public double getSupportLowerBound() {
- return a;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is equal to the upper limit parameter
- * {@code b} of the distribution.
- *
- * @return upper bound of the support
- */
- @Override
- public double getSupportUpperBound() {
- return b;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-
- /** {@inheritDoc} */
- @Override
- public double inverseCumulativeProbability(double p) {
- if (p < 0 ||
- p > 1) {
- throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
- }
- if (p == 0) {
- return a;
- }
- if (p == 1) {
- return b;
- }
- if (p < (c - a) / (b - a)) {
- return a + Math.sqrt(p * (b - a) * (c - a));
- }
- return b - Math.sqrt((1 - p) * (b - a) * (b - c));
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
deleted file mode 100644
index 2bb9a0b..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformContinuousDistribution.java
+++ /dev/null
@@ -1,168 +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.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
-import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Uniform_distribution_(continuous)">uniform distribution</a>.
- */
-public class UniformContinuousDistribution extends AbstractContinuousDistribution {
- private final double lower;
- /** Upper bound of this distribution (exclusive). */
- private final double upper;
-
- /**
- * Create a standard uniform real distribution with lower bound (inclusive)
- * equal to zero and upper bound (exclusive) equal to one.
- */
- public UniformContinuousDistribution() {
- this(0, 1);
- }
-
- /**
- * Creates a uniform distribution.
- *
- * @param lower Lower bound of this distribution (inclusive).
- * @param upper Upper bound of this distribution (exclusive).
- * @throws IllegalArgumentException if {@code lower >= upper}.
- */
- public UniformContinuousDistribution(double lower,
- double upper) {
- if (lower >= upper) {
- throw new DistributionException(DistributionException.TOO_LARGE,
- lower, upper);
- }
-
- this.lower = lower;
- this.upper = upper;
- }
-
- /** {@inheritDoc} */
- @Override
- public double density(double x) {
- if (x < lower ||
- x > upper) {
- return 0;
- }
- return 1 / (upper - lower);
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(double x) {
- if (x <= lower) {
- return 0;
- }
- if (x >= upper) {
- return 1;
- }
- return (x - lower) / (upper - lower);
- }
-
- /** {@inheritDoc} */
- @Override
- public double inverseCumulativeProbability(final double p) {
- if (p < 0 ||
- p > 1) {
- throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
- }
- return p * (upper - lower) + lower;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower bound {@code lower} and upper bound {@code upper}, the mean is
- * {@code 0.5 * (lower + upper)}.
- */
- @Override
- public double getNumericalMean() {
- return 0.5 * (lower + upper);
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower bound {@code lower} and upper bound {@code upper}, the
- * variance is {@code (upper - lower)^2 / 12}.
- */
- @Override
- public double getNumericalVariance() {
- double ul = upper - lower;
- return ul * ul / 12;
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is equal to the lower bound parameter
- * of the distribution.
- *
- * @return lower bound of the support
- */
- @Override
- public double getSupportLowerBound() {
- return lower;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is equal to the upper bound parameter
- * of the distribution.
- *
- * @return upper bound of the support
- */
- @Override
- public double getSupportUpperBound() {
- return upper;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-
- /** {@inheritDoc} */
- @Override
- public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
- return new ContinuousDistribution.Sampler() {
- /**
- * Uniform distribution sampler.
- */
- private final ContinuousSampler sampler =
- new ContinuousUniformSampler(rng, lower, upper);
-
- /**{@inheritDoc} */
- @Override
- public double sample() {
- return sampler.sample();
- }
- };
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
deleted file mode 100644
index 41df2bd..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/UniformDiscreteDistribution.java
+++ /dev/null
@@ -1,159 +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.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-import org.apache.commons.rng.sampling.distribution.DiscreteUniformSampler;
-
-/**
- * Implementation of the <a href="http://en.wikipedia.org/wiki/Uniform_distribution_(discrete)">
- * uniform integer distribution</a>.
- */
-public class UniformDiscreteDistribution extends AbstractDiscreteDistribution {
- /** 1 / 12 **/
- private static final double ONE_TWELFTH = 1 / 12d;
- /** Lower bound (inclusive) of this distribution. */
- private final int lower;
- /** Upper bound (inclusive) of this distribution. */
- private final int upper;
- /** "upper" + "lower" (to avoid overflow). */
- private final double upperPlusLower;
- /** "upper" - "lower" (to avoid overflow). */
- private final double upperMinusLower;
-
- /**
- * Creates a new uniform integer distribution using the given lower and
- * upper bounds (both inclusive).
- *
- * @param lower Lower bound (inclusive) of this distribution.
- * @param upper Upper bound (inclusive) of this distribution.
- * @throws IllegalArgumentException if {@code lower > upper}.
- */
- public UniformDiscreteDistribution(int lower,
- int upper) {
- if (lower > upper) {
- throw new DistributionException(DistributionException.TOO_LARGE,
- lower, upper);
- }
- this.lower = lower;
- this.upper = upper;
- upperPlusLower = (double) upper + (double) lower;
- upperMinusLower = (double) upper - (double) lower;
- }
-
- /** {@inheritDoc} */
- @Override
- public double probability(int x) {
- if (x < lower || x > upper) {
- return 0;
- }
- return 1 / (upperMinusLower + 1);
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(int x) {
- if (x < lower) {
- return 0;
- }
- if (x > upper) {
- return 1;
- }
- return (x - lower + 1) / (upperMinusLower + 1);
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower bound {@code lower} and upper bound {@code upper}, the mean is
- * {@code 0.5 * (lower + upper)}.
- */
- @Override
- public double getNumericalMean() {
- return 0.5 * upperPlusLower;
- }
-
- /**
- * {@inheritDoc}
- *
- * For lower bound {@code lower} and upper bound {@code upper}, and
- * {@code n = upper - lower + 1}, the variance is {@code (n^2 - 1) / 12}.
- */
- @Override
- public double getNumericalVariance() {
- double n = upperMinusLower + 1;
- return ONE_TWELFTH * (n * n - 1);
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is equal to the lower bound parameter
- * of the distribution.
- *
- * @return lower bound of the support
- */
- @Override
- public int getSupportLowerBound() {
- return lower;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is equal to the upper bound parameter
- * of the distribution.
- *
- * @return upper bound of the support
- */
- @Override
- public int getSupportUpperBound() {
- return upper;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-
- /**{@inheritDoc} */
- @Override
- public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
- return new DiscreteDistribution.Sampler() {
- /**
- * Discrete uniform distribution sampler.
- */
- private final DiscreteSampler sampler =
- new DiscreteUniformSampler(rng, lower, upper);
-
- /**{@inheritDoc} */
- @Override
- public int sample() {
- return sampler.sample();
- }
- };
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
deleted file mode 100644
index 0e396d8..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/WeibullDistribution.java
+++ /dev/null
@@ -1,220 +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.statistics.distribution;
-
-import org.apache.commons.numbers.gamma.LogGamma;
-
-/**
- * Implementation of the Weibull distribution. This implementation uses the
- * two parameter form of the distribution defined by
- * <a href="http://mathworld.wolfram.com/WeibullDistribution.html">
- * Weibull Distribution</a>, equations (1) and (2).
- *
- * @see <a href="http://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
- * @see <a href="http://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
- *
- * @since 1.1
- */
-public class WeibullDistribution extends AbstractContinuousDistribution {
- /** The shape parameter. */
- private final double shape;
- /** The scale parameter. */
- private final double scale;
-
- /**
- * Creates a distribution.
- *
- * @param alpha Shape parameter.
- * @param beta Scale parameter.
- * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}.
- */
- public WeibullDistribution(double alpha,
- double beta) {
- if (alpha <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE,
- alpha);
- }
- if (beta <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE,
- beta);
- }
- scale = beta;
- shape = alpha;
- }
-
- /**
- * Access the shape parameter, {@code alpha}.
- *
- * @return the shape parameter, {@code alpha}.
- */
- public double getShape() {
- return shape;
- }
-
- /**
- * Access the scale parameter, {@code beta}.
- *
- * @return the scale parameter, {@code beta}.
- */
- public double getScale() {
- return scale;
- }
-
- /** {@inheritDoc} */
- @Override
- public double density(double x) {
- if (x < 0) {
- return 0;
- }
-
- final double xscale = x / scale;
- final double xscalepow = Math.pow(xscale, shape - 1);
-
- /*
- * Math.pow(x / scale, shape) =
- * Math.pow(xscale, shape) =
- * Math.pow(xscale, shape - 1) * xscale
- */
- final double xscalepowshape = xscalepow * xscale;
-
- return (shape / scale) * xscalepow * Math.exp(-xscalepowshape);
- }
-
- /** {@inheritDoc} */
- @Override
- public double logDensity(double x) {
- if (x < 0) {
- return Double.NEGATIVE_INFINITY;
- }
-
- final double xscale = x / scale;
- final double logxscalepow = Math.log(xscale) * (shape - 1);
-
- /*
- * Math.pow(x / scale, shape) =
- * Math.pow(xscale, shape) =
- * Math.pow(xscale, shape - 1) * xscale
- */
- final double xscalepowshape = Math.exp(logxscalepow) * xscale;
-
- return Math.log(shape / scale) + logxscalepow - xscalepowshape;
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(double x) {
- double ret;
- if (x <= 0.0) {
- ret = 0.0;
- } else {
- ret = 1.0 - Math.exp(-Math.pow(x / scale, shape));
- }
- 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) {
- double ret;
- if (p < 0 ||
- p > 1) {
- throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
- } else if (p == 0) {
- ret = 0.0;
- } else if (p == 1) {
- ret = Double.POSITIVE_INFINITY;
- } else {
- ret = scale * Math.pow(-Math.log1p(-p), 1.0 / shape);
- }
- return ret;
- }
-
- /**
- * {@inheritDoc}
- *
- * The mean is {@code scale * Gamma(1 + (1 / shape))}, where {@code Gamma()}
- * is the Gamma-function.
- */
- @Override
- public double getNumericalMean() {
- final double sh = getShape();
- final double sc = getScale();
-
- return sc * Math.exp(LogGamma.value(1 + (1 / sh)));
- }
-
- /**
- * {@inheritDoc}
- *
- * The variance is {@code scale^2 * Gamma(1 + (2 / shape)) - mean^2}
- * where {@code Gamma()} is the Gamma-function.
- */
- @Override
- public double getNumericalVariance() {
- final double sh = getShape();
- final double sc = getScale();
- final double mn = getNumericalMean();
-
- return (sc * sc) * Math.exp(LogGamma.value(1 + (2 / sh))) -
- (mn * mn);
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is always 0 no matter the parameters.
- *
- * @return lower bound of the support (always 0)
- */
- @Override
- 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
- * {@code Double.POSITIVE_INFINITY})
- */
- @Override
- public double getSupportUpperBound() {
- return Double.POSITIVE_INFINITY;
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-}
-
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
deleted file mode 100644
index f0d54a1..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ZipfDistribution.java
+++ /dev/null
@@ -1,236 +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.statistics.distribution;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
-import org.apache.commons.rng.sampling.distribution.RejectionInversionZipfSampler;
-
-/**
- * Implementation of the <a href="https://en.wikipedia.org/wiki/Zipf's_law">Zipf distribution</a>.
- * <p>
- * <strong>Parameters:</strong>
- * For a random variable {@code X} whose values are distributed according to this
- * distribution, the probability mass function is given by
- * <pre>
- * P(X = k) = H(N,s) * 1 / k^s for {@code k = 1,2,...,N}.
- * </pre>
- * {@code H(N,s)} is the normalizing constant
- * which corresponds to the generalized harmonic number of order N of s.
- * <ul>
- * <li>{@code N} is the number of elements</li>
- * <li>{@code s} is the exponent</li>
- * </ul>
- */
-public class ZipfDistribution extends AbstractDiscreteDistribution {
- /** Number of elements. */
- private final int numberOfElements;
- /** Exponent parameter of the distribution. */
- private final double exponent;
- /** Cached values of the nth generalized harmonic. */
- private final double nthHarmonic;
-
- /**
- * Creates a distribution.
- *
- * @param numberOfElements Number of elements.
- * @param exponent Exponent.
- * @exception IllegalArgumentException if {@code numberOfElements <= 0}
- * or {@code exponent <= 0}.
- */
- public ZipfDistribution(int numberOfElements,
- double exponent) {
- if (numberOfElements <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE,
- numberOfElements);
- }
- if (exponent <= 0) {
- throw new DistributionException(DistributionException.NEGATIVE,
- exponent);
- }
-
- this.numberOfElements = numberOfElements;
- this.exponent = exponent;
- this.nthHarmonic = generalizedHarmonic(numberOfElements, exponent);
- }
-
- /**
- * Get the number of elements (e.g. corpus size) for the distribution.
- *
- * @return the number of elements
- */
- public int getNumberOfElements() {
- return numberOfElements;
- }
-
- /**
- * Get the exponent characterizing the distribution.
- *
- * @return the exponent
- */
- public double getExponent() {
- return exponent;
- }
-
- /** {@inheritDoc} */
- @Override
- public double probability(final int x) {
- if (x <= 0 || x > numberOfElements) {
- return 0;
- }
-
- return (1 / Math.pow(x, exponent)) / nthHarmonic;
- }
-
- /** {@inheritDoc} */
- @Override
- public double logProbability(int x) {
- if (x <= 0 || x > numberOfElements) {
- return Double.NEGATIVE_INFINITY;
- }
-
- return -Math.log(x) * exponent - Math.log(nthHarmonic);
- }
-
- /** {@inheritDoc} */
- @Override
- public double cumulativeProbability(final int x) {
- if (x <= 0) {
- return 0;
- } else if (x >= numberOfElements) {
- return 1;
- }
-
- return generalizedHarmonic(x, exponent) / nthHarmonic;
- }
-
- /**
- * {@inheritDoc}
- *
- * For number of elements {@code N} and exponent {@code s}, the mean is
- * {@code Hs1 / Hs}, where
- * <ul>
- * <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
- * <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
- * </ul>
- */
- @Override
- public double getNumericalMean() {
- final int N = getNumberOfElements();
- final double s = getExponent();
-
- final double Hs1 = generalizedHarmonic(N, s - 1);
- final double Hs = nthHarmonic;
-
- return Hs1 / Hs;
- }
-
- /**
- * {@inheritDoc}
- *
- * For number of elements {@code N} and exponent {@code s}, the mean is
- * {@code (Hs2 / Hs) - (Hs1^2 / Hs^2)}, where
- * <ul>
- * <li>{@code Hs2 = generalizedHarmonic(N, s - 2)},</li>
- * <li>{@code Hs1 = generalizedHarmonic(N, s - 1)},</li>
- * <li>{@code Hs = generalizedHarmonic(N, s)}.</li>
- * </ul>
- */
- @Override
- public double getNumericalVariance() {
- final int N = getNumberOfElements();
- final double s = getExponent();
-
- final double Hs2 = generalizedHarmonic(N, s - 2);
- final double Hs1 = generalizedHarmonic(N, s - 1);
- final double Hs = nthHarmonic;
-
- return (Hs2 / Hs) - ((Hs1 * Hs1) / (Hs * Hs));
- }
-
- /**
- * Calculates the Nth generalized harmonic number. See
- * <a href="http://mathworld.wolfram.com/HarmonicSeries.html">Harmonic
- * Series</a>.
- *
- * @param n Term in the series to calculate (must be larger than 1)
- * @param m Exponent (special case {@code m = 1} is the harmonic series).
- * @return the n<sup>th</sup> generalized harmonic number.
- */
- private double generalizedHarmonic(final int n, final double m) {
- double value = 0;
- for (int k = n; k > 0; --k) {
- value += 1 / Math.pow(k, m);
- }
- return value;
- }
-
- /**
- * {@inheritDoc}
- *
- * The lower bound of the support is always 1 no matter the parameters.
- *
- * @return lower bound of the support (always 1)
- */
- @Override
- public int getSupportLowerBound() {
- return 1;
- }
-
- /**
- * {@inheritDoc}
- *
- * The upper bound of the support is the number of elements.
- *
- * @return upper bound of the support
- */
- @Override
- public int getSupportUpperBound() {
- return getNumberOfElements();
- }
-
- /**
- * {@inheritDoc}
- *
- * The support of this distribution is connected.
- *
- * @return {@code true}
- */
- @Override
- public boolean isSupportConnected() {
- return true;
- }
-
- /**{@inheritDoc} */
- @Override
- public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
- return new DiscreteDistribution.Sampler() {
- /**
- * Zipf distribution sampler.
- */
- private final DiscreteSampler sampler =
- new RejectionInversionZipfSampler(rng, numberOfElements, exponent);
-
- /**{@inheritDoc} */
- @Override
- public int sample() {
- return sampler.sample();
- }
- };
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
deleted file mode 100644
index 98315c6..0000000
--- a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/package-info.java
+++ /dev/null
@@ -1,20 +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.
- */
-/**
- * Implementations of common discrete and continuous distributions.
- */
-package org.apache.commons.statistics.distribution;
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
new file mode 100644
index 0000000..1d6b254
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractContinuousDistribution.java
@@ -0,0 +1,453 @@
+/*
+ * 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.statistics.distribution;
+
+import java.util.function.DoubleUnaryOperator;
+import org.apache.commons.numbers.core.Precision;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.InverseTransformContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.ContinuousInverseCumulativeProbabilityFunction;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+
+/**
+ * Base class for probability distributions on the reals.
+ * Default implementations are provided for some of the methods
+ * that do not vary from distribution to distribution.
+ *
+ * This base class provides a default factory method for creating
+ * a {@link ContinuousDistribution.Sampler sampler instance} that uses the
+ * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
+ * inversion method</a> for generating random samples that follow the
+ * distribution.
+ */
+public abstract class AbstractContinuousDistribution
+ implements ContinuousDistribution {
+ /**
+ * For a random variable {@code X} whose values are distributed according
+ * to this distribution, this method returns {@code P(x0 < X <= x1)}.
+ *
+ * @param x0 Lower bound (excluded).
+ * @param x1 Upper bound (included).
+ * @return the probability that a random variable with this distribution
+ * takes a value between {@code x0} and {@code x1}, excluding the lower
+ * and including the upper endpoint.
+ * @throws IllegalArgumentException if {@code x0 > x1}.
+ *
+ * The default implementation uses the identity
+ * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
+ */
+ @Override
+ public double probability(double x0,
+ double x1) {
+ if (x0 > x1) {
+ throw new DistributionException(DistributionException.TOO_LARGE, x0, x1);
+ }
+ return cumulativeProbability(x1) - cumulativeProbability(x0);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The default implementation returns
+ * <ul>
+ * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
+ * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
+ * </ul>
+ */
+ @Override
+ public double inverseCumulativeProbability(final double p) {
+ /*
+ * IMPLEMENTATION NOTES
+ * --------------------
+ * Where applicable, use is made of the one-sided Chebyshev inequality
+ * to bracket the root. This inequality states that
+ * P(X - mu >= k * sig) <= 1 / (1 + k^2),
+ * mu: mean, sig: standard deviation. Equivalently
+ * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
+ * F(mu + k * sig) >= k^2 / (1 + k^2).
+ *
+ * For k = sqrt(p / (1 - p)), we find
+ * F(mu + k * sig) >= p,
+ * and (mu + k * sig) is an upper-bound for the root.
+ *
+ * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
+ * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
+ * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
+ * P(X <= mu - k * sig) <= 1 / (1 + k^2),
+ * F(mu - k * sig) <= 1 / (1 + k^2).
+ *
+ * For k = sqrt((1 - p) / p), we find
+ * F(mu - k * sig) <= p,
+ * and (mu - k * sig) is a lower-bound for the root.
+ *
+ * In cases where the Chebyshev inequality does not apply, geometric
+ * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
+ * the root.
+ */
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+
+ double lowerBound = getSupportLowerBound();
+ if (p == 0) {
+ return lowerBound;
+ }
+
+ double upperBound = getSupportUpperBound();
+ if (p == 1) {
+ return upperBound;
+ }
+
+ final double mu = getNumericalMean();
+ final double sig = Math.sqrt(getNumericalVariance());
+ final boolean chebyshevApplies;
+ chebyshevApplies = !(Double.isInfinite(mu) ||
+ Double.isNaN(mu) ||
+ Double.isInfinite(sig) ||
+ Double.isNaN(sig));
+
+ if (lowerBound == Double.NEGATIVE_INFINITY) {
+ if (chebyshevApplies) {
+ lowerBound = mu - sig * Math.sqrt((1 - p) / p);
+ } else {
+ lowerBound = -1;
+ while (cumulativeProbability(lowerBound) >= p) {
+ lowerBound *= 2;
+ }
+ }
+ }
+
+ if (upperBound == Double.POSITIVE_INFINITY) {
+ if (chebyshevApplies) {
+ upperBound = mu + sig * Math.sqrt(p / (1 - p));
+ } else {
+ upperBound = 1;
+ while (cumulativeProbability(upperBound) < p) {
+ upperBound *= 2;
+ }
+ }
+ }
+
+ // XXX Values copied from defaults in class
+ // "o.a.c.math4.analysis.solvers.BaseAbstractUnivariateSolver"
+ final double solverRelativeAccuracy = 1e-14;
+ final double solverAbsoluteAccuracy = 1e-9;
+ final double solverFunctionValueAccuracy = 1e-15;
+
+ double x = new BrentSolver(solverRelativeAccuracy,
+ solverAbsoluteAccuracy,
+ solverFunctionValueAccuracy)
+ .solve((arg) -> cumulativeProbability(arg) - p,
+ lowerBound,
+ 0.5 * (lowerBound + upperBound),
+ upperBound);
+
+ if (!isSupportConnected()) {
+ /* Test for plateau. */
+ final double dx = solverAbsoluteAccuracy;
+ if (x - dx >= getSupportLowerBound()) {
+ double px = cumulativeProbability(x);
+ if (cumulativeProbability(x - dx) == px) {
+ upperBound = x;
+ while (upperBound - lowerBound > dx) {
+ final double midPoint = 0.5 * (lowerBound + upperBound);
+ if (cumulativeProbability(midPoint) < px) {
+ lowerBound = midPoint;
+ } else {
+ upperBound = midPoint;
+ }
+ }
+ return upperBound;
+ }
+ }
+ }
+ return x;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * @return zero.
+ */
+ @Override
+ public double probability(double x) {
+ return 0;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The default implementation computes the logarithm of {@code density(x)}.
+ */
+ @Override
+ public double logDensity(double x) {
+ return Math.log(density(x));
+ }
+
+ /**
+ * Utility function for allocating an array and filling it with {@code n}
+ * samples generated by the given {@code sampler}.
+ *
+ * @param n Number of samples.
+ * @param sampler Sampler.
+ * @return an array of size {@code n}.
+ */
+ public static double[] sample(int n,
+ ContinuousDistribution.Sampler sampler) {
+ final double[] samples = new double[n];
+ for (int i = 0; i < n; i++) {
+ samples[i] = sampler.sample();
+ }
+ return samples;
+ }
+
+ /**{@inheritDoc} */
+ @Override
+ public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new ContinuousDistribution.Sampler() {
+ /**
+ * Inversion method distribution sampler.
+ */
+ private final ContinuousSampler sampler =
+ new InverseTransformContinuousSampler(rng, createICPF());
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+
+ /**
+ * @return an instance for use by {@link #createSampler(UniformRandomProvider)}
+ */
+ private ContinuousInverseCumulativeProbabilityFunction createICPF() {
+ return new ContinuousInverseCumulativeProbabilityFunction() {
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ return AbstractContinuousDistribution.this.inverseCumulativeProbability(p);
+ }
+ };
+ }
+
+ /**
+ * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
+ * Brent algorithm</a> for finding zeros of real univariate functions.
+ * The function should be continuous but not necessarily smooth.
+ * The {@code solve} method returns a zero {@code x} of the function {@code f}
+ * in the given interval {@code [a, b]} to within a tolerance
+ * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
+ * {@code t} is the absolute accuracy.
+ * <p>The given interval must bracket the root.</p>
+ * <p>
+ * The reference implementation is given in chapter 4 of
+ * <blockquote>
+ * <b>Algorithms for Minimization Without Derivatives</b>,
+ * <em>Richard P. Brent</em>,
+ * Dover, 2002
+ * </blockquote>
+ *
+ * Used by {@link #inverseCumulativeProbability(double)}.
+ */
+ private static class BrentSolver {
+ /** Relative accuracy. */
+ private final double relativeAccuracy;
+ /** Absolutee accuracy. */
+ private final double absoluteAccuracy;
+ /** Function accuracy. */
+ private final double functionValueAccuracy;
+
+ /**
+ * Construct a solver.
+ *
+ * @param relativeAccuracy Relative accuracy.
+ * @param absoluteAccuracy Absolute accuracy.
+ * @param functionValueAccuracy Function value accuracy.
+ */
+ BrentSolver(double relativeAccuracy,
+ double absoluteAccuracy,
+ double functionValueAccuracy) {
+ this.relativeAccuracy = relativeAccuracy;
+ this.absoluteAccuracy = absoluteAccuracy;
+ this.functionValueAccuracy = functionValueAccuracy;
+ }
+
+ /**
+ * @param func Function to solve.
+ * @param min Lower bound.
+ * @param init Initial guess.
+ * @param max Upper bound.
+ * @return the root.
+ */
+ double solve(DoubleUnaryOperator func,
+ double min,
+ double initial,
+ double max) {
+ if (min > max) {
+ throw new DistributionException(DistributionException.TOO_LARGE, min, max);
+ }
+ if (initial < min ||
+ initial > max) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, initial, min, max);
+ }
+
+ // Return the initial guess if it is good enough.
+ double yInitial = func.applyAsDouble(initial);
+ if (Math.abs(yInitial) <= functionValueAccuracy) {
+ return initial;
+ }
+
+ // Return the first endpoint if it is good enough.
+ double yMin = func.applyAsDouble(min);
+ if (Math.abs(yMin) <= functionValueAccuracy) {
+ return min;
+ }
+
+ // Reduce interval if min and initial bracket the root.
+ if (yInitial * yMin < 0) {
+ return brent(func, min, initial, yMin, yInitial);
+ }
+
+ // Return the second endpoint if it is good enough.
+ double yMax = func.applyAsDouble(max);
+ if (Math.abs(yMax) <= functionValueAccuracy) {
+ return max;
+ }
+
+ // Reduce interval if initial and max bracket the root.
+ if (yInitial * yMax < 0) {
+ return brent(func, initial, max, yInitial, yMax);
+ }
+
+ throw new DistributionException(DistributionException.BRACKETING, min, yMin, max, yMax);
+ }
+
+ /**
+ * Search for a zero inside the provided interval.
+ * This implementation is based on the algorithm described at page 58 of
+ * the book
+ * <blockquote>
+ * <b>Algorithms for Minimization Without Derivatives</b>,
+ * <it>Richard P. Brent</it>,
+ * Dover 0-486-41998-3
+ * </blockquote>
+ *
+ * @param func Function to solve.
+ * @param lo Lower bound of the search interval.
+ * @param hi Higher bound of the search interval.
+ * @param fLo Function value at the lower bound of the search interval.
+ * @param fHi Function value at the higher bound of the search interval.
+ * @return the value where the function is zero.
+ */
+ private double brent(DoubleUnaryOperator func,
+ double lo, double hi,
+ double fLo, double fHi) {
+ double a = lo;
+ double fa = fLo;
+ double b = hi;
+ double fb = fHi;
+ double c = a;
+ double fc = fa;
+ double d = b - a;
+ double e = d;
+
+ final double t = absoluteAccuracy;
+ final double eps = relativeAccuracy;
+
+ while (true) {
+ if (Math.abs(fc) < Math.abs(fb)) {
+ a = b;
+ b = c;
+ c = a;
+ fa = fb;
+ fb = fc;
+ fc = fa;
+ }
+
+ final double tol = 2 * eps * Math.abs(b) + t;
+ final double m = 0.5 * (c - b);
+
+ if (Math.abs(m) <= tol ||
+ Precision.equals(fb, 0)) {
+ return b;
+ }
+ if (Math.abs(e) < tol ||
+ Math.abs(fa) <= Math.abs(fb)) {
+ // Force bisection.
+ d = m;
+ e = d;
+ } else {
+ double s = fb / fa;
+ double p;
+ double q;
+ // The equality test (a == c) is intentional,
+ // it is part of the original Brent's method and
+ // it should NOT be replaced by proximity test.
+ if (a == c) {
+ // Linear interpolation.
+ p = 2 * m * s;
+ q = 1 - s;
+ } else {
+ // Inverse quadratic interpolation.
+ q = fa / fc;
+ final double r = fb / fc;
+ p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
+ q = (q - 1) * (r - 1) * (s - 1);
+ }
+ if (p > 0) {
+ q = -q;
+ } else {
+ p = -p;
+ }
+ s = e;
+ e = d;
+ if (p >= 1.5 * m * q - Math.abs(tol * q) ||
+ p >= Math.abs(0.5 * s * q)) {
+ // Inverse quadratic interpolation gives a value
+ // in the wrong direction, or progress is slow.
+ // Fall back to bisection.
+ d = m;
+ e = d;
+ } else {
+ d = p / q;
+ }
+ }
+ a = b;
+ fa = fb;
+
+ if (Math.abs(d) > tol) {
+ b += d;
+ } else if (m > 0) {
+ b += tol;
+ } else {
+ b -= tol;
+ }
+ fb = func.applyAsDouble(b);
+ if ((fb > 0 && fc > 0) ||
+ (fb <= 0 && fc <= 0)) {
+ c = a;
+ fc = fa;
+ d = b - a;
+ e = d;
+ }
+ }
+ }
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/05626d01/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
new file mode 100644
index 0000000..faef96c
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/AbstractDiscreteDistribution.java
@@ -0,0 +1,220 @@
+/*
+ * 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.statistics.distribution;
+
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.InverseTransformDiscreteSampler;
+import org.apache.commons.rng.sampling.distribution.DiscreteInverseCumulativeProbabilityFunction;
+import org.apache.commons.rng.sampling.distribution.DiscreteSampler;
+
+/**
+ * Base class for integer-valued discrete distributions. Default
+ * implementations are provided for some of the methods that do not vary
+ * from distribution to distribution.
+ */
+public abstract class AbstractDiscreteDistribution
+ implements DiscreteDistribution {
+ /**
+ * {@inheritDoc}
+ *
+ * The default implementation uses the identity
+ * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
+ */
+ @Override
+ public double probability(int x0,
+ int x1) {
+ if (x1 < x0) {
+ throw new DistributionException(DistributionException.TOO_SMALL,
+ x1, x0);
+ }
+ return cumulativeProbability(x1) - cumulativeProbability(x0);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The default implementation returns
+ * <ul>
+ * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
+ * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li>
+ * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for
+ * {@code 0 < p < 1}.</li>
+ * </ul>
+ */
+ @Override
+ public int inverseCumulativeProbability(final double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+
+ int lower = getSupportLowerBound();
+ if (p == 0.0) {
+ return lower;
+ }
+ if (lower == Integer.MIN_VALUE) {
+ if (checkedCumulativeProbability(lower) >= p) {
+ return lower;
+ }
+ } else {
+ lower -= 1; // this ensures cumulativeProbability(lower) < p, which
+ // is important for the solving step
+ }
+
+ int upper = getSupportUpperBound();
+ if (p == 1.0) {
+ return upper;
+ }
+
+ // use the one-sided Chebyshev inequality to narrow the bracket
+ // cf. AbstractRealDistribution.inverseCumulativeProbability(double)
+ final double mu = getNumericalMean();
+ final double sigma = Math.sqrt(getNumericalVariance());
+ final boolean chebyshevApplies = !(Double.isInfinite(mu) ||
+ Double.isNaN(mu) ||
+ Double.isInfinite(sigma) ||
+ Double.isNaN(sigma) ||
+ sigma == 0.0);
+ if (chebyshevApplies) {
+ double k = Math.sqrt((1.0 - p) / p);
+ double tmp = mu - k * sigma;
+ if (tmp > lower) {
+ lower = ((int) Math.ceil(tmp)) - 1;
+ }
+ k = 1.0 / k;
+ tmp = mu + k * sigma;
+ if (tmp < upper) {
+ upper = ((int) Math.ceil(tmp)) - 1;
+ }
+ }
+
+ return solveInverseCumulativeProbability(p, lower, upper);
+ }
+
+ /**
+ * This is a utility function used by {@link
+ * #inverseCumulativeProbability(double)}. It assumes {@code 0 < p < 1} and
+ * that the inverse cumulative probability lies in the bracket {@code
+ * (lower, upper]}. The implementation does simple bisection to find the
+ * smallest {@code p}-quantile {@code inf{x in Z | P(X <= x) >= p}}.
+ *
+ * @param p Cumulative probability.
+ * @param lower Value satisfying {@code cumulativeProbability(lower) < p}.
+ * @param upper Value satisfying {@code p <= cumulativeProbability(upper)}.
+ * @return the smallest {@code p}-quantile of this distribution.
+ */
+ private int solveInverseCumulativeProbability(final double p,
+ int lower,
+ int upper) {
+ while (lower + 1 < upper) {
+ int xm = (lower + upper) / 2;
+ if (xm < lower || xm > upper) {
+ /*
+ * Overflow.
+ * There will never be an overflow in both calculation methods
+ * for xm at the same time
+ */
+ xm = lower + (upper - lower) / 2;
+ }
+
+ double pm = checkedCumulativeProbability(xm);
+ if (pm >= p) {
+ upper = xm;
+ } else {
+ lower = xm;
+ }
+ }
+ return upper;
+ }
+
+ /**
+ * Computes the cumulative probability function and checks for {@code NaN}
+ * values returned. Throws {@code MathInternalError} if the value is
+ * {@code NaN}. Rethrows any exception encountered evaluating the cumulative
+ * probability function. Throws {@code MathInternalError} if the cumulative
+ * probability function returns {@code NaN}.
+ *
+ * @param argument Input value.
+ * @return the cumulative probability.
+ * @throws IllegalStateException if the cumulative probability is {@code NaN}.
+ */
+ private double checkedCumulativeProbability(int argument) {
+ final double result = cumulativeProbability(argument);
+ if (Double.isNaN(result)) {
+ throw new IllegalStateException("Internal error");
+ }
+ return result;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The default implementation simply computes the logarithm of {@code probability(x)}.
+ */
+ @Override
+ public double logProbability(int x) {
+ return Math.log(probability(x));
+ }
+
+ /**
+ * Utility function for allocating an array and filling it with {@code n}
+ * samples generated by the given {@code sampler}.
+ *
+ * @param n Number of samples.
+ * @param sampler Sampler.
+ * @return an array of size {@code n}.
+ */
+ public static int[] sample(int n,
+ DiscreteDistribution.Sampler sampler) {
+ final int[] samples = new int[n];
+ for (int i = 0; i < n; i++) {
+ samples[i] = sampler.sample();
+ }
+ return samples;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public DiscreteDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new DiscreteDistribution.Sampler() {
+ /**
+ * Inversion method distribution sampler.
+ */
+ private final DiscreteSampler sampler =
+ new InverseTransformDiscreteSampler(rng, createICPF());
+
+ /** {@inheritDoc} */
+ @Override
+ public int sample() {
+ return sampler.sample();
+ }
+ };
+ }
+
+ /**
+ * @return an instance for use by {@link #createSampler(UniformRandomProvider)}.
+ */
+ private DiscreteInverseCumulativeProbabilityFunction createICPF() {
+ return new DiscreteInverseCumulativeProbabilityFunction() {
+ /** {@inheritDoc} */
+ @Override
+ public int inverseCumulativeProbability(double p) {
+ return AbstractDiscreteDistribution.this.inverseCumulativeProbability(p);
+ }
+ };
+ }
+}