/*
* 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;
}
}