1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.statistics;
17
18 /**
19 * Contains methods to work with normal (i.e. Gaussian) distributions.
20 * Methods of this class use methods of the Erf class.
21 * This class is based in code of Numerical Recipes 3rd ed. section 6.14.1.
22 */
23 public class NormalDist {
24 /**
25 * Square root of 2.
26 */
27 private static final double SQRT2 = Math.sqrt(2.0);
28
29 /**
30 * Square root of 2 divided by 2.
31 */
32 private static final double HALF_SQRT2 = SQRT2 / 2.0;
33
34 /**
35 * Term to normalize Gaussian so that its integral from -infinity to infinity is one.
36 */
37 private static final double GAUSSIAN_NORM = 1.0 / Math.sqrt(2.0 * Math.PI);
38
39 /**
40 * Mean value of Gaussian distribution.
41 */
42 private double mu;
43
44 /**
45 * Standard deviation of Gaussian distribution.
46 */
47 private double sig;
48
49 /**
50 * Constructor. Initializes a Gaussian distribution with zero mean and
51 * unitary standard deviation (i.e. N(0,1)).
52 */
53 public NormalDist() {
54 mu = 0.0;
55 sig = 1.0;
56 }
57
58 /**
59 * Constructor with mean and standard deviation.
60 *
61 * @param mu mean value of Gaussian distribution.
62 * @param sig standard deviation of Gaussian distribution.
63 * @throws IllegalArgumentException if provided standard deviation is zero
64 * or negative.
65 */
66 public NormalDist(final double mu, final double sig) {
67 setStandardDeviation(sig);
68 setMean(mu);
69 }
70
71 /**
72 * Gets mean value of Gaussian distribution.
73 *
74 * @return mean value of Gaussian distribution.
75 */
76 public double getMean() {
77 return mu;
78 }
79
80 /**
81 * Sets mean value of Gaussian distribution.
82 *
83 * @param mu mean value of Gaussian distribution.
84 */
85 public final void setMean(final double mu) {
86 this.mu = mu;
87 }
88
89 /**
90 * Gets standard deviation of Gaussian distribution.
91 *
92 * @return standard deviation of Gaussian distribution.
93 */
94 public double getStandardDeviation() {
95 return sig;
96 }
97
98 /**
99 * Sets standard deviation of Gaussian distribution.
100 *
101 * @param sig standard deviation to be set.
102 * @throws IllegalArgumentException if provided standard deviation is zero
103 * or negative.
104 */
105 public final void setStandardDeviation(final double sig) {
106 if (sig <= 0.0) {
107 throw new IllegalArgumentException();
108 }
109 this.sig = sig;
110 }
111
112 /**
113 * Gets variance of Gaussian distribution.
114 *
115 * @return variance of Gaussian distribution.
116 */
117 public double getVariance() {
118 return sig * sig;
119 }
120
121 /**
122 * Sets variance of Gaussian distribution.
123 *
124 * @param variance variance of Gaussian distribution.
125 * @throws IllegalArgumentException if provided variance is zero or
126 * negative.
127 */
128 public void setVariance(final double variance) {
129 if (variance <= 0.0) {
130 throw new IllegalArgumentException(
131 "variance must be greater than zero");
132 }
133 sig = Math.sqrt(variance);
134 }
135
136 /**
137 * Evaluates the probability density function (p.d.f.) of a Gaussian
138 * distribution having mean mu and standard deviation sig at provided point
139 * x.
140 *
141 * @param x point where p.d.f. is evaluated.
142 * @param mu mean of Gaussian distribution.
143 * @param sig standard deviation of Gaussian distribution.
144 * @return evaluation of p.d.f.
145 * @throws IllegalArgumentException if provided standard deviation is zero
146 * or negative.
147 */
148 public static double p(final double x, final double mu, final double sig) {
149 if (sig <= 0.0) {
150 throw new IllegalArgumentException();
151 }
152
153 return internalP(x, mu, sig);
154 }
155
156 /**
157 * Evaluates the probability density function (p.d.f.) of a Gaussian
158 * distribution having the mean and standard deviation of this instance at
159 * provided point x.
160 *
161 * @param x point where p.d.f. is evaluated.
162 * @return evaluation of p.d.f.
163 */
164 public double p(final double x) {
165 return internalP(x, mu, sig);
166 }
167
168 /**
169 * Evaluates the cumulative distribution function (c.d.f.) of a Gaussian
170 * distribution having mean mu and standard deviation sig at provided point
171 * x.
172 * The c.d.f is equivalent to the probability of the Gaussian distribution
173 * of having a value less than x, and it is computed as the integral from
174 * -infinity to x of the Gaussian p.d.f.
175 *
176 * @param x point where c.d.f. is evaluated.
177 * @param mu mean of Gaussian distribution.
178 * @param sig standard deviation of Gaussian distribution.
179 * @return evaluation of c.d.f.
180 * @throws IllegalArgumentException if provided standard deviation is zero
181 * or negative.
182 */
183 public static double cdf(final double x, final double mu, final double sig) {
184 if (sig <= 0.0) {
185 throw new IllegalArgumentException();
186 }
187
188 return internalCdf(x, mu, sig);
189 }
190
191 /**
192 * Evaluates the cumulative distribution function (c.d.f.) of a Gaussian
193 * distribution having the mean and standard deviation of this instance at
194 * provided point x.
195 * The c.d.f is equivalent to the probability of the Gaussian distribution
196 * of having a value less than x, and it is computed as the integral from
197 * -infinity to x of the Gaussian p.d.f.
198 * Because the c.d.f is a probability, it always returns values between 0.0
199 * and 1.0.
200 *
201 * @param x point where c.d.f. is evaluated.
202 * @return evaluation of c.d.f.
203 */
204 public double cdf(final double x) {
205 return internalCdf(x, mu, sig);
206 }
207
208 /**
209 * Evaluates the inverse cumulative distribution function of a Gaussian
210 * distribution having mean mu and standard deviation sig at provided point
211 * p.
212 * Because the c.d.f is a monotonically increasing function with values
213 * between 0.0 and 1.0, its inverse is uniquely defined between such range
214 * of values.
215 *
216 * @param p value to evaluate the inverse c.d.f. at. This value is
217 * equivalent to a probability and must be between 0.0 and 1.0.
218 * @param mu mean of Gaussian distribution.
219 * @param sig standard deviation of Gaussian distribution.
220 * @return the value x for which the c.d.f. has value p.
221 * @throws IllegalArgumentException if provided standard deviation is zero
222 * or negative, or if provided probability value is not between 0.0 and 1.0.
223 */
224 public static double invcdf(final double p, final double mu, final double sig) {
225 if (sig <= 0.0) {
226 throw new IllegalArgumentException("standard deviation must be greater than zero");
227 }
228
229 return internalInvcdf(p, mu, sig);
230 }
231
232 /**
233 * Evaluates the inverse cumulative distribution function of a Gaussian
234 * distribution having the mean and standard deviation of this instance at
235 * provided point p.
236 * Because the c.d.f is a monotonically increasing function with values
237 * between 0.0 and 1.0, its inverse is uniquely defined between such range
238 * of values.
239 *
240 * @param p value to evaluate the inverse c.d.f. at. This value is
241 * equivalent to a probability and must be between 0.0 and 1.0.
242 * @return the value x for which the c.d.f. has value p.
243 * @throws IllegalArgumentException if provided probability value is not
244 * between 0.0 and 1.0.
245 */
246 public double invcdf(final double p) {
247 return internalInvcdf(p, mu, sig);
248 }
249
250 /**
251 * Computes the Mahalanobis distance of provided point x for provided
252 * mean and standard deviation values.
253 *
254 * @param x point where Mahalanobis distance is evaluated.
255 * @param mu mean of Gaussian distribution.
256 * @param sig standard deviation of Gaussian distribution.
257 * @return Mahalanobis distance of provided point respect to mean.
258 * @throws IllegalArgumentException if provided standard deviation is zero
259 * or negative.
260 */
261 public static double mahalanobisDistance(final double x, final double mu, final double sig) {
262 if (sig <= 0.0) {
263 throw new IllegalArgumentException("standard deviation must be greater than zero");
264 }
265
266 return internalMahalanobisDistance(x, mu, sig);
267 }
268
269 /**
270 * Computes the Mahalanobis distance of provided point x for current mean
271 * and standard deviation values.
272 *
273 * @param x point where Mahalanobis distance is evaluated.
274 * @return Mahalanobis distance of provided point respect to current mean.
275 */
276 public double mahalanobisDistance(final double x) {
277 return internalMahalanobisDistance(x, mu, sig);
278 }
279
280 /**
281 * Evaluates the probability density function (p.d.f.) of a Gaussian
282 * distribution having mean mu and standard deviation sig at provided point
283 * x.
284 * This method is used internally.
285 *
286 * @param x point where p.d.f. is evaluated.
287 * @param mu mean of Gaussian distribution.
288 * @param sig standard deviation of Gaussian distribution.
289 * @return evaluation of p.d.f.
290 * @throws IllegalArgumentException if provided standard deviation is zero
291 * or negative.
292 */
293 private static double internalP(final double x, final double mu, final double sig) {
294 return (GAUSSIAN_NORM / sig) * Math.exp(-0.5 * Math.pow((x - mu) / sig, 2.0));
295 }
296
297 /**
298 * Evaluates the cumulative distribution function (c.d.f.) of a Gaussian
299 * distribution having mean mu and standard deviation sig at provided point
300 * x.
301 * The c.d.f is equivalent to the probability of the Gaussian distribution
302 * of having a value less than x, and it is computed as the integral from
303 * -infinity to x of the Gaussian p.d.f.
304 * This method is used internally.
305 *
306 * @param x point where c.d.f. is evaluated.
307 * @param mu mean of Gaussian distribution.
308 * @param sig standard deviation of Gaussian distribution.
309 * @return evaluation of c.d.f.
310 * @throws IllegalArgumentException if provided standard deviation is zero
311 * or negative.
312 */
313 private static double internalCdf(final double x, final double mu, final double sig) {
314 return 0.5 * Erf.erfc(-HALF_SQRT2 * (x - mu) / sig);
315 }
316
317 /**
318 * Evaluates the inverse cumulative distribution function of a Gaussian
319 * distribution having mean mu and standard deviation sig at provided point
320 * p.
321 * Because the c.d.f is a monotonically increasing function with values
322 * between 0.0 and 1.0, its inverse is uniquely defined between such range
323 * of values.
324 * This method is used internally.
325 *
326 * @param p value to evaluate the inverse c.d.f. at. This value is
327 * equivalent to a probability and must be between 0.0 and 1.0.
328 * @param mu mean of Gaussian distribution.
329 * @param sig standard deviation of Gaussian distribution.
330 * @return the value x for which the c.d.f. has value p.
331 * @throws IllegalArgumentException if provided probability value is not
332 * between 0.0 and 1.0.
333 */
334 private static double internalInvcdf(final double p, final double mu, final double sig)
335 throws IllegalArgumentException {
336 if (p <= 0.0 || p >= 1.0) {
337 throw new IllegalArgumentException("probability value must be between 0.0 and 1.0");
338 }
339 return -SQRT2 * sig * Erf.inverfc(2.0 * p) + mu;
340 }
341
342 /**
343 * Computes the Mahalanobis distance of provided point x for provided
344 * mean and standard deviation values.
345 *
346 * @param x point where Mahalanobis distance is evaluated.
347 * @param mu mean of Gaussian distribution.
348 * @param sig standard deviation of Gaussian distribution.
349 * @return Mahalanobis distance of provided point respect to mean.
350 */
351 private static double internalMahalanobisDistance(final double x, final double mu, final double sig) {
352 return Math.abs(x - mu) / sig;
353 }
354
355 /**
356 * Evaluates the derivative and a 1D function at a certain mean point and
357 * computes the non-linear propagation of Gaussian uncertainty through such
358 * function at such point.
359 *
360 * @param evaluator interface to evaluate derivative of a function at a
361 * certain point.
362 * @param mean mean of original Gaussian distribution to be propagated.
363 * @param standardDeviation standard deviation of original Gaussian
364 * distribution to be propagated.
365 * @param result instance where propagated Gaussian distribution will be
366 * stored.
367 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
368 */
369 public static void propagate(
370 final DerivativeEvaluator evaluator, final double mean, final double standardDeviation,
371 final NormalDist result) {
372 final var evaluation = evaluator.evaluate(mean);
373 final var derivative = evaluator.evaluateDerivative(mean);
374 result.setMean(evaluation);
375 result.setStandardDeviation(Math.abs(derivative * standardDeviation));
376 }
377
378 /**
379 * Evaluates the derivative and a 1D function at a certain mean point and
380 * computes the non-linear propagation of Gaussian uncertainty through such
381 * function at such point.
382 *
383 * @param evaluator interface to evaluate derivative of a function at a
384 * certain point.
385 * @param mean mean of original Gaussian distribution to be propagated.
386 * @param standardDeviation standard deviation of original Gaussian
387 * distribution to be propagated.
388 * @return a new propagated Gaussian distribution.
389 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
390 */
391 public static NormalDist propagate(
392 final DerivativeEvaluator evaluator, final double mean, final double standardDeviation) {
393 final var result = new NormalDist();
394 propagate(evaluator, mean, standardDeviation, result);
395 return result;
396 }
397
398 /**
399 * Evaluates the derivative and a 1D function at a certain mean point and
400 * computes the non-linear propagation of Gaussian uncertainty through such
401 * function at such point.
402 *
403 * @param evaluator interface to evaluate derivative of a function at a
404 * certain point.
405 * @param dist 1D Gaussian distribution to be propagated.
406 * @param result instance where propagated Gaussian distribution will be
407 * stored.
408 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
409 */
410 public static void propagate(final DerivativeEvaluator evaluator, final NormalDist dist, final NormalDist result) {
411 propagate(evaluator, dist.getMean(), dist.getStandardDeviation(), result);
412 }
413
414 /**
415 * Evaluates the derivative and a 1D function at a certain mean point and
416 * computes the non-linear propagation of Gaussian uncertainty through such
417 * function at such point.
418 *
419 * @param evaluator interface to evaluate derivative of a function at a
420 * certain point.
421 * @param dist 1D Gaussian distribution to be propagated.
422 * @return a new propagated Gaussian distribution.
423 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
424 */
425 public static NormalDist propagate(final DerivativeEvaluator evaluator, final NormalDist dist) {
426 final var result = new NormalDist();
427 propagate(evaluator, dist, result);
428 return result;
429 }
430
431 /**
432 * Evaluates the derivative and a 1D function at the mean point of this
433 * normal distribution and computes the non-linear propagation of Gaussian
434 * uncertainty through such function at such point.
435 *
436 * @param evaluator interface to evaluate derivative of a function at the
437 * mean point of this normal distribution.
438 * @param result instance where propagated Gaussian distribution will be
439 * stored.
440 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
441 */
442 public void propagateThisDistribution(final DerivativeEvaluator evaluator, final NormalDist result) {
443 propagate(evaluator, this, result);
444 }
445
446 /**
447 * Evaluates the derivative and a 1D function at the mean point of this
448 * normal distribution and computes the non-linear propagation of Gaussian
449 * uncertainty through such function at such point.
450 *
451 * @param evaluator interface to evaluate derivative of a function at the
452 * mean point of this normal distribution.
453 * @return a new propagated Gaussian distribution.
454 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
455 */
456 public NormalDist propagateThisDistribution(final DerivativeEvaluator evaluator) {
457 final var result = new NormalDist();
458 propagateThisDistribution(evaluator, result);
459 return result;
460 }
461
462 /**
463 * Interface to evaluate a one dimensional function at point x and to obtain
464 * its derivative at such point.
465 */
466 public interface DerivativeEvaluator {
467
468 /**
469 * Evaluates function at point x.
470 *
471 * @param x point x where derivative is evaluated.
472 * @return evaluation of function at point x.
473 */
474 double evaluate(final double x);
475
476 /**
477 * Evaluates derivative of a function at point x.
478 *
479 * @param x point x where derivative is evaluated.
480 * @return derivative of one dimensional function at point x.
481 */
482 double evaluateDerivative(final double x);
483 }
484 }