source: src/main/java/agents/anac/y2019/harddealer/math3/distribution/LogNormalDistribution.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 12.4 KB
Line 
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 * http://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
18package agents.anac.y2019.harddealer.math3.distribution;
19
20import agents.anac.y2019.harddealer.math3.exception.NotStrictlyPositiveException;
21import agents.anac.y2019.harddealer.math3.exception.NumberIsTooLargeException;
22import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
23import agents.anac.y2019.harddealer.math3.random.RandomGenerator;
24import agents.anac.y2019.harddealer.math3.random.Well19937c;
25import agents.anac.y2019.harddealer.math3.special.Erf;
26import agents.anac.y2019.harddealer.math3.util.FastMath;
27
28/**
29 * Implementation of the log-normal (gaussian) distribution.
30 *
31 * <p>
32 * <strong>Parameters:</strong>
33 * {@code X} is log-normally distributed if its natural logarithm {@code log(X)}
34 * is normally distributed. The probability distribution function of {@code X}
35 * is given by (for {@code x > 0})
36 * </p>
37 * <p>
38 * {@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
39 * </p>
40 * <ul>
41 * <li>{@code m} is the <em>scale</em> parameter: this is the mean of the
42 * normally distributed natural logarithm of this distribution,</li>
43 * <li>{@code s} is the <em>shape</em> parameter: this is the standard
44 * deviation of the normally distributed natural logarithm of this
45 * distribution.
46 * </ul>
47 *
48 * @see <a href="http://en.wikipedia.org/wiki/Log-normal_distribution">
49 * Log-normal distribution (Wikipedia)</a>
50 * @see <a href="http://mathworld.wolfram.com/LogNormalDistribution.html">
51 * Log Normal distribution (MathWorld)</a>
52 *
53 * @since 3.0
54 */
55public class LogNormalDistribution extends AbstractRealDistribution {
56 /** Default inverse cumulative probability accuracy. */
57 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
58
59 /** Serializable version identifier. */
60 private static final long serialVersionUID = 20120112;
61
62 /** &radic;(2 &pi;) */
63 private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI);
64
65 /** &radic;(2) */
66 private static final double SQRT2 = FastMath.sqrt(2.0);
67
68 /** The scale parameter of this distribution. */
69 private final double scale;
70
71 /** The shape parameter of this distribution. */
72 private final double shape;
73 /** The value of {@code log(shape) + 0.5 * log(2*PI)} stored for faster computation. */
74 private final double logShapePlusHalfLog2Pi;
75
76 /** Inverse cumulative probability accuracy. */
77 private final double solverAbsoluteAccuracy;
78
79 /**
80 * Create a log-normal distribution, where the mean and standard deviation
81 * of the {@link NormalDistribution normally distributed} natural
82 * logarithm of the log-normal distribution are equal to zero and one
83 * respectively. In other words, the scale of the returned distribution is
84 * {@code 0}, while its shape is {@code 1}.
85 * <p>
86 * <b>Note:</b> this constructor will implicitly create an instance of
87 * {@link Well19937c} as random generator to be used for sampling only (see
88 * {@link #sample()} and {@link #sample(int)}). In case no sampling is
89 * needed for the created distribution, it is advised to pass {@code null}
90 * as random generator via the appropriate constructors to avoid the
91 * additional initialisation overhead.
92 */
93 public LogNormalDistribution() {
94 this(0, 1);
95 }
96
97 /**
98 * Create a log-normal distribution using the specified scale and shape.
99 * <p>
100 * <b>Note:</b> this constructor will implicitly create an instance of
101 * {@link Well19937c} as random generator to be used for sampling only (see
102 * {@link #sample()} and {@link #sample(int)}). In case no sampling is
103 * needed for the created distribution, it is advised to pass {@code null}
104 * as random generator via the appropriate constructors to avoid the
105 * additional initialisation overhead.
106 *
107 * @param scale the scale parameter of this distribution
108 * @param shape the shape parameter of this distribution
109 * @throws NotStrictlyPositiveException if {@code shape <= 0}.
110 */
111 public LogNormalDistribution(double scale, double shape)
112 throws NotStrictlyPositiveException {
113 this(scale, shape, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
114 }
115
116 /**
117 * Create a log-normal distribution using the specified scale, shape and
118 * inverse cumulative distribution accuracy.
119 * <p>
120 * <b>Note:</b> this constructor will implicitly create an instance of
121 * {@link Well19937c} as random generator to be used for sampling only (see
122 * {@link #sample()} and {@link #sample(int)}). In case no sampling is
123 * needed for the created distribution, it is advised to pass {@code null}
124 * as random generator via the appropriate constructors to avoid the
125 * additional initialisation overhead.
126 *
127 * @param scale the scale parameter of this distribution
128 * @param shape the shape parameter of this distribution
129 * @param inverseCumAccuracy Inverse cumulative probability accuracy.
130 * @throws NotStrictlyPositiveException if {@code shape <= 0}.
131 */
132 public LogNormalDistribution(double scale, double shape, double inverseCumAccuracy)
133 throws NotStrictlyPositiveException {
134 this(new Well19937c(), scale, shape, inverseCumAccuracy);
135 }
136
137 /**
138 * Creates a log-normal distribution.
139 *
140 * @param rng Random number generator.
141 * @param scale Scale parameter of this distribution.
142 * @param shape Shape parameter of this distribution.
143 * @throws NotStrictlyPositiveException if {@code shape <= 0}.
144 * @since 3.3
145 */
146 public LogNormalDistribution(RandomGenerator rng, double scale, double shape)
147 throws NotStrictlyPositiveException {
148 this(rng, scale, shape, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
149 }
150
151 /**
152 * Creates a log-normal distribution.
153 *
154 * @param rng Random number generator.
155 * @param scale Scale parameter of this distribution.
156 * @param shape Shape parameter of this distribution.
157 * @param inverseCumAccuracy Inverse cumulative probability accuracy.
158 * @throws NotStrictlyPositiveException if {@code shape <= 0}.
159 * @since 3.1
160 */
161 public LogNormalDistribution(RandomGenerator rng,
162 double scale,
163 double shape,
164 double inverseCumAccuracy)
165 throws NotStrictlyPositiveException {
166 super(rng);
167
168 if (shape <= 0) {
169 throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape);
170 }
171
172 this.scale = scale;
173 this.shape = shape;
174 this.logShapePlusHalfLog2Pi = FastMath.log(shape) + 0.5 * FastMath.log(2 * FastMath.PI);
175 this.solverAbsoluteAccuracy = inverseCumAccuracy;
176 }
177
178 /**
179 * Returns the scale parameter of this distribution.
180 *
181 * @return the scale parameter
182 */
183 public double getScale() {
184 return scale;
185 }
186
187 /**
188 * Returns the shape parameter of this distribution.
189 *
190 * @return the shape parameter
191 */
192 public double getShape() {
193 return shape;
194 }
195
196 /**
197 * {@inheritDoc}
198 *
199 * For scale {@code m}, and shape {@code s} of this distribution, the PDF
200 * is given by
201 * <ul>
202 * <li>{@code 0} if {@code x <= 0},</li>
203 * <li>{@code exp(-0.5 * ((ln(x) - m) / s)^2) / (s * sqrt(2 * pi) * x)}
204 * otherwise.</li>
205 * </ul>
206 */
207 public double density(double x) {
208 if (x <= 0) {
209 return 0;
210 }
211 final double x0 = FastMath.log(x) - scale;
212 final double x1 = x0 / shape;
213 return FastMath.exp(-0.5 * x1 * x1) / (shape * SQRT2PI * x);
214 }
215
216 /** {@inheritDoc}
217 *
218 * See documentation of {@link #density(double)} for computation details.
219 */
220 @Override
221 public double logDensity(double x) {
222 if (x <= 0) {
223 return Double.NEGATIVE_INFINITY;
224 }
225 final double logX = FastMath.log(x);
226 final double x0 = logX - scale;
227 final double x1 = x0 / shape;
228 return -0.5 * x1 * x1 - (logShapePlusHalfLog2Pi + logX);
229 }
230
231 /**
232 * {@inheritDoc}
233 *
234 * For scale {@code m}, and shape {@code s} of this distribution, the CDF
235 * is given by
236 * <ul>
237 * <li>{@code 0} if {@code x <= 0},</li>
238 * <li>{@code 0} if {@code ln(x) - m < 0} and {@code m - ln(x) > 40 * s}, as
239 * in these cases the actual value is within {@code Double.MIN_VALUE} of 0,
240 * <li>{@code 1} if {@code ln(x) - m >= 0} and {@code ln(x) - m > 40 * s},
241 * as in these cases the actual value is within {@code Double.MIN_VALUE} of
242 * 1,</li>
243 * <li>{@code 0.5 + 0.5 * erf((ln(x) - m) / (s * sqrt(2))} otherwise.</li>
244 * </ul>
245 */
246 public double cumulativeProbability(double x) {
247 if (x <= 0) {
248 return 0;
249 }
250 final double dev = FastMath.log(x) - scale;
251 if (FastMath.abs(dev) > 40 * shape) {
252 return dev < 0 ? 0.0d : 1.0d;
253 }
254 return 0.5 + 0.5 * Erf.erf(dev / (shape * SQRT2));
255 }
256
257 /**
258 * {@inheritDoc}
259 *
260 * @deprecated See {@link RealDistribution#cumulativeProbability(double,double)}
261 */
262 @Override@Deprecated
263 public double cumulativeProbability(double x0, double x1)
264 throws NumberIsTooLargeException {
265 return probability(x0, x1);
266 }
267
268 /** {@inheritDoc} */
269 @Override
270 public double probability(double x0,
271 double x1)
272 throws NumberIsTooLargeException {
273 if (x0 > x1) {
274 throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
275 x0, x1, true);
276 }
277 if (x0 <= 0 || x1 <= 0) {
278 return super.probability(x0, x1);
279 }
280 final double denom = shape * SQRT2;
281 final double v0 = (FastMath.log(x0) - scale) / denom;
282 final double v1 = (FastMath.log(x1) - scale) / denom;
283 return 0.5 * Erf.erf(v0, v1);
284 }
285
286 /** {@inheritDoc} */
287 @Override
288 protected double getSolverAbsoluteAccuracy() {
289 return solverAbsoluteAccuracy;
290 }
291
292 /**
293 * {@inheritDoc}
294 *
295 * For scale {@code m} and shape {@code s}, the mean is
296 * {@code exp(m + s^2 / 2)}.
297 */
298 public double getNumericalMean() {
299 double s = shape;
300 return FastMath.exp(scale + (s * s / 2));
301 }
302
303 /**
304 * {@inheritDoc}
305 *
306 * For scale {@code m} and shape {@code s}, the variance is
307 * {@code (exp(s^2) - 1) * exp(2 * m + s^2)}.
308 */
309 public double getNumericalVariance() {
310 final double s = shape;
311 final double ss = s * s;
312 return (FastMath.expm1(ss)) * FastMath.exp(2 * scale + ss);
313 }
314
315 /**
316 * {@inheritDoc}
317 *
318 * The lower bound of the support is always 0 no matter the parameters.
319 *
320 * @return lower bound of the support (always 0)
321 */
322 public double getSupportLowerBound() {
323 return 0;
324 }
325
326 /**
327 * {@inheritDoc}
328 *
329 * The upper bound of the support is always positive infinity
330 * no matter the parameters.
331 *
332 * @return upper bound of the support (always
333 * {@code Double.POSITIVE_INFINITY})
334 */
335 public double getSupportUpperBound() {
336 return Double.POSITIVE_INFINITY;
337 }
338
339 /** {@inheritDoc} */
340 public boolean isSupportLowerBoundInclusive() {
341 return true;
342 }
343
344 /** {@inheritDoc} */
345 public boolean isSupportUpperBoundInclusive() {
346 return false;
347 }
348
349 /**
350 * {@inheritDoc}
351 *
352 * The support of this distribution is connected.
353 *
354 * @return {@code true}
355 */
356 public boolean isSupportConnected() {
357 return true;
358 }
359
360 /** {@inheritDoc} */
361 @Override
362 public double sample() {
363 final double n = random.nextGaussian();
364 return FastMath.exp(scale + shape * n);
365 }
366}
Note: See TracBrowser for help on using the repository browser.