| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| package org.apache.commons.statistics.distribution; |
| |
| import org.apache.commons.numbers.gamma.Erf; |
| import org.apache.commons.numbers.gamma.ErfDifference; |
| import org.apache.commons.rng.UniformRandomProvider; |
| import org.apache.commons.rng.sampling.distribution.LogNormalSampler; |
| import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; |
| |
| /** |
| * Implementation of the <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">log-normal distribution</a>. |
| * |
| * <p> |
| * <strong>Parameters:</strong> |
| * {@code X} is log-normally distributed if its natural logarithm {@code log(X)} |
| * is normally distributed. The probability distribution function of {@code X} |
| * is given by (for {@code x > 0}) |
| * </p> |
| * <p> |
| * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)} |
| * </p> |
| * <ul> |
| * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the |
| * normally distributed natural logarithm of this distribution,</li> |
| * <li>{@code s} is the <em>shape</em> parameter: this is the standard |
| * deviation of the normally distributed natural logarithm of this |
| * distribution. |
| * </ul> |
| */ |
| public class LogNormalDistribution extends AbstractContinuousDistribution { |
| /** √(2 π). */ |
| private static final double SQRT2PI = Math.sqrt(2 * Math.PI); |
| /** √(2). */ |
| private static final double SQRT2 = Math.sqrt(2); |
| /** The scale parameter of this distribution. */ |
| private final double scale; |
| /** The shape parameter of this distribution. */ |
| private final double shape; |
| /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */ |
| private final double logShapePlusHalfLog2Pi; |
| |
| /** |
| * Creates a log-normal distribution. |
| * |
| * @param scale Scale parameter of this distribution. |
| * @param shape Shape parameter of this distribution. |
| * @throws IllegalArgumentException if {@code shape <= 0}. |
| */ |
| public LogNormalDistribution(double scale, |
| double shape) { |
| if (shape <= 0) { |
| throw new DistributionException(DistributionException.NEGATIVE, shape); |
| } |
| |
| this.scale = scale; |
| this.shape = shape; |
| this.logShapePlusHalfLog2Pi = Math.log(shape) + 0.5 * Math.log(2 * Math.PI); |
| } |
| |
| /** |
| * Returns the scale parameter of this distribution. |
| * |
| * @return the scale parameter |
| */ |
| public double getScale() { |
| return scale; |
| } |
| |
| /** |
| * Returns the shape parameter of this distribution. |
| * |
| * @return the shape parameter |
| */ |
| public double getShape() { |
| return shape; |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * For scale {@code m}, and shape {@code s} of this distribution, the PDF |
| * is given by |
| * <ul> |
| * <li>{@code 0} if {@code x <= 0},</li> |
| * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)} |
| * otherwise.</li> |
| * </ul> |
| */ |
| @Override |
| public double density(double x) { |
| if (x <= 0) { |
| return 0; |
| } |
| final double x0 = Math.log(x) - scale; |
| final double x1 = x0 / shape; |
| return Math.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x); |
| } |
| |
| /** {@inheritDoc} |
| * |
| * See documentation of {@link #density(double)} for computation details. |
| */ |
| @Override |
| public double logDensity(double x) { |
| if (x <= 0) { |
| return Double.NEGATIVE_INFINITY; |
| } |
| final double logX = Math.log(x); |
| final double x0 = logX - scale; |
| final double x1 = x0 / shape; |
| return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX); |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * For scale {@code m}, and shape {@code s} of this distribution, the CDF |
| * is given by |
| * <ul> |
| * <li>{@code 0} if {@code x <= 0},</li> |
| * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as |
| * in these cases the actual value is within {@code Double.MIN_VALUE} of 0, |
| * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s}, |
| * as in these cases the actual value is within {@code Double.MIN_VALUE} of |
| * 1,</li> |
| * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li> |
| * </ul> |
| */ |
| @Override |
| public double cumulativeProbability(double x) { |
| if (x <= 0) { |
| return 0; |
| } |
| final double dev = Math.log(x) - scale; |
| if (Math.abs(dev) > 40 * shape) { |
| return dev < 0 ? 0.0d : 1.0d; |
| } |
| return 0.5 + 0.5 * Erf.value(dev / (shape * SQRT2)); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public double probability(double x0, |
| double x1) { |
| if (x0 > x1) { |
| throw new DistributionException(DistributionException.TOO_LARGE, |
| x0, x1); |
| } |
| if (x0 <= 0 || x1 <= 0) { |
| return super.probability(x0, x1); |
| } |
| final double denom = shape * SQRT2; |
| final double v0 = (Math.log(x0) - scale) / denom; |
| final double v1 = (Math.log(x1) - scale) / denom; |
| return 0.5 * ErfDifference.value(v0, v1); |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * For scale {@code m} and shape {@code s}, the mean is |
| * {@code exp(m + s^2 / 2)}. |
| */ |
| @Override |
| public double getMean() { |
| final double s = shape; |
| return Math.exp(scale + (s * s / 2)); |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * For scale {@code m} and shape {@code s}, the variance is |
| * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}. |
| */ |
| @Override |
| public double getVariance() { |
| final double s = shape; |
| final double ss = s * s; |
| return (Math.expm1(ss)) * Math.exp(2 * scale + ss); |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * The lower bound of the support is always 0 no matter the parameters. |
| * |
| * @return lower bound of the support (always 0) |
| */ |
| @Override |
| public double getSupportLowerBound() { |
| return 0; |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * The upper bound of the support is always positive infinity |
| * no matter the parameters. |
| * |
| * @return upper bound of the support (always |
| * {@code Double.POSITIVE_INFINITY}) |
| */ |
| @Override |
| public double getSupportUpperBound() { |
| return Double.POSITIVE_INFINITY; |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * The support of this distribution is connected. |
| * |
| * @return {@code true} |
| */ |
| @Override |
| public boolean isSupportConnected() { |
| return true; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) { |
| // Log normal distribution sampler. |
| return new LogNormalSampler(new ZigguratNormalizedGaussianSampler(rng), scale, shape)::sample; |
| } |
| } |