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    }