Added the Lévy distribution.

The patch from Andrew Waterman has been added, with several fixes and
extensions.

JIRA: MATH-460

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1456907 13f79535-47bb-0310-9956-ffa450edef68

pom.xml

@@ -276,6 +276,9 @@
     <contributor>
       <name>Kim van der Linde</name>
     </contributor>
+    <contributor>
+      <name>Andrew Waterman</name>
+    </contributor>
     <contributor>
       <name>J&#246;rg Weimar</name>
     </contributor>

src/changes/changes.xml

@@ -55,7 +55,10 @@ This is a minor release: It combines bug fixes and new features.
   Changes to existing features were made in a backwards-compatible
   way such as to allow drop-in replacement of the v3.1[.1] JAR file.
 ">
-      <action dev="luc" type="update"  >
+      <action dev="luc" type="add" issue="MATH-460" due-to="Andrew Waterman">
+        Added the Lévy distribution.
+      </action>
+        <action dev="luc" type="update"  >
         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,

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 <a href="http://en.wikipedia.org/wiki/L%C3%A9vy_distribution">
+ * L&eacute;vy distribution</a>.
+ *
+ * @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}
+    * <p>
+    * From Wikipedia: The probability density function of the L&eacute;vy distribution
+    * over the domain is
+    * </p>
+    * <pre>
+    * f(x; &mu;, c) = &radic;(c / 2&pi;) * e<sup>-c / 2 (x - &mu;)</sup> / (x - &mu;)<sup>3/2</sup>
+    * </pre>
+    * <p>
+    * For this distribution, {@code X}, this method returns {@code P(X < x)}.
+    * If {@code x} is less than location parameter &mu;, {@code Double.NaN} is
+    * returned, as in these cases the distribution is not defined.
+    * </p>
+    */
+    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}
+     * <p>
+     * From Wikipedia: the cumulative distribution function is
+     * </p>
+     * <pre>
+     * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
+     * </pre>
+     */
+    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;
+    }
+
+}

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("<name-of-this-file>")
+#
+# 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)

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")

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
+        };
+    }
+
+}