001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.statistics.distribution;
018
019import org.apache.commons.numbers.gamma.Erfc;
020import org.apache.commons.numbers.gamma.InverseErfc;
021
022/**
023 * This class implements the <a href="http://en.wikipedia.org/wiki/L%C3%A9vy_distribution">
024 * L&eacute;vy distribution</a>.
025 */
026public class LevyDistribution extends AbstractContinuousDistribution {
027    /** Location parameter. */
028    private final double mu;
029    /** Scale parameter. */
030    private final double c;
031    /** Half of c (for calculations). */
032    private final double halfC;
033
034    /**
035     * Creates a distribution.
036     *
037     * @param mu location
038     * @param c scale parameter
039     */
040    public LevyDistribution(final double mu,
041                            final double c) {
042        this.mu = mu;
043        this.c = c;
044        this.halfC = 0.5 * c;
045    }
046
047    /** {@inheritDoc}
048    * <p>
049    * From Wikipedia: The probability density function of the L&eacute;vy distribution
050    * over the domain is
051    * </p>
052    * <div style="white-space: pre"><code>
053    * f(x; &mu;, c) = &radic;(c / 2&pi;) * e<sup>-c / 2 (x - &mu;)</sup> / (x - &mu;)<sup>3/2</sup>
054    * </code></div>
055    * <p>
056    * For this distribution, {@code X}, this method returns {@code P(X < x)}.
057    * If {@code x} is less than location parameter &mu;, {@code Double.NaN} is
058    * returned, as in these cases the distribution is not defined.
059    * </p>
060    */
061    @Override
062    public double density(final double x) {
063        if (x < mu) {
064            return 0;
065        }
066
067        final double delta = x - mu;
068        final double f = halfC / delta;
069        return Math.sqrt(f / Math.PI) * Math.exp(-f) / delta;
070    }
071
072    /** {@inheritDoc}
073     *
074     * See documentation of {@link #density(double)} for computation details.
075     */
076    @Override
077    public double logDensity(double x) {
078        if (x < mu) {
079            return Double.NEGATIVE_INFINITY;
080        }
081
082        final double delta = x - mu;
083        final double f     = halfC / delta;
084        return 0.5 * Math.log(f / Math.PI) - f - Math.log(delta);
085    }
086
087    /** {@inheritDoc}
088     * <p>
089     * From Wikipedia: the cumulative distribution function is
090     * </p>
091     * <pre>
092     * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
093     * </pre>
094     */
095    @Override
096    public double cumulativeProbability(final double x) {
097        if (x < mu) {
098            return 0;
099        }
100        return Erfc.value(Math.sqrt(halfC / (x - mu)));
101    }
102
103    /** {@inheritDoc} */
104    @Override
105    public double inverseCumulativeProbability(final double p) {
106        if (p < 0 ||
107            p > 1) {
108            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
109        }
110        final double t = InverseErfc.value(p);
111        return mu + halfC / (t * t);
112    }
113
114    /**
115     * Gets the scale parameter of the distribution.
116     *
117     * @return scale parameter of the distribution
118     */
119    public double getScale() {
120        return c;
121    }
122
123    /**
124     * Gets the location parameter of the distribution.
125     *
126     * @return location parameter of the distribution
127     */
128    public double getLocation() {
129        return mu;
130    }
131
132    /** {@inheritDoc} */
133    @Override
134    public double getMean() {
135        return Double.POSITIVE_INFINITY;
136    }
137
138    /** {@inheritDoc} */
139    @Override
140    public double getVariance() {
141        return Double.POSITIVE_INFINITY;
142    }
143
144    /** {@inheritDoc} */
145    @Override
146    public double getSupportLowerBound() {
147        return getLocation();
148    }
149
150    /** {@inheritDoc} */
151    @Override
152    public double getSupportUpperBound() {
153        return Double.POSITIVE_INFINITY;
154    }
155
156    /** {@inheritDoc} */
157    @Override
158    public boolean isSupportConnected() {
159        return true;
160    }
161}