View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  
23  package org.hipparchus.distribution.continuous;
24  
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathIllegalArgumentException;
27  import org.hipparchus.special.Erf;
28  import org.hipparchus.util.FastMath;
29  
30  /**
31   * Implementation of the log-normal (gaussian) distribution.
32   * <p>
33   * <strong>Parameters:</strong>
34   * {@code X} is log-normally distributed if its natural logarithm {@code log(X)}
35   * is normally distributed. The probability distribution function of {@code X}
36   * is given by (for {@code x > 0})
37   * <p>
38   * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
39   * <ul>
40   * <li>{@code m} is the <em>location</em> parameter: this is the mean of the
41   * normally distributed natural logarithm of this distribution,</li>
42   * <li>{@code s} is the <em>shape</em> parameter: this is the standard
43   * deviation of the normally distributed natural logarithm of this
44   * distribution.
45   * </ul>
46   *
47   * @see <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">
48   * Log-normal distribution (Wikipedia)</a>
49   * @see <a href="http://mathworld.wolfram.com/LogNormalDistribution.html">
50   * Log Normal distribution (MathWorld)</a>
51   */
52  public class LogNormalDistribution extends AbstractRealDistribution {
53  
54      /** Serializable version identifier. */
55      private static final long serialVersionUID = 20120112;
56  
57      /** &radic;(2 &pi;) */
58      private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI);
59  
60      /** &radic;(2) */
61      private static final double SQRT2 = FastMath.sqrt(2.0);
62  
63      /** The location parameter of this distribution (named m in MathWorld and ยต in Wikipedia). */
64      private final double location;
65  
66      /** The shape parameter of this distribution. */
67      private final double shape;
68      /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */
69      private final double logShapePlusHalfLog2Pi;
70  
71      /**
72       * Create a log-normal distribution, where the mean and standard deviation
73       * of the {@link NormalDistribution normally distributed} natural
74       * logarithm of the log-normal distribution are equal to zero and one
75       * respectively. In other words, the location of the returned distribution is
76       * {@code 0}, while its shape is {@code 1}.
77       */
78      public LogNormalDistribution() {
79          this(0, 1);
80      }
81  
82      /**
83       * Create a log-normal distribution using the specified location and shape.
84       *
85       * @param location the location parameter of this distribution
86       * @param shape the shape parameter of this distribution
87       * @throws MathIllegalArgumentException if {@code shape <= 0}.
88       */
89      public LogNormalDistribution(double location, double shape)
90          throws MathIllegalArgumentException {
91          this(location, shape, DEFAULT_SOLVER_ABSOLUTE_ACCURACY);
92      }
93  
94  
95      /**
96       * Creates a log-normal distribution.
97       *
98       * @param location Location parameter of this distribution.
99       * @param shape Shape parameter of this distribution.
100      * @param inverseCumAccuracy Inverse cumulative probability accuracy.
101      * @throws MathIllegalArgumentException if {@code shape <= 0}.
102      */
103     public LogNormalDistribution(double location,
104                                  double shape,
105                                  double inverseCumAccuracy)
106         throws MathIllegalArgumentException {
107         super(inverseCumAccuracy);
108 
109         if (shape <= 0) {
110             throw new MathIllegalArgumentException(LocalizedCoreFormats.SHAPE, shape);
111         }
112 
113         this.location = location;
114         this.shape = shape;
115         this.logShapePlusHalfLog2Pi = FastMath.log(shape) + 0.5 * FastMath.log(2 * FastMath.PI);
116     }
117 
118     /**
119      * Returns the location parameter of this distribution.
120      *
121      * @return the location parameter
122      * @since 1.4
123      */
124     public double getLocation() {
125         return location;
126     }
127 
128     /**
129      * Returns the shape parameter of this distribution.
130      *
131      * @return the shape parameter
132      */
133     public double getShape() {
134         return shape;
135     }
136 
137     /**
138      * {@inheritDoc}
139      *
140      * For location {@code m}, and shape {@code s} of this distribution, the PDF
141      * is given by
142      * <ul>
143      * <li>{@code 0} if {@code x <= 0},</li>
144      * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
145      * otherwise.</li>
146      * </ul>
147      */
148     @Override
149     public double density(double x) {
150         if (x <= 0) {
151             return 0;
152         }
153         final double x0 = FastMath.log(x) - location;
154         final double x1 = x0 / shape;
155         return FastMath.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
156     }
157 
158     /** {@inheritDoc}
159      *
160      * See documentation of {@link #density(double)} for computation details.
161      */
162     @Override
163     public double logDensity(double x) {
164         if (x <= 0) {
165             return Double.NEGATIVE_INFINITY;
166         }
167         final double logX = FastMath.log(x);
168         final double x0 = logX - location;
169         final double x1 = x0 / shape;
170         return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX);
171     }
172 
173     /**
174      * {@inheritDoc}
175      *
176      * For location {@code m}, and shape {@code s} of this distribution, the CDF
177      * is given by
178      * <ul>
179      * <li>{@code 0} if {@code x <= 0},</li>
180      * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as
181      * in these cases the actual value is within {@code Double.MIN_VALUE} of 0,
182      * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s},
183      * as in these cases the actual value is within {@code Double.MIN_VALUE} of 1,</li>
184      * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li>
185      * </ul>
186      */
187     @Override
188     public double cumulativeProbability(double x)  {
189         if (x <= 0) {
190             return 0;
191         }
192         final double dev = FastMath.log(x) - location;
193         if (FastMath.abs(dev) > 40 * shape) {
194             return dev < 0 ? 0.0d : 1.0d;
195         }
196         return 0.5 + 0.5 * Erf.erf(dev / (shape * SQRT2));
197     }
198 
199     /** {@inheritDoc} */
200     @Override
201     public double probability(double x0,
202                               double x1)
203         throws MathIllegalArgumentException {
204         if (x0 > x1) {
205             throw new MathIllegalArgumentException(LocalizedCoreFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
206                                                 x0, x1, true);
207         }
208         if (x0 <= 0 || x1 <= 0) {
209             return super.probability(x0, x1);
210         }
211         final double denom = shape * SQRT2;
212         final double v0 = (FastMath.log(x0) - location) / denom;
213         final double v1 = (FastMath.log(x1) - location) / denom;
214         return 0.5 * Erf.erf(v0, v1);
215     }
216 
217     /**
218      * {@inheritDoc}
219      *
220      * For location {@code m} and shape {@code s}, the mean is
221      * {@code exp(m + s^2 / 2)}.
222      */
223     @Override
224     public double getNumericalMean() {
225         double s = shape;
226         return FastMath.exp(location + (s * s / 2));
227     }
228 
229     /**
230      * {@inheritDoc}
231      *
232      * For location {@code m} and shape {@code s}, the variance is
233      * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}.
234      */
235     @Override
236     public double getNumericalVariance() {
237         final double s = shape;
238         final double ss = s * s;
239         return (FastMath.expm1(ss)) * FastMath.exp(2 * location + ss);
240     }
241 
242     /**
243      * {@inheritDoc}
244      *
245      * The lower bound of the support is always 0 no matter the parameters.
246      *
247      * @return lower bound of the support (always 0)
248      */
249     @Override
250     public double getSupportLowerBound() {
251         return 0;
252     }
253 
254     /**
255      * {@inheritDoc}
256      *
257      * The upper bound of the support is always positive infinity
258      * no matter the parameters.
259      *
260      * @return upper bound of the support (always
261      * {@code Double.POSITIVE_INFINITY})
262      */
263     @Override
264     public double getSupportUpperBound() {
265         return Double.POSITIVE_INFINITY;
266     }
267 
268     /**
269      * {@inheritDoc}
270      *
271      * The support of this distribution is connected.
272      *
273      * @return {@code true}
274      */
275     @Override
276     public boolean isSupportConnected() {
277         return true;
278     }
279 }