Normal distribution now uses a direct implementation of the
inverse error function to compute inverse cumulative probability
instead of relying on a numerical solver. This is much faster,
diff --git a/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java b/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
new file mode 100644
index 000000000..01bf03a48
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/distribution/LevyDistribution.java
@@ -0,0 +1,143 @@
+package org.apache.commons.math3.distribution;
+
+import org.apache.commons.math3.exception.OutOfRangeException;
+import org.apache.commons.math3.random.RandomGenerator;
+import org.apache.commons.math3.special.Erf;
+import org.apache.commons.math3.util.FastMath;
+
+/**
+ * This class implements the
+ * Lévy distribution.
+ *
+ * @version $Id$
+ * @since 3.2
+ */
+public class LevyDistribution extends AbstractRealDistribution {
+
+ /** Serializable UID. */
+ private static final long serialVersionUID = 20130314L;
+
+ /** Location parameter. */
+ private final double mu;
+
+ /** Scale parameter. */
+ private final double c; // Setting this to 1 returns a cumProb of 1.0
+
+ /** Half of c (for calculations). */
+ private final double halfC;
+
+ /**
+ * Creates a LevyDistribution.
+ * @param rng random generator to be used for sampling
+ * @param mu location
+ * @param c scale parameter
+ */
+ public LevyDistribution(final RandomGenerator rng, final double mu, final double c) {
+ super(rng);
+ this.mu = mu;
+ this.c = c;
+ this.halfC = 0.5 * c;
+ }
+
+ /** {@inheritDoc}
+ *
+ * From Wikipedia: The probability density function of the Lévy distribution
+ * over the domain is
+ *
+ *
+ * f(x; μ, c) = √(c / 2π) * e-c / 2 (x - μ) / (x - μ)3/2
+ *
+ *
+ * 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.
+ *
+ */
+ public double density(final double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+
+ final double delta = x - mu;
+ final double f = halfC / delta;
+ return FastMath.sqrt(f / FastMath.PI) * FastMath.exp(-f) /delta;
+ }
+
+ /** {@inheritDoc}
+ *
+ * From Wikipedia: the cumulative distribution function is
+ *
+ *
+ * f(x; u, c) = erfc (√ (c / 2 (x - u )))
+ *
+ */
+ public double cumulativeProbability(final double x) {
+ if (x < mu) {
+ return Double.NaN;
+ }
+ return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
+ if (p < 0.0 || p > 1.0) {
+ throw new OutOfRangeException(p, 0, 1);
+ }
+ final double t = Erf.erfcInv(p);
+ return mu + halfC / (t * t);
+ }
+
+ /** Get the scale parameter of the distribution.
+ * @return scale parameter of the distribution
+ */
+ public double getScale() {
+ return c;
+ }
+
+ /** Get the location parameter of the distribution.
+ * @return location parameter of the distribution
+ */
+ public double getLocation() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ public double getNumericalMean() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ public double getNumericalVariance() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ public double getSupportLowerBound() {
+ return mu;
+ }
+
+ /** {@inheritDoc} */
+ public double getSupportUpperBound() {
+ return Double.POSITIVE_INFINITY;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isSupportLowerBoundInclusive() {
+ // there is a division by x-mu in the computation, so density
+ // is not finite at lower bound, bound must be excluded
+ return false;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isSupportUpperBoundInclusive() {
+ // upper bound is infinite, so it must be excluded
+ return false;
+ }
+
+ /** {@inheritDoc} */
+ public boolean isSupportConnected() {
+ return true;
+ }
+
+}
diff --git a/src/test/R/LevyDistributionTestCases.R b/src/test/R/LevyDistributionTestCases.R
new file mode 100644
index 000000000..7a6651f4b
--- /dev/null
+++ b/src/test/R/LevyDistributionTestCases.R
@@ -0,0 +1,89 @@
+# 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.
+#
+#------------------------------------------------------------------------------
+# R source file to validate Lévy distribution tests in
+# org.apache.commons.math3.distribution.LevyDistributionTest
+#
+# To run the test, install R, put this file and testFunctions
+# into the same directory, launch R from this directory and then enter
+# source("")
+#
+# R functions used
+# dlevy(q, m=0, s=1, log=FALSE)
+# plevy(q, m=0, s=1)
+# qlevy(p, m=0, s=1)
+#-----------------------------------------------------------------------------
+tol <- 1E-9
+
+# Function definitions
+
+source("testFunctions") # utility test functions
+library(rmutil)
+
+# function to verify distribution computations
+
+verifyDistribution <- function(points, expected, m, s, tol) {
+ rDistValues <- rep(0, length(points))
+ i <- 0
+ for (point in points) {
+ i <- i + 1
+ rDistValues[i] <- plevy(point, m, s)
+ }
+ output <- c("Distribution test m = ",m,", s = ", s)
+ if (assertEquals(expected, rDistValues, tol, "Distribution Values")) {
+ displayPadded(output, SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(output, FAILED, WIDTH)
+ }
+}
+
+# function to verify density computations
+
+verifyDensity <- function(points, expected, m, s, tol) {
+ rDensityValues <- rep(0, length(points))
+ i <- 0
+ for (point in points) {
+ i <- i + 1
+ rDensityValues[i] <- dlevy(point, m, s, log=FALSE)
+ }
+ output <- c("Density test m = ",m,", s = ", s)
+ if (assertEquals(expected, rDensityValues, tol, "Density Values")) {
+ displayPadded(output, SUCCEEDED, WIDTH)
+ } else {
+ displayPadded(output, FAILED, WIDTH)
+ }
+}
+
+#--------------------------------------------------------------------------
+cat("Levy test cases\n")
+
+m <- 1.2
+s <- 0.4
+distributionPoints <- c(1.2001, 1.21, 1.225, 1.25, 1.3, 1.9, 3.4, 5.6)
+densityValues <- c(0.0, 5.200563737654472E-7, 0.021412836122382383, 0.4133397070818418, 1.0798193302637613, 0.3237493191610873, 0.07060325500936372, 0.026122839883975738)
+distributionValues <- c(0.0, 2.539628589470901E-10, 6.334248366624259E-5, 0.004677734981047284, 0.04550026389635843, 0.4496917979688907, 0.6698153575994166, 0.763024600552995)
+verifyDistribution(distributionPoints, distributionValues, m, s, tol)
+verifyDensity(distributionPoints, densityValues, m, s, tol)
+
+m <- 5
+s <- 1.3
+distributionPoints <- c(5.0001, 6, 7, 8, 9, 10, 11, 12, 13, 14)
+densityValues <- c(0.0, 0.23745992633364185, 0.1161959636020616, 0.07048597672583455, 0.04833023442399538, 0.03572468867742048, 0.02777194506550441, 0.022382435270909086, 0.018533623436073274, 0.0156730047506865)
+distributionValues <- c(0.0, 0.25421322360396437, 0.42011267955064, 0.5103578488686935, 0.5686182086635944, 0.6101201547975077, 0.6415915735304425, 0.6665077778509312, 0.6868651803414656, 0.7039020091632311)
+verifyDistribution(distributionPoints, distributionValues, m, s, tol)
+verifyDensity(distributionPoints, densityValues, m, s, tol)
+
+displayDashes(WIDTH)
diff --git a/src/test/R/testAll b/src/test/R/testAll
index ba5d22388..ac3cfb670 100644
--- a/src/test/R/testAll
+++ b/src/test/R/testAll
@@ -40,6 +40,7 @@ source("FDistributionTestCases.R")
source("GammaDistributionTestCases.R")
source("WeibullDistributionTestCases.R")
source("ChiSquareDistributionTestCases.R")
+source("LevyDistributionTestCases.R")
# regression
source("regressionTestCases")
diff --git a/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java b/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
new file mode 100644
index 000000000..12a0a338a
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/distribution/LevyDistributionTest.java
@@ -0,0 +1,51 @@
+package org.apache.commons.math3.distribution;
+
+import org.apache.commons.math3.random.Well19937a;
+import org.apache.commons.math3.util.Precision;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class LevyDistributionTest extends RealDistributionAbstractTest {
+
+ @Test
+ public void testParameters() {
+ LevyDistribution d = makeDistribution();
+ Assert.assertEquals(1.2, d.getLocation(), Precision.EPSILON);
+ Assert.assertEquals(0.4, d.getScale(), Precision.EPSILON);
+ }
+
+ @Test
+ public void testSupport() {
+ LevyDistribution d = makeDistribution();
+ Assert.assertEquals(d.getLocation(), d.getSupportLowerBound(), Precision.EPSILON);
+ Assert.assertTrue(Double.isInfinite(d.getSupportUpperBound()));
+ Assert.assertTrue(d.isSupportConnected());
+ }
+
+ public LevyDistribution makeDistribution() {
+ return new LevyDistribution(new Well19937a(0xc5a5506bbb17e57al), 1.2, 0.4);
+ }
+
+ public double[] makeCumulativeTestPoints() {
+ return new double[] {
+ 1.2001, 1.21, 1.225, 1.25, 1.3, 1.9, 3.4, 5.6
+ };
+ }
+
+ public double[] makeCumulativeTestValues() {
+ // values computed with R and function plevy from rmutil package
+ return new double[] {
+ 0, 2.53962850749e-10, 6.33424836662e-05, 0.00467773498105,
+ 0.0455002638964, 0.449691797969, 0.669815357599, 0.763024600553
+ };
+ }
+
+ public double[] makeDensityTestValues() {
+ // values computed with R and function dlevy from rmutil package
+ return new double[] {
+ 0, 5.20056373765e-07, 0.0214128361224, 0.413339707082, 1.07981933026,
+ 0.323749319161, 0.0706032550094, 0.026122839884
+ };
+ }
+
+}