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/21 14:05:47 UTC
[12/16] commons-statistics git commit: STATISTICS-2: Migrate
"o.a.c.math4.distribution" from Commons Math.
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java
new file mode 100644
index 0000000..1cd4a35
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DiscreteDistribution.java
@@ -0,0 +1,163 @@
+/*
+ * 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;
+
+/**
+ * Interface for distributions on the integers.
+ */
+public interface DiscreteDistribution {
+
+ /**
+ * For a random variable {@code X} whose values are distributed according to
+ * this distribution, this method returns {@code log(P(X = x))}, where
+ * {@code log} is the natural logarithm. In other words, this method
+ * represents the logarithm of the probability mass function (PMF) for the
+ * distribution. Note that due to the floating point precision and
+ * under/overflow issues, this method will for some distributions be more
+ * precise and faster than computing the logarithm of
+ * {@link #probability(int)}.
+ *
+ * @param x the point at which the PMF is evaluated
+ * @return the logarithm of the value of the probability mass function at {@code x}
+ */
+ double logProbability(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 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(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 IllegalArgumentException if {@code x0 > x1}
+ */
+ double probability(int x0, int x1);
+
+ /**
+ * 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);
+
+ /**
+ * 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}} for {@code 0 < p <= 1},</li>
+ * <li>{@code inf{x in Z | P(X<=x) > 0}} 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 IllegalArgumentException if {@code p < 0} or {@code p > 1}
+ */
+ int inverseCumulativeProbability(double p);
+
+ /**
+ * 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}}.</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}}.</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();
+
+ /**
+ * Creates a sampler.
+ *
+ * @param rng Generator of uniformly distributed numbers.
+ * @return a sampler that produces random numbers according this
+ * distribution.
+ */
+ Sampler createSampler(UniformRandomProvider rng);
+
+ /**
+ * Sampling functionality.
+ */
+ interface Sampler {
+ /**
+ * Generates a random value sampled from this distribution.
+ *
+ * @return a random value.
+ */
+ int sample();
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java
new file mode 100644
index 0000000..31d7adc
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/DistributionException.java
@@ -0,0 +1,61 @@
+/*
+ * 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.text.MessageFormat;
+
+/**
+ * Package private exception class with constants for frequently used messages.
+ */
+class DistributionException extends IllegalArgumentException {
+ /** Error message for "too large" condition. */
+ static final String TOO_LARGE = "{0} > {1}";
+ /** Error message for "too small" condition. */
+ static final String TOO_SMALL = "{0} < {1}";
+ /** Error message for "out of range" condition. */
+ static final String OUT_OF_RANGE = "Number {0} is out of range [{1}, {2}]";
+ /** Error message for "out of range" condition. */
+ static final String NEGATIVE = "Number {0} is negative";
+ /** Error message for "mismatch" condition. */
+ static final String MISMATCH = "Expected {1} but was {0}";
+ /** Error message for "failed bracketing" condition. */
+ static final String BRACKETING = "No bracketing: f({0})={1}, f({2})={3}";
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 20180119L;
+
+ /** Arguments for formatting the message. */
+ private Object[] formatArguments;
+
+ /**
+ * Create an exception where the message is constructed by applying
+ * the {@code format()} method from {@code java.text.MessageFormat}.
+ *
+ * @param message the exception message with replaceable parameters
+ * @param formatArguments the arguments for formatting the message
+ */
+ DistributionException(String message, Object... formatArguments) {
+ super(message);
+ this.formatArguments = formatArguments;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public String getMessage() {
+ return MessageFormat.format(super.getMessage(), formatArguments);
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java
new file mode 100644
index 0000000..091c6a4
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/ExponentialDistribution.java
@@ -0,0 +1,197 @@
+/*
+ * 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.AhrensDieterExponentialSampler;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Exponential_distribution">exponential distribution</a>.
+ */
+public class ExponentialDistribution extends AbstractContinuousDistribution {
+ /** The mean of this distribution. */
+ private final double mean;
+ /** The logarithm of the mean, stored to reduce computing time. */
+ private final double logMean;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mean Mean of this distribution.
+ * @throws IllegalArgumentException if {@code mean <= 0}.
+ */
+ public ExponentialDistribution(double mean) {
+ if (mean <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, mean);
+ }
+ this.mean = mean;
+ logMean = Math.log(mean);
+ }
+
+ /**
+ * Access the mean.
+ *
+ * @return the mean.
+ */
+ public double getMean() {
+ return mean;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ final double logDensity = logDensity(x);
+ return logDensity == Double.NEGATIVE_INFINITY ? 0 : Math.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>
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ double ret;
+ if (x <= 0) {
+ ret = 0;
+ } else {
+ ret = 1 - Math.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) {
+ double ret;
+
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ } else if (p == 1) {
+ ret = Double.POSITIVE_INFINITY;
+ } else {
+ ret = -mean * Math.log(1 - p);
+ }
+
+ return ret;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For mean parameter {@code k}, the mean is {@code k}.
+ */
+ @Override
+ public double getNumericalMean() {
+ return getMean();
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For mean parameter {@code k}, the variance is {@code k^2}.
+ */
+ @Override
+ 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)
+ */
+ @Override
+ 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)
+ */
+ @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;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>Sampling algorithm 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>
+ */
+ @Override
+ public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new ContinuousDistribution.Sampler() {
+ /**
+ * Exponential distribution sampler.
+ */
+ private final ContinuousSampler sampler =
+ new AhrensDieterExponentialSampler(rng, mean);
+
+ /**{@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java
new file mode 100644
index 0000000..a8b9890
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/FDistribution.java
@@ -0,0 +1,209 @@
+/*
+ * 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.LogBeta;
+import org.apache.commons.numbers.gamma.RegularizedBeta;
+
+/**
+ * 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 AbstractContinuousDistribution {
+ /** The numerator degrees of freedom. */
+ private final double numeratorDegreesOfFreedom;
+ /** The numerator degrees of freedom. */
+ private final double denominatorDegreesOfFreedom;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
+ * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
+ * @throws IllegalArgumentException if {@code numeratorDegreesOfFreedom <= 0} or
+ * {@code denominatorDegreesOfFreedom <= 0}.
+ */
+ public FDistribution(double numeratorDegreesOfFreedom,
+ double denominatorDegreesOfFreedom) {
+ if (numeratorDegreesOfFreedom <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ numeratorDegreesOfFreedom);
+ }
+ if (denominatorDegreesOfFreedom <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ denominatorDegreesOfFreedom);
+ }
+ this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
+ this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
+ }
+
+ /**
+ * {@inheritDoc}
+ */
+ @Override
+ public double density(double x) {
+ return Math.exp(logDensity(x));
+ }
+
+ /** {@inheritDoc} **/
+ @Override
+ public double logDensity(double x) {
+ final double nhalf = numeratorDegreesOfFreedom / 2;
+ final double mhalf = denominatorDegreesOfFreedom / 2;
+ final double logx = Math.log(x);
+ final double logn = Math.log(numeratorDegreesOfFreedom);
+ final double logm = Math.log(denominatorDegreesOfFreedom);
+ final double lognxm = Math.log(numeratorDegreesOfFreedom * x +
+ denominatorDegreesOfFreedom);
+ return nhalf * logn + nhalf * logx - logx +
+ mhalf * logm - nhalf * lognxm - mhalf * lognxm -
+ LogBeta.value(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>
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ double ret;
+ if (x <= 0) {
+ ret = 0;
+ } else {
+ double n = numeratorDegreesOfFreedom;
+ double m = denominatorDegreesOfFreedom;
+
+ ret = RegularizedBeta.value((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}
+ *
+ * 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>
+ */
+ @Override
+ 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>
+ */
+ @Override
+ public double getNumericalVariance() {
+ 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)
+ */
+ @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 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/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java
new file mode 100644
index 0000000..c13db5b
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GammaDistribution.java
@@ -0,0 +1,354 @@
+/*
+ * 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.LanczosApproximation;
+import org.apache.commons.numbers.gamma.RegularizedGamma;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.AhrensDieterMarsagliaTsangGammaSampler;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Gamma_distribution">Gamma distribution</a>.
+ */
+public class GammaDistribution extends AbstractContinuousDistribution {
+ /** Lanczos constant. */
+ private static final double LANCZOS_G = LanczosApproximation.g();
+ /** 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 LanczosApproximation#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 LanczosApproximation#value(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 LanczosApproximation#value(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 LanczosApproximation#value(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 LanczosApproximation#value(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;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param shape the shape parameter
+ * @param scale the scale parameter
+ * @throws IllegalArgumentException if {@code shape <= 0} or
+ * {@code scale <= 0}.
+ */
+ public GammaDistribution(double shape,
+ double scale) {
+ if (shape <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, shape);
+ }
+ if (scale <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, scale);
+ }
+
+ this.shape = shape;
+ this.scale = scale;
+ this.shiftedShape = shape + LANCZOS_G + 0.5;
+ final double aux = Math.E / (2.0 * Math.PI * shiftedShape);
+ this.densityPrefactor2 = shape * Math.sqrt(aux) / LanczosApproximation.value(shape);
+ this.logDensityPrefactor2 = Math.log(shape) + 0.5 * Math.log(aux) -
+ Math.log(LanczosApproximation.value(shape));
+ this.densityPrefactor1 = this.densityPrefactor2 / scale *
+ Math.pow(shiftedShape, -shape) * // XXX FastMath vs Math
+ Math.exp(shape + LANCZOS_G);
+ this.logDensityPrefactor1 = this.logDensityPrefactor2 - Math.log(scale) -
+ Math.log(shiftedShape) * shape +
+ shape + LANCZOS_G;
+ this.minY = shape + LANCZOS_G - Math.log(Double.MAX_VALUE);
+ this.maxLogY = Math.log(Double.MAX_VALUE) / (shape - 1.0);
+ }
+
+ /**
+ * Returns the shape parameter of {@code this} distribution.
+ *
+ * @return the shape parameter
+ */
+ public double getShape() {
+ return shape;
+ }
+
+ /**
+ * Returns the scale parameter of {@code this} distribution.
+ *
+ * @return the scale parameter
+ */
+ public double getScale() {
+ return scale;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ 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) || (Math.log(y) >= maxLogY)) {
+ /*
+ * Overflow.
+ */
+ final double aux1 = (y - shiftedShape) / shiftedShape;
+ final double aux2 = shape * (Math.log1p(aux1) - aux1); // XXX FastMath vs Math
+ final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
+ return densityPrefactor2 / x * Math.exp(aux3);
+ }
+ /*
+ * Natural calculation.
+ */
+ return densityPrefactor1 * Math.exp(-y) * Math.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) || (Math.log(y) >= maxLogY)) {
+ /*
+ * Overflow.
+ */
+ final double aux1 = (y - shiftedShape) / shiftedShape;
+ final double aux2 = shape * (Math.log1p(aux1) - aux1);
+ final double aux3 = -y * (LANCZOS_G + 0.5) / shiftedShape + LANCZOS_G + aux2;
+ return logDensityPrefactor2 - Math.log(x) + aux3;
+ }
+ /*
+ * Natural calculation.
+ */
+ return logDensityPrefactor1 - y + Math.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>
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ double ret;
+
+ if (x <= 0) {
+ ret = 0;
+ } else {
+ ret = RegularizedGamma.P.value(shape, x / scale);
+ }
+
+ return ret;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For shape parameter {@code alpha} and scale parameter {@code beta}, the
+ * mean is {@code alpha * beta}.
+ */
+ @Override
+ 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}
+ */
+ @Override
+ 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)
+ */
+ @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 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;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>
+ * Sampling algorithms:
+ * <ul>
+ * <li>
+ * For {@code 0 < shape < 1}:
+ * <blockquote>
+ * Ahrens, J. H. and Dieter, U.,
+ * <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
+ * Computing, 12, 223-246, 1974.
+ * </blockquote>
+ * </li>
+ * <li>
+ * For {@code shape >= 1}:
+ * <blockquote>
+ * Marsaglia and Tsang, <i>A Simple Method for Generating
+ * Gamma Variables.</i> ACM Transactions on Mathematical Software,
+ * Volume 26 Issue 3, September, 2000.
+ * </blockquote>
+ * </li>
+ * </ul>
+ */
+ @Override
+ public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new ContinuousDistribution.Sampler() {
+ /**
+ * Gamma distribution sampler.
+ */
+ private final ContinuousSampler sampler =
+ new AhrensDieterMarsagliaTsangGammaSampler(rng, scale, shape);
+
+ /**{@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java
new file mode 100644
index 0000000..fba2580
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GeometricDistribution.java
@@ -0,0 +1,160 @@
+/*
+ * 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 <a href="http://en.wikipedia.org/wiki/Geometric_distribution">geometric distribution</a>.
+ */
+public class GeometricDistribution extends AbstractDiscreteDistribution {
+ /** The probability of success. */
+ private final double probabilityOfSuccess;
+ /** {@code log(p)} where p is the probability of success. */
+ private final double logProbabilityOfSuccess;
+ /** {@code log(1 - p)} where p is the probability of success. */
+ private final double log1mProbabilityOfSuccess;
+
+ /**
+ * Creates a geometric distribution.
+ *
+ * @param p Probability of success.
+ * @throws IllegalArgumentException if {@code p <= 0} or {@code p > 1}.
+ */
+ public GeometricDistribution(double p) {
+ if (p <= 0 || p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+
+ probabilityOfSuccess = p;
+ logProbabilityOfSuccess = Math.log(p);
+ log1mProbabilityOfSuccess = Math.log1p(-p);
+ }
+
+ /**
+ * Access the probability of success for this distribution.
+ *
+ * @return the probability of success.
+ */
+ public double getProbabilityOfSuccess() {
+ return probabilityOfSuccess;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double probability(int x) {
+ if (x < 0) {
+ return 0.0;
+ } else {
+ return Math.exp(log1mProbabilityOfSuccess * x) * probabilityOfSuccess;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double logProbability(int x) {
+ if (x < 0) {
+ return Double.NEGATIVE_INFINITY;
+ } else {
+ return x * log1mProbabilityOfSuccess + logProbabilityOfSuccess;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(int x) {
+ if (x < 0) {
+ return 0.0;
+ } else {
+ return -Math.expm1(log1mProbabilityOfSuccess * (x + 1));
+ }
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For probability parameter {@code p}, the mean is {@code (1 - p) / p}.
+ */
+ @Override
+ public double getNumericalMean() {
+ return (1 - probabilityOfSuccess) / probabilityOfSuccess;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For probability parameter {@code p}, the variance is
+ * {@code (1 - p) / (p * p)}.
+ */
+ @Override
+ public double getNumericalVariance() {
+ return (1 - probabilityOfSuccess) / (probabilityOfSuccess * probabilityOfSuccess);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * The lower bound of the support is always 0.
+ *
+ * @return lower bound of the support (always 0)
+ */
+ @Override
+ 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)
+ */
+ @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 int inverseCumulativeProbability(double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+ if (p == 1) {
+ return Integer.MAX_VALUE;
+ }
+ if (p == 0) {
+ return 0;
+ }
+ return Math.max(0, (int) Math.ceil(Math.log1p(-p)/log1mProbabilityOfSuccess-1));
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java
new file mode 100644
index 0000000..8898b5e
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/GumbelDistribution.java
@@ -0,0 +1,128 @@
+/*
+ * 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;
+
+/**
+ * This class implements the <a href="http://en.wikipedia.org/wiki/Gumbel_distribution">Gumbel distribution</a>.
+ */
+public class GumbelDistribution extends AbstractContinuousDistribution {
+ /** π<sup>2</sup>/6. */
+ private static final double PI_SQUARED_OVER_SIX = Math.PI * Math.PI / 6;
+ /**
+ * <a href="http://mathworld.wolfram.com/Euler-MascheroniConstantApproximations.html">
+ * Approximation of Euler's constant</a>.
+ */
+ private static final double EULER = Math.PI / (2 * Math.E);
+ /** Location parameter. */
+ private final double mu;
+ /** Scale parameter. */
+ private final double beta;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mu location parameter
+ * @param beta scale parameter (must be positive)
+ * @throws IllegalArgumenException if {@code beta <= 0}
+ */
+ public GumbelDistribution(double mu,
+ double beta) {
+ if (beta <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, beta);
+ }
+
+ this.beta = beta;
+ this.mu = mu;
+ }
+
+ /**
+ * Gets the location parameter.
+ *
+ * @return the location parameter.
+ */
+ public double getLocation() {
+ return mu;
+ }
+
+ /**
+ * Gets the scale parameter.
+ *
+ * @return the scale parameter.
+ */
+ public double getScale() {
+ return beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ final double z = (x - mu) / beta;
+ final double t = Math.exp(-z);
+ return Math.exp(-z - t) / beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(double x) {
+ final double z = (x - mu) / beta;
+ return Math.exp(-Math.exp(-z));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ if (p < 0 || p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ } else if (p == 0) {
+ return Double.NEGATIVE_INFINITY;
+ } else if (p == 1) {
+ return Double.POSITIVE_INFINITY;
+ }
+ return mu - Math.log(-Math.log(p)) * beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalMean() {
+ return mu + EULER * beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalVariance() {
+ return PI_SQUARED_OVER_SIX * beta * beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportLowerBound() {
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public boolean isSupportConnected() {
+ return true;
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java
new file mode 100644
index 0000000..732a253
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/HypergeometricDistribution.java
@@ -0,0 +1,293 @@
+/*
+ * 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 <a href="http://en.wikipedia.org/wiki/Hypergeometric_distribution">hypergeometric distribution</a>.
+ */
+public class HypergeometricDistribution extends AbstractDiscreteDistribution {
+ /** 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;
+
+ /**
+ * Creates a new hypergeometric distribution.
+ *
+ * @param populationSize Population size.
+ * @param numberOfSuccesses Number of successes in the population.
+ * @param sampleSize Sample size.
+ * @throws IllegalArgumentException if {@code numberOfSuccesses < 0}, or
+ * {@code populationSize <= 0} or {@code numberOfSuccesses > populationSize},
+ * or {@code sampleSize > populationSize}.
+ */
+ public HypergeometricDistribution(int populationSize,
+ int numberOfSuccesses,
+ int sampleSize) {
+ if (populationSize <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ populationSize);
+ }
+ if (numberOfSuccesses < 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ numberOfSuccesses);
+ }
+ if (sampleSize < 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ sampleSize);
+ }
+
+ if (numberOfSuccesses > populationSize) {
+ throw new DistributionException(DistributionException.TOO_LARGE,
+ numberOfSuccesses, populationSize);
+ }
+ if (sampleSize > populationSize) {
+ throw new DistributionException(DistributionException.TOO_LARGE,
+ sampleSize, populationSize);
+ }
+
+ this.numberOfSuccesses = numberOfSuccesses;
+ this.populationSize = populationSize;
+ this.sampleSize = sampleSize;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ 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 Math.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 Math.min(k, m);
+ }
+
+ /** {@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;
+
+ 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}.
+ */
+ @Override
+ 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))}.
+ */
+ @Override
+ public double getNumericalVariance() {
+ 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
+ */
+ @Override
+ public int getSupportLowerBound() {
+ return Math.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
+ */
+ @Override
+ public int getSupportUpperBound() {
+ return Math.min(getNumberOfSuccesses(), getSampleSize());
+ }
+
+ /**
+ * {@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/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java
new file mode 100644
index 0000000..0d1a8bf
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LaplaceDistribution.java
@@ -0,0 +1,132 @@
+/*
+ * 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;
+
+/**
+ * This class implements the Laplace distribution.
+ *
+ * @see <a href="http://en.wikipedia.org/wiki/Laplace_distribution">Laplace distribution (Wikipedia)</a>
+ *
+ * @since 3.4
+ */
+public class LaplaceDistribution extends AbstractContinuousDistribution {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 20160311L;
+
+ /** The location parameter. */
+ private final double mu;
+ /** The scale parameter. */
+ private final double beta;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mu location parameter
+ * @param beta scale parameter (must be positive)
+ * @throws IllegalArgumentException if {@code beta <= 0}
+ */
+ public LaplaceDistribution(double mu,
+ double beta) {
+ if (beta <= 0.0) {
+ throw new DistributionException(DistributionException.NEGATIVE, beta);
+ }
+
+ this.mu = mu;
+ this.beta = beta;
+ }
+
+ /**
+ * 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} */
+ @Override
+ public double density(double x) {
+ return Math.exp(-Math.abs(x - mu) / beta) / (2.0 * beta);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(double x) {
+ if (x <= mu) {
+ return Math.exp((x - mu) / beta) / 2.0;
+ } else {
+ return 1.0 - Math.exp((mu - x) / beta) / 2.0;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ } else if (p == 0) {
+ return Double.NEGATIVE_INFINITY;
+ } else if (p == 1) {
+ return Double.POSITIVE_INFINITY;
+ }
+ double x = (p > 0.5) ? -Math.log(2.0 - 2.0 * p) : Math.log(2.0 * p);
+ return mu + beta * x;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalMean() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalVariance() {
+ return 2.0 * beta * beta;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportLowerBound() {
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public boolean isSupportConnected() {
+ return true;
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java
new file mode 100644
index 0000000..d16da8d
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LevyDistribution.java
@@ -0,0 +1,161 @@
+/*
+ * 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.Erfc;
+import org.apache.commons.numbers.gamma.InverseErfc;
+
+/**
+ * This class implements the <a href="http://en.wikipedia.org/wiki/L%C3%A9vy_distribution">
+ * Lévy distribution</a>.
+ */
+public class LevyDistribution extends AbstractContinuousDistribution {
+ /** Location parameter. */
+ private final double mu;
+ /** Scale parameter. */
+ private final double c;
+ /** Half of c (for calculations). */
+ private final double halfC;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mu location
+ * @param c scale parameter
+ */
+ public LevyDistribution(final double mu,
+ final double c) {
+ this.mu = mu;
+ this.c = c;
+ this.halfC = 0.5 * c;
+ }
+
+ /** {@inheritDoc}
+ * <p>
+ * From Wikipedia: The probability density function of the Lévy distribution
+ * over the domain is
+ * </p>
+ * <div style="white-space: pre"><code>
+ * f(x; μ, c) = √(c / 2π) * e<sup>-c / 2 (x - μ)</sup> / (x - μ)<sup>3/2</sup>
+ * </code></div>
+ * <p>
+ * For this distribution, {@code X}, this method returns {@code P(X < x)}.
+ * If {@code x} is less than location parameter μ, {@code Double.NaN} is
+ * returned, as in these cases the distribution is not defined.
+ * </p>
+ */
+ @Override
+ public double density(final double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+
+ final double delta = x - mu;
+ final double f = halfC / delta;
+ return Math.sqrt(f / Math.PI) * Math.exp(-f) /delta;
+ }
+
+ /** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+
+ final double delta = x - mu;
+ final double f = halfC / delta;
+ return 0.5 * Math.log(f / Math.PI) - f - Math.log(delta);
+ }
+
+ /** {@inheritDoc}
+ * <p>
+ * From Wikipedia: the cumulative distribution function is
+ * </p>
+ * <pre>
+ * f(x; u, c) = erfc (√ (c / 2 (x - u )))
+ * </pre>
+ */
+ @Override
+ public double cumulativeProbability(final double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+ return Erfc.value(Math.sqrt(halfC / (x - mu)));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(final double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ }
+ final double t = InverseErfc.value(p);
+ return mu + halfC / (t * t);
+ }
+
+ /**
+ * Gets the scale parameter of the distribution.
+ *
+ * @return scale parameter of the distribution
+ */
+ public double getScale() {
+ return c;
+ }
+
+ /**
+ * Gets the location parameter of the distribution.
+ *
+ * @return location parameter of the distribution
+ */
+ public double getLocation() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalMean() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalVariance() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportLowerBound() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public boolean isSupportConnected() {
+ return true;
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java
new file mode 100644
index 0000000..25bdd33
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogNormalDistribution.java
@@ -0,0 +1,266 @@
+/*
+ * 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.Erf;
+import org.apache.commons.numbers.gamma.ErfDifference;
+import org.apache.commons.rng.UniformRandomProvider;
+import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
+import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
+import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
+
+/**
+ * Implementation of the <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">log-normal distribution</a>.
+ *
+ * <p>
+ * <strong>Parameters:</strong>
+ * {@code X} is log-normally distributed if its natural logarithm {@code log(X)}
+ * is normally distributed. The probability distribution function of {@code X}
+ * is given by (for {@code x > 0})
+ * </p>
+ * <p>
+ * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
+ * </p>
+ * <ul>
+ * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the
+ * normally distributed natural logarithm of this distribution,</li>
+ * <li>{@code s} is the <em>shape</em> parameter: this is the standard
+ * deviation of the normally distributed natural logarithm of this
+ * distribution.
+ * </ul>
+ */
+public class LogNormalDistribution extends AbstractContinuousDistribution {
+ /** √(2 π) */
+ private static final double SQRT2PI = Math.sqrt(2 * Math.PI);
+ /** √(2) */
+ private static final double SQRT2 = Math.sqrt(2);
+ /** The scale parameter of this distribution. */
+ private final double scale;
+ /** The shape parameter of this distribution. */
+ private final double shape;
+ /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */
+ private final double logShapePlusHalfLog2Pi;
+
+ /**
+ * Creates a log-normal distribution, where the mean and standard deviation
+ * of the {@link NormalDistribution normally distributed} natural
+ * logarithm of the log-normal distribution are equal to zero and one
+ * respectively. In other words, the scale of the returned distribution is
+ * {@code 0}, while its shape is {@code 1}.
+ */
+ public LogNormalDistribution() {
+ this(0, 1);
+ }
+
+ /**
+ * Creates a log-normal distribution.
+ *
+ * @param scale Scale parameter of this distribution.
+ * @param shape Shape parameter of this distribution.
+ * @throws IllegalArgumentException if {@code shape <= 0}.
+ */
+ public LogNormalDistribution(double scale,
+ double shape) {
+ if (shape <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, shape);
+ }
+
+ this.scale = scale;
+ this.shape = shape;
+ this.logShapePlusHalfLog2Pi = Math.log(shape) + 0.5 * Math.log(2 * Math.PI);
+ }
+
+ /**
+ * Returns the scale parameter of this distribution.
+ *
+ * @return the scale parameter
+ */
+ public double getScale() {
+ return scale;
+ }
+
+ /**
+ * Returns the shape parameter of this distribution.
+ *
+ * @return the shape parameter
+ */
+ public double getShape() {
+ return shape;
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For scale {@code m}, and shape {@code s} of this distribution, the PDF
+ * is given by
+ * <ul>
+ * <li>{@code 0} if {@code x <= 0},</li>
+ * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
+ * otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double density(double x) {
+ if (x <= 0) {
+ return 0;
+ }
+ final double x0 = Math.log(x) - scale;
+ final double x1 = x0 / shape;
+ return Math.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
+ }
+
+ /** {@inheritDoc}
+ *
+ * See documentation of {@link #density(double)} for computation details.
+ */
+ @Override
+ public double logDensity(double x) {
+ if (x <= 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
+ final double logX = Math.log(x);
+ final double x0 = logX - scale;
+ final double x1 = x0 / shape;
+ return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For scale {@code m}, and shape {@code s} of this distribution, the CDF
+ * is given by
+ * <ul>
+ * <li>{@code 0} if {@code x <= 0},</li>
+ * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as
+ * in these cases the actual value is within {@code Double.MIN_VALUE} of 0,
+ * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s},
+ * as in these cases the actual value is within {@code Double.MIN_VALUE} of
+ * 1,</li>
+ * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li>
+ * </ul>
+ */
+ @Override
+ public double cumulativeProbability(double x) {
+ if (x <= 0) {
+ return 0;
+ }
+ final double dev = Math.log(x) - scale;
+ if (Math.abs(dev) > 40 * shape) {
+ return dev < 0 ? 0.0d : 1.0d;
+ }
+ return 0.5 + 0.5 * Erf.value(dev / (shape * SQRT2));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double probability(double x0,
+ double x1) {
+ if (x0 > x1) {
+ throw new DistributionException(DistributionException.TOO_LARGE,
+ x0, x1);
+ }
+ if (x0 <= 0 || x1 <= 0) {
+ return super.probability(x0, x1);
+ }
+ final double denom = shape * SQRT2;
+ final double v0 = (Math.log(x0) - scale) / denom;
+ final double v1 = (Math.log(x1) - scale) / denom;
+ return 0.5 * ErfDifference.value(v0, v1);
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For scale {@code m} and shape {@code s}, the mean is
+ * {@code exp(m + s^2 / 2)}.
+ */
+ @Override
+ public double getNumericalMean() {
+ double s = shape;
+ return Math.exp(scale + (s * s / 2));
+ }
+
+ /**
+ * {@inheritDoc}
+ *
+ * For scale {@code m} and shape {@code s}, the variance is
+ * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}.
+ */
+ @Override
+ public double getNumericalVariance() {
+ final double s = shape;
+ final double ss = s * s;
+ return (Math.expm1(ss)) * Math.exp(2 * scale + ss);
+ }
+
+ /**
+ * {@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;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new ContinuousDistribution.Sampler() {
+ /**
+ * Log normal distribution sampler.
+ */
+ private final ContinuousSampler sampler =
+ new LogNormalSampler(new ZigguratNormalizedGaussianSampler(rng), scale, shape);
+
+ /**{@inheritDoc} */
+ @Override
+ public double sample() {
+ return sampler.sample();
+ }
+ };
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java
new file mode 100644
index 0000000..28a6657
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/LogisticDistribution.java
@@ -0,0 +1,128 @@
+/*
+ * 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 <a href="http://en.wikipedia.org/wiki/Logistic_distribution">Logistic distribution</a>.
+ */
+public class LogisticDistribution extends AbstractContinuousDistribution {
+ /** π<sup>2</sup>/3. */
+ private static final double PI_SQUARED_OVER_THREE = Math.PI * Math.PI / 3;
+ /** Location parameter. */
+ private final double mu;
+ /** Scale parameter. */
+ private final double scale;
+ /** Inverse of "scale". */
+ private final double oneOverScale;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mu Location parameter.
+ * @param scale Scale parameter (must be positive).
+ * @throws IllegalArgumentException if {@code scale <= 0}.
+ */
+ public LogisticDistribution(double mu,
+ double scale) {
+ if (scale <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE,
+ scale);
+ }
+
+ this.mu = mu;
+ this.scale = scale;
+ this.oneOverScale = 1 / scale;
+ }
+
+ /**
+ * Gets the location parameter.
+ *
+ * @return the location parameter.
+ */
+ public double getLocation() {
+ return mu;
+ }
+
+ /**
+ * Gets the scale parameter.
+ *
+ * @return the scale parameter.
+ */
+ public double getScale() {
+ return scale;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ final double z = oneOverScale * (x - mu);
+ final double v = Math.exp(-z);
+ return oneOverScale * v / ((1 + v) * (1 + v));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(double x) {
+ final double z = oneOverScale * (x - mu);
+ return 1 / (1 + Math.exp(-z));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(double p) {
+ if (p < 0 ||
+ p > 1) {
+ throw new DistributionException(DistributionException.OUT_OF_RANGE, p, 0, 1);
+ } else if (p == 0) {
+ return 0;
+ } else if (p == 1) {
+ return Double.POSITIVE_INFINITY;
+ } else {
+ return scale * Math.log(p / (1 - p)) + mu;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalMean() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalVariance() {
+ return oneOverScale * oneOverScale * PI_SQUARED_OVER_THREE;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportLowerBound() {
+ return Double.NEGATIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public boolean isSupportConnected() {
+ return true;
+ }
+}
http://git-wip-us.apache.org/repos/asf/commons-statistics/blob/9c794a15/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java
----------------------------------------------------------------------
diff --git a/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java
new file mode 100644
index 0000000..9bf7d2f
--- /dev/null
+++ b/commons-statistics-distribution/src/main/java/commons/statistics/distribution/NakagamiDistribution.java
@@ -0,0 +1,117 @@
+/*
+ * 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.Gamma;
+import org.apache.commons.numbers.gamma.RegularizedGamma;
+
+/**
+ * This class implements the <a href="http://en.wikipedia.org/wiki/Nakagami_distribution">Nakagami distribution</a>.
+ */
+public class NakagamiDistribution extends AbstractContinuousDistribution {
+ /** The shape parameter. */
+ private final double mu;
+ /** The scale parameter. */
+ private final double omega;
+
+ /**
+ * Creates a distribution.
+ *
+ * @param mu shape parameter
+ * @param omega scale parameter (must be positive)
+ * @throws IllegalArgumentException if {@code mu < 0.5} or if
+ * {@code omega <= 0}.
+ */
+ public NakagamiDistribution(double mu,
+ double omega) {
+ if (mu < 0.5) {
+ throw new DistributionException(DistributionException.TOO_SMALL, mu, 0.5);
+ }
+ if (omega <= 0) {
+ throw new DistributionException(DistributionException.NEGATIVE, omega);
+ }
+
+ this.mu = mu;
+ this.omega = omega;
+ }
+
+ /**
+ * Access the shape parameter, {@code mu}.
+ *
+ * @return the shape parameter.
+ */
+ public double getShape() {
+ return mu;
+ }
+
+ /**
+ * Access the scale parameter, {@code omega}.
+ *
+ * @return the scale parameter.
+ */
+ public double getScale() {
+ return omega;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double density(double x) {
+ if (x <= 0) {
+ return 0.0;
+ }
+ return 2.0 * Math.pow(mu, mu) / (Gamma.value(mu) * Math.pow(omega, mu)) *
+ Math.pow(x, 2 * mu - 1) * Math.exp(-mu * x * x / omega);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double cumulativeProbability(double x) {
+ return RegularizedGamma.P.value(mu, mu * x * x / omega);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalMean() {
+ return Gamma.value(mu + 0.5) / Gamma.value(mu) * Math.sqrt(omega / mu);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getNumericalVariance() {
+ double v = Gamma.value(mu + 0.5) / Gamma.value(mu);
+ return omega * (1 - 1 / mu * v * v);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportLowerBound() {
+ return 0;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public boolean isSupportConnected() {
+ return true;
+ }
+
+}