    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
    018    package org.apache.commons.math3.distribution;
    019
    020    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    021    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    022    import org.apache.commons.math3.exception.util.LocalizedFormats;
    023    import org.apache.commons.math3.special.Erf;
    024    import org.apache.commons.math3.util.FastMath;
    025    import org.apache.commons.math3.random.RandomGenerator;
    026    import org.apache.commons.math3.random.Well19937c;
    027
    028    /**
    029     * Implementation of the log-normal (gaussian) distribution.
    030     *
    031     * <p>
    032     * <strong>Parameters:</strong>
    033     * {@code X} is log-normally distributed if its natural logarithm {@code log(X)}
    034     * is normally distributed. The probability distribution function of {@code X}
    035     * is given by (for {@code x > 0})
    036     * </p>
    037     * <p>
    038     * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
    039     * </p>
    040     * <ul>
    041     * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the
    042     * normally distributed natural logarithm of this distribution,</li>
    043     * <li>{@code s} is the <em>shape</em> parameter: this is the standard
    044     * deviation of the normally distributed natural logarithm of this
    045     * distribution.
    046     * </ul>
    047     *
    048     * @see <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">
    049     * Log-normal distribution (Wikipedia)</a>
    050     * @see <a href="http://mathworld.wolfram.com/LogNormalDistribution.html">
    051     * Log Normal distribution (MathWorld)</a>
    052     *
    053     * @version $Id: LogNormalDistribution.java 1422195 2012-12-15 06:45:18Z psteitz $
    054     * @since 3.0
    055     */
    056    public class LogNormalDistribution extends AbstractRealDistribution {
    057        /** Default inverse cumulative probability accuracy. */
    058        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    059
    060        /** Serializable version identifier. */
    061        private static final long serialVersionUID = 20120112;
    062
    063        /** &radic;(2 &pi;) */
    064        private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI);
    065
    066        /** &radic;(2) */
    067        private static final double SQRT2 = FastMath.sqrt(2.0);
    068
    069        /** The scale parameter of this distribution. */
    070        private final double scale;
    071
    072        /** The shape parameter of this distribution. */
    073        private final double shape;
    074
    075        /** Inverse cumulative probability accuracy. */
    076        private final double solverAbsoluteAccuracy;
    077
    078        /**
    079         * Create a log-normal distribution, where the mean and standard deviation
    080         * of the {@link NormalDistribution normally distributed} natural
    081         * logarithm of the log-normal distribution are equal to zero and one
    082         * respectively. In other words, the scale of the returned distribution is
    083         * {@code 0}, while its shape is {@code 1}.
    084         */
    085        public LogNormalDistribution() {
    086            this(0, 1);
    087        }
    088
    089        /**
    090         * Create a log-normal distribution using the specified scale and shape.
    091         *
    092         * @param scale the scale parameter of this distribution
    093         * @param shape the shape parameter of this distribution
    094         * @throws NotStrictlyPositiveException if {@code shape <= 0}.
    095         */
    096        public LogNormalDistribution(double scale, double shape)
    097            throws NotStrictlyPositiveException {
    098            this(scale, shape, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    099        }
    100
    101        /**
    102         * Create a log-normal distribution using the specified scale, shape and
    103         * inverse cumulative distribution accuracy.
    104         *
    105         * @param scale the scale parameter of this distribution
    106         * @param shape the shape parameter of this distribution
    107         * @param inverseCumAccuracy Inverse cumulative probability accuracy.
    108         * @throws NotStrictlyPositiveException if {@code shape <= 0}.
    109         */
    110        public LogNormalDistribution(double scale, double shape, double inverseCumAccuracy)
    111            throws NotStrictlyPositiveException {
    112            this(new Well19937c(), scale, shape, inverseCumAccuracy);
    113        }
    114
    115        /**
    116         * Creates a log-normal distribution.
    117         *
    118         * @param rng Random number generator.
    119         * @param scale Scale parameter of this distribution.
    120         * @param shape Shape parameter of this distribution.
    121         * @param inverseCumAccuracy Inverse cumulative probability accuracy.
    122         * @throws NotStrictlyPositiveException if {@code shape <= 0}.
    123         * @since 3.1
    124         */
    125        public LogNormalDistribution(RandomGenerator rng,
    126                                     double scale,
    127                                     double shape,
    128                                     double inverseCumAccuracy)
    129            throws NotStrictlyPositiveException {
    130            super(rng);
    131
    132            if (shape <= 0) {
    133                throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape);
    134            }
    135
    136            this.scale = scale;
    137            this.shape = shape;
    138            this.solverAbsoluteAccuracy = inverseCumAccuracy;
    139        }
    140
    141        /**
    142         * Returns the scale parameter of this distribution.
    143         *
    144         * @return the scale parameter
    145         */
    146        public double getScale() {
    147            return scale;
    148        }
    149
    150        /**
    151         * Returns the shape parameter of this distribution.
    152         *
    153         * @return the shape parameter
    154         */
    155        public double getShape() {
    156            return shape;
    157        }
    158
    159        /**
    160         * {@inheritDoc}
    161         *
    162         * For scale {@code m}, and shape {@code s} of this distribution, the PDF
    163         * is given by
    164         * <ul>
    165         * <li>{@code 0} if {@code x <= 0},</li>
    166         * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
    167         * otherwise.</li>
    168         * </ul>
    169         */
    170        public double density(double x) {
    171            if (x <= 0) {
    172                return 0;
    173            }
    174            final double x0 = FastMath.log(x) - scale;
    175            final double x1 = x0 / shape;
    176            return FastMath.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
    177        }
    178
    179        /**
    180         * {@inheritDoc}
    181         *
    182         * For scale {@code m}, and shape {@code s} of this distribution, the CDF
    183         * is given by
    184         * <ul>
    185         * <li>{@code 0} if {@code x <= 0},</li>
    186         * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as
    187         * in these cases the actual value is within {@code Double.MIN_VALUE} of 0,
    188         * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s},
    189         * as in these cases the actual value is within {@code Double.MIN_VALUE} of
    190         * 1,</li>
    191         * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li>
    192         * </ul>
    193         */
    194        public double cumulativeProbability(double x)  {
    195            if (x <= 0) {
    196                return 0;
    197            }
    198            final double dev = FastMath.log(x) - scale;
    199            if (FastMath.abs(dev) > 40 * shape) {
    200                return dev < 0 ? 0.0d : 1.0d;
    201            }
    202            return 0.5 + 0.5 * Erf.erf(dev / (shape * SQRT2));
    203        }
    204
    205        /**
    206         * {@inheritDoc}
    207         *
    208         * @deprecated See {@link RealDistribution#cumulativeProbability(double,double)}
    209         */
    210        @Override@Deprecated
    211        public double cumulativeProbability(double x0, double x1)
    212            throws NumberIsTooLargeException {
    213            return probability(x0, x1);
    214        }
    215
    216        /** {@inheritDoc} */
    217        @Override
    218        public double probability(double x0,
    219                                  double x1)
    220            throws NumberIsTooLargeException {
    221            if (x0 > x1) {
    222                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
    223                                                    x0, x1, true);
    224            }
    225            if (x0 <= 0 || x1 <= 0) {
    226                return super.probability(x0, x1);
    227            }
    228            final double denom = shape * SQRT2;
    229            final double v0 = (FastMath.log(x0) - scale) / denom;
    230            final double v1 = (FastMath.log(x1) - scale) / denom;
    231            return 0.5 * Erf.erf(v0, v1);
    232        }
    233
    234        /** {@inheritDoc} */
    235        @Override
    236        protected double getSolverAbsoluteAccuracy() {
    237            return solverAbsoluteAccuracy;
    238        }
    239
    240        /**
    241         * {@inheritDoc}
    242         *
    243         * For scale {@code m} and shape {@code s}, the mean is
    244         * {@code exp(m + s^2 / 2)}.
    245         */
    246        public double getNumericalMean() {
    247            double s = shape;
    248            return FastMath.exp(scale + (s * s / 2));
    249        }
    250
    251        /**
    252         * {@inheritDoc}
    253         *
    254         * For scale {@code m} and shape {@code s}, the variance is
    255         * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}.
    256         */
    257        public double getNumericalVariance() {
    258            final double s = shape;
    259            final double ss = s * s;
    260            return (FastMath.exp(ss) - 1) * FastMath.exp(2 * scale + ss);
    261        }
    262
    263        /**
    264         * {@inheritDoc}
    265         *
    266         * The lower bound of the support is always 0 no matter the parameters.
    267         *
    268         * @return lower bound of the support (always 0)
    269         */
    270        public double getSupportLowerBound() {
    271            return 0;
    272        }
    273
    274        /**
    275         * {@inheritDoc}
    276         *
    277         * The upper bound of the support is always positive infinity
    278         * no matter the parameters.
    279         *
    280         * @return upper bound of the support (always
    281         * {@code Double.POSITIVE_INFINITY})
    282         */
    283        public double getSupportUpperBound() {
    284            return Double.POSITIVE_INFINITY;
    285        }
    286
    287        /** {@inheritDoc} */
    288        public boolean isSupportLowerBoundInclusive() {
    289            return true;
    290        }
    291
    292        /** {@inheritDoc} */
    293        public boolean isSupportUpperBoundInclusive() {
    294            return false;
    295        }
    296
    297        /**
    298         * {@inheritDoc}
    299         *
    300         * The support of this distribution is connected.
    301         *
    302         * @return {@code true}
    303         */
    304        public boolean isSupportConnected() {
    305            return true;
    306        }
    307
    308        /** {@inheritDoc} */
    309        @Override
    310        public double sample()  {
    311            final double n = random.nextGaussian();
    312            return FastMath.exp(scale + shape * n);
    313        }
    314    }