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 */
017
018package org.apache.commons.statistics.distribution;
019
020import org.apache.commons.numbers.gamma.Erf;
021import org.apache.commons.numbers.gamma.ErfDifference;
022import org.apache.commons.rng.UniformRandomProvider;
023import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
024import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
025
026/**
027 * Implementation of the <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">log-normal distribution</a>.
028 *
029 * <p>
030 * <strong>Parameters:</strong>
031 * {@code X} is log-normally distributed if its natural logarithm {@code log(X)}
032 * is normally distributed. The probability distribution function of {@code X}
033 * is given by (for {@code x > 0})
034 * </p>
035 * <p>
036 * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
037 * </p>
038 * <ul>
039 * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the
040 * normally distributed natural logarithm of this distribution,</li>
041 * <li>{@code s} is the <em>shape</em> parameter: this is the standard
042 * deviation of the normally distributed natural logarithm of this
043 * distribution.
044 * </ul>
045 */
046public class LogNormalDistribution extends AbstractContinuousDistribution {
047    /** &radic;(2 &pi;). */
048    private static final double SQRT2PI = Math.sqrt(2 * Math.PI);
049    /** &radic;(2). */
050    private static final double SQRT2 = Math.sqrt(2);
051    /** The scale parameter of this distribution. */
052    private final double scale;
053    /** The shape parameter of this distribution. */
054    private final double shape;
055    /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */
056    private final double logShapePlusHalfLog2Pi;
057
058    /**
059     * Creates a log-normal distribution.
060     *
061     * @param scale Scale parameter of this distribution.
062     * @param shape Shape parameter of this distribution.
063     * @throws IllegalArgumentException if {@code shape <= 0}.
064     */
065    public LogNormalDistribution(double scale,
066                                 double shape) {
067        if (shape <= 0) {
068            throw new DistributionException(DistributionException.NEGATIVE, shape);
069        }
070
071        this.scale = scale;
072        this.shape = shape;
073        this.logShapePlusHalfLog2Pi = Math.log(shape) + 0.5 * Math.log(2 * Math.PI);
074    }
075
076    /**
077     * Returns the scale parameter of this distribution.
078     *
079     * @return the scale parameter
080     */
081    public double getScale() {
082        return scale;
083    }
084
085    /**
086     * Returns the shape parameter of this distribution.
087     *
088     * @return the shape parameter
089     */
090    public double getShape() {
091        return shape;
092    }
093
094    /**
095     * {@inheritDoc}
096     *
097     * For scale {@code m}, and shape {@code s} of this distribution, the PDF
098     * is given by
099     * <ul>
100     * <li>{@code 0} if {@code x <= 0},</li>
101     * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
102     * otherwise.</li>
103     * </ul>
104     */
105    @Override
106    public double density(double x) {
107        if (x <= 0) {
108            return 0;
109        }
110        final double x0 = Math.log(x) - scale;
111        final double x1 = x0 / shape;
112        return Math.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
113    }
114
115    /** {@inheritDoc}
116     *
117     * See documentation of {@link #density(double)} for computation details.
118     */
119    @Override
120    public double logDensity(double x) {
121        if (x <= 0) {
122            return Double.NEGATIVE_INFINITY;
123        }
124        final double logX = Math.log(x);
125        final double x0 = logX - scale;
126        final double x1 = x0 / shape;
127        return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX);
128    }
129
130    /**
131     * {@inheritDoc}
132     *
133     * For scale {@code m}, and shape {@code s} of this distribution, the CDF
134     * is given by
135     * <ul>
136     * <li>{@code 0} if {@code x <= 0},</li>
137     * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as
138     * in these cases the actual value is within {@code Double.MIN_VALUE} of 0,
139     * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s},
140     * as in these cases the actual value is within {@code Double.MIN_VALUE} of
141     * 1,</li>
142     * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li>
143     * </ul>
144     */
145    @Override
146    public double cumulativeProbability(double x)  {
147        if (x <= 0) {
148            return 0;
149        }
150        final double dev = Math.log(x) - scale;
151        if (Math.abs(dev) > 40 * shape) {
152            return dev < 0 ? 0.0d : 1.0d;
153        }
154        return 0.5 + 0.5 * Erf.value(dev / (shape * SQRT2));
155    }
156
157    /** {@inheritDoc} */
158    @Override
159    public double probability(double x0,
160                              double x1) {
161        if (x0 > x1) {
162            throw new DistributionException(DistributionException.TOO_LARGE,
163                                            x0, x1);
164        }
165        if (x0 <= 0 || x1 <= 0) {
166            return super.probability(x0, x1);
167        }
168        final double denom = shape * SQRT2;
169        final double v0 = (Math.log(x0) - scale) / denom;
170        final double v1 = (Math.log(x1) - scale) / denom;
171        return 0.5 * ErfDifference.value(v0, v1);
172    }
173
174    /**
175     * {@inheritDoc}
176     *
177     * For scale {@code m} and shape {@code s}, the mean is
178     * {@code exp(m + s^2 / 2)}.
179     */
180    @Override
181    public double getMean() {
182        final double s = shape;
183        return Math.exp(scale + (s * s / 2));
184    }
185
186    /**
187     * {@inheritDoc}
188     *
189     * For scale {@code m} and shape {@code s}, the variance is
190     * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}.
191     */
192    @Override
193    public double getVariance() {
194        final double s = shape;
195        final double ss = s * s;
196        return (Math.expm1(ss)) * Math.exp(2 * scale + ss);
197    }
198
199    /**
200     * {@inheritDoc}
201     *
202     * The lower bound of the support is always 0 no matter the parameters.
203     *
204     * @return lower bound of the support (always 0)
205     */
206    @Override
207    public double getSupportLowerBound() {
208        return 0;
209    }
210
211    /**
212     * {@inheritDoc}
213     *
214     * The upper bound of the support is always positive infinity
215     * no matter the parameters.
216     *
217     * @return upper bound of the support (always
218     * {@code Double.POSITIVE_INFINITY})
219     */
220    @Override
221    public double getSupportUpperBound() {
222        return Double.POSITIVE_INFINITY;
223    }
224
225    /**
226     * {@inheritDoc}
227     *
228     * The support of this distribution is connected.
229     *
230     * @return {@code true}
231     */
232    @Override
233    public boolean isSupportConnected() {
234        return true;
235    }
236
237    /** {@inheritDoc} */
238    @Override
239    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
240        // Log normal distribution sampler.
241        return new LogNormalSampler(new ZigguratNormalizedGaussianSampler(rng), scale, shape)::sample;
242    }
243}