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.ErfDifference;
021import org.apache.commons.numbers.gamma.Erfc;
022import org.apache.commons.numbers.gamma.InverseErfc;
023import org.apache.commons.rng.UniformRandomProvider;
024import org.apache.commons.rng.sampling.distribution.LogNormalSampler;
025import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
026
027/**
028 * Implementation of the log-normal distribution.
029 *
030 * <p>\( X \) is log-normally distributed if its natural logarithm \( \ln(x) \)
031 * is normally distributed. The probability density function of \( X \) is:
032 *
033 * <p>\[ f(x; \mu, \sigma) = \frac 1 {x\sigma\sqrt{2\pi\,}} e^{-{\frac 1 2}\left( \frac{\ln x-\mu}{\sigma} \right)^2 } \]
034 *
035 * <p>for \( \mu \) the mean of the normally distributed natural logarithm of this distribution,
036 * \( \sigma &gt; 0 \) the standard deviation of the normally distributed natural logarithm of this
037 * distribution, and
038 * \( x \in (0, \infty) \).
039 *
040 * @see <a href="https://en.wikipedia.org/wiki/Log-normal_distribution">Log-normal distribution (Wikipedia)</a>
041 * @see <a href="https://mathworld.wolfram.com/LogNormalDistribution.html">Log-normal distribution (MathWorld)</a>
042 */
043public final class LogNormalDistribution extends AbstractContinuousDistribution {
044    /** &radic;(2 &pi;). */
045    private static final double SQRT2PI = Math.sqrt(2 * Math.PI);
046    /** The mu parameter of this distribution. */
047    private final double mu;
048    /** The sigma parameter of this distribution. */
049    private final double sigma;
050    /** The value of {@code log(sigma) + 0.5 * log(2*PI)} stored for faster computation. */
051    private final double logSigmaPlusHalfLog2Pi;
052    /** Sigma multiplied by sqrt(2). */
053    private final double sigmaSqrt2;
054    /** Sigma multiplied by sqrt(2 * pi). */
055    private final double sigmaSqrt2Pi;
056
057    /**
058     * @param mu Mean of the natural logarithm of the distribution values.
059     * @param sigma Standard deviation of the natural logarithm of the distribution values.
060     */
061    private LogNormalDistribution(double mu,
062                                  double sigma) {
063        this.mu = mu;
064        this.sigma = sigma;
065        logSigmaPlusHalfLog2Pi = Math.log(sigma) + Constants.HALF_LOG_TWO_PI;
066        sigmaSqrt2 = ExtendedPrecision.sqrt2xx(sigma);
067        sigmaSqrt2Pi = sigma * SQRT2PI;
068    }
069
070    /**
071     * Creates a log-normal distribution.
072     *
073     * @param mu Mean of the natural logarithm of the distribution values.
074     * @param sigma Standard deviation of the natural logarithm of the distribution values.
075     * @return the distribution
076     * @throws IllegalArgumentException if {@code sigma <= 0}.
077     */
078    public static LogNormalDistribution of(double mu,
079                                           double sigma) {
080        if (sigma <= 0) {
081            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, sigma);
082        }
083        return new LogNormalDistribution(mu, sigma);
084    }
085
086    /**
087     * Gets the {@code mu} parameter of this distribution.
088     * This is the mean of the natural logarithm of the distribution values,
089     * not the mean of distribution.
090     *
091     * @return the mu parameter.
092     */
093    public double getMu() {
094        return mu;
095    }
096
097    /**
098     * Gets the {@code sigma} parameter of this distribution.
099     * This is the standard deviation of the natural logarithm of the distribution values,
100     * not the standard deviation of distribution.
101     *
102     * @return the sigma parameter.
103     */
104    public double getSigma() {
105        return sigma;
106    }
107
108    /**
109     * {@inheritDoc}
110     *
111     * <p>For {@code mu}, and sigma {@code s} of this distribution, the PDF
112     * is given by
113     * <ul>
114     * <li>{@code 0} if {@code x <= 0},</li>
115     * <li>{@code exp(-0.5 * ((ln(x) - mu) / s)^2) / (s * sqrt(2 * pi) * x)}
116     * otherwise.</li>
117     * </ul>
118     */
119    @Override
120    public double density(double x) {
121        if (x <= 0) {
122            return 0;
123        }
124        final double x0 = Math.log(x) - mu;
125        final double x1 = x0 / sigma;
126        return Math.exp(-0.5 * x1 * x1) / (sigmaSqrt2Pi * x);
127    }
128
129    /** {@inheritDoc} */
130    @Override
131    public double probability(double x0,
132                              double x1) {
133        if (x0 > x1) {
134            throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GT_HIGH,
135                                            x0, x1);
136        }
137        if (x0 <= 0) {
138            return cumulativeProbability(x1);
139        }
140        // Assumes x1 >= x0 && x0 > 0
141        final double v0 = (Math.log(x0) - mu) / sigmaSqrt2;
142        final double v1 = (Math.log(x1) - mu) / sigmaSqrt2;
143        return 0.5 * ErfDifference.value(v0, v1);
144    }
145
146    /** {@inheritDoc}
147     *
148     * <p>See documentation of {@link #density(double)} for computation details.
149     */
150    @Override
151    public double logDensity(double x) {
152        if (x <= 0) {
153            return Double.NEGATIVE_INFINITY;
154        }
155        final double logX = Math.log(x);
156        final double x0 = logX - mu;
157        final double x1 = x0 / sigma;
158        return -0.5 * x1 * x1 - (logSigmaPlusHalfLog2Pi + logX);
159    }
160
161    /**
162     * {@inheritDoc}
163     *
164     * <p>For {@code mu}, and sigma {@code s} of this distribution, the CDF
165     * is given by
166     * <ul>
167     * <li>{@code 0} if {@code x <= 0},</li>
168     * <li>{@code 0} if {@code ln(x) - mu < 0} and {@code mu - ln(x) > 40 * s}, as
169     * in these cases the actual value is within {@link Double#MIN_VALUE} of 0,
170     * <li>{@code 1} if {@code ln(x) - mu >= 0} and {@code ln(x) - mu > 40 * s},
171     * as in these cases the actual value is within {@link Double#MIN_VALUE} of
172     * 1,</li>
173     * <li>{@code 0.5 + 0.5 * erf((ln(x) - mu) / (s * sqrt(2))} otherwise.</li>
174     * </ul>
175     */
176    @Override
177    public double cumulativeProbability(double x)  {
178        if (x <= 0) {
179            return 0;
180        }
181        final double dev = Math.log(x) - mu;
182        return 0.5 * Erfc.value(-dev / sigmaSqrt2);
183    }
184
185    /** {@inheritDoc} */
186    @Override
187    public double survivalProbability(double x)  {
188        if (x <= 0) {
189            return 1;
190        }
191        final double dev = Math.log(x) - mu;
192        return 0.5 * Erfc.value(dev / sigmaSqrt2);
193    }
194
195    /** {@inheritDoc} */
196    @Override
197    public double inverseCumulativeProbability(double p) {
198        ArgumentUtils.checkProbability(p);
199        return Math.exp(mu - sigmaSqrt2 * InverseErfc.value(2 * p));
200    }
201
202    /** {@inheritDoc} */
203    @Override
204    public double inverseSurvivalProbability(double p) {
205        ArgumentUtils.checkProbability(p);
206        return Math.exp(mu + sigmaSqrt2 * InverseErfc.value(2 * p));
207    }
208
209    /**
210     * {@inheritDoc}
211     *
212     * <p>For \( \mu \) the mean of the normally distributed natural logarithm of
213     * this distribution, \( \sigma &gt; 0 \) the standard deviation of the normally
214     * distributed natural logarithm of this distribution, the mean is:
215     *
216     * <p>\[ \exp(\mu + \frac{\sigma^2}{2}) \]
217     */
218    @Override
219    public double getMean() {
220        final double s = sigma;
221        return Math.exp(mu + (s * s / 2));
222    }
223
224    /**
225     * {@inheritDoc}
226     *
227     * <p>For \( \mu \) the mean of the normally distributed natural logarithm of
228     * this distribution, \( \sigma &gt; 0 \) the standard deviation of the normally
229     * distributed natural logarithm of this distribution, the variance is:
230     *
231     * <p>\[ [\exp(\sigma^2) - 1)] \exp(2 \mu + \sigma^2) \]
232     */
233    @Override
234    public double getVariance() {
235        final double s = sigma;
236        final double ss = s * s;
237        return Math.expm1(ss) * Math.exp(2 * mu + ss);
238    }
239
240    /**
241     * {@inheritDoc}
242     *
243     * <p>The lower bound of the support is always 0.
244     *
245     * @return 0.
246     */
247    @Override
248    public double getSupportLowerBound() {
249        return 0;
250    }
251
252    /**
253     * {@inheritDoc}
254     *
255     * <p>The upper bound of the support is always positive infinity.
256     *
257     * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
258     */
259    @Override
260    public double getSupportUpperBound() {
261        return Double.POSITIVE_INFINITY;
262    }
263
264    /** {@inheritDoc} */
265    @Override
266    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
267        // Log normal distribution sampler.
268        final ZigguratSampler.NormalizedGaussian gaussian = ZigguratSampler.NormalizedGaussian.of(rng);
269        return LogNormalSampler.of(gaussian, mu, sigma)::sample;
270    }
271}