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 /** √(2 π) */ 064 private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI); 065 066 /** √(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 }