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 import com.irurueta.algebra.*;
19
20 import java.util.Arrays;
21
22 /**
23 * Contains methods to work with multivariate normal (i.e. Gaussian)
24 * distributions.
25 */
26 public class MultivariateNormalDist {
27
28 /**
29 * Mean value of Gaussian distribution.
30 */
31 private double[] mu;
32
33 /**
34 * Covariance of Gaussian distribution.
35 */
36 private Matrix cov;
37
38 /**
39 * Basis in which the covariance matrix is expressed.
40 * This value is obtained after decomposition.
41 */
42 private Matrix covBasis;
43
44 /**
45 * Variances on each direction of the basis.
46 */
47 private double[] variances;
48
49 /**
50 * Constructor.
51 * Creates a normal distribution of 1 dimension.
52 */
53 public MultivariateNormalDist() {
54 this(1);
55 }
56
57 /**
58 * Constructor.
59 * Creates a multivariate normal distribution having the provided
60 * number of dimensions, with zero mean and unitary independent variances.
61 *
62 * @param dims number of dimensions. Must be greater than zero.
63 * @throws IllegalArgumentException if provided number of dimensions is
64 * zero or less.
65 */
66 public MultivariateNormalDist(final int dims) {
67 if (dims <= 0) {
68 throw new IllegalArgumentException("number of dimensions must be greater than zero");
69 }
70
71 mu = new double[dims];
72 try {
73 cov = Matrix.identity(dims, dims);
74 } catch (final WrongSizeException e) {
75 throw new IllegalArgumentException("number of dimensions must be greater than zero", e);
76 }
77 }
78
79 /**
80 * Constructor.
81 * Creates a multivariate normal distribution having provided mean and
82 * covariance.
83 *
84 * @param mean array containing mean. Must have the same number of rows as
85 * provided covariance matrix
86 * @param covariance matrix containing covariance. Must be square, symmetric
87 * and positive definite (i.e. non-singular) and must have the same number
88 * of rows as provided mean.
89 * @throws IllegalArgumentException if provided mean array has length
90 * smaller than 1 or if length of mean array is not the same as the number
91 * of rows of covariance matrix.
92 * @throws InvalidCovarianceMatrixException if provided covariance matrix is
93 * not square, symmetric and positive definite (i.e. non singular).
94 */
95 public MultivariateNormalDist(final double[] mean, final Matrix covariance)
96 throws InvalidCovarianceMatrixException {
97 setMeanAndCovariance(mean, covariance);
98 }
99
100 /**
101 * Constructor.
102 * Creates a multivariate normal distribution having provided mean and
103 * covariance.
104 *
105 * @param mean array containing mean. Must have the same number of rows as
106 * provided covariance matrix.
107 * @param covariance matrix containing covariance. Must be square, symmetric
108 * and positive definite (i.e. non-singular) and must have the same number
109 * of rows as provided mean.
110 * @param validateSymmetricPositiveDefinite true if covariance matrix must
111 * be validated to be positive definite, false to skip validation.
112 * @throws IllegalArgumentException if provided mean array has length
113 * smaller than 1 or if length of mean array is not the same as the number
114 * of rows of covariance matrix.
115 * @throws InvalidCovarianceMatrixException if provided matrix is not
116 * valid (nor square or symmetric positive definite if validation is
117 * enabled).
118 */
119 public MultivariateNormalDist(
120 final double[] mean, final Matrix covariance, final boolean validateSymmetricPositiveDefinite)
121 throws InvalidCovarianceMatrixException {
122 setMeanAndCovariance(mean, covariance, validateSymmetricPositiveDefinite);
123 }
124
125 /**
126 * Gets array containing mean of this multivariate Gaussian distribution.
127 *
128 * @return mean of multivariate Gaussian distribution.
129 */
130 public double[] getMean() {
131 return mu;
132 }
133
134 /**
135 * Sets mean of this multivariate Gaussian distribution.
136 * Length of provided mean must be equal to the number of rows of provided
137 * covariance, otherwise instance won't be ready.
138 *
139 * @param mu mean of multivariate Gaussian distribution.
140 * @throws IllegalArgumentException if provided array has a length smaller
141 * than 1.
142 */
143 public void setMean(final double[] mu) {
144 if (mu.length == 0) {
145 throw new IllegalArgumentException("length of mean array must be greater than zero");
146 }
147 this.mu = mu;
148 }
149
150 /**
151 * Gets matrix containing covariance of this multivariate Gaussian
152 * distribution.
153 *
154 * @return covariance of multivariate Gaussian distribution.
155 */
156 public Matrix getCovariance() {
157 return new Matrix(cov);
158 }
159
160 /**
161 * Gets matrix containing covariance of this multivariate Gaussian
162 * distribution.
163 *
164 * @param result instance where covariance of multivariate Gaussian
165 * distribution will be stored.
166 */
167 public void getCovariance(final Matrix result) {
168 cov.copyTo(result);
169 }
170
171 /**
172 * Sets covariance of this multivariate Gaussian distribution.
173 *
174 * @param cov covariance of this multivariate Gaussian distribution.
175 * @throws InvalidCovarianceMatrixException if provided matrix is not valid
176 * (not square or symmetric positive definite).
177 */
178 public void setCovariance(final Matrix cov) throws InvalidCovarianceMatrixException {
179 setCovariance(cov, true);
180 }
181
182 /**
183 * Sets covariance of this multivariate Gaussian distribution.
184 *
185 * @param cov covariance of this multivariate Gaussian distribution.
186 * @param validateSymmetricPositiveDefinite true if matrix must be
187 * validated to be positive definite, false to skip validation.
188 * @throws InvalidCovarianceMatrixException if provided matrix is not
189 * valid (nor square or symmetric positive definite if validation is
190 * enabled).
191 */
192 public void setCovariance(final Matrix cov, final boolean validateSymmetricPositiveDefinite)
193 throws InvalidCovarianceMatrixException {
194 if (cov.getRows() != cov.getColumns()) {
195 throw new InvalidCovarianceMatrixException("covariance matrix must be square");
196 }
197
198 try {
199 if (validateSymmetricPositiveDefinite) {
200 final var decomposer = new CholeskyDecomposer(cov);
201 decomposer.decompose();
202 if (!decomposer.isSPD()) {
203 throw new InvalidCovarianceMatrixException(
204 "covariance matrix must be symmetric positive definite (non singular)");
205 }
206 }
207
208 this.cov = new Matrix(cov);
209 covBasis = null;
210 variances = null;
211 } catch (final AlgebraException e) {
212 throw new InvalidCovarianceMatrixException("covariance matrix must be square", e);
213 }
214 }
215
216 /**
217 * Sets mean and covariance of this multivariate Gaussian distribution.
218 *
219 * @param mu array containing mean. Must have the same number of rows as
220 * provided covariance matrix
221 * @param cov matrix containing covariance. Must be square, symmetric
222 * and positive definite (i.e. non-singular) and must have the same number
223 * of rows as provided mean.
224 * @throws IllegalArgumentException if provided mean array has length
225 * smaller than 1 or if length of mean array is not the same as the number
226 * of rows of covariance matrix.
227 * @throws InvalidCovarianceMatrixException if provided covariance matrix is
228 * not square, symmetric and positive definite (i.e. non singular).
229 */
230 public final void setMeanAndCovariance(final double[] mu, final Matrix cov)
231 throws InvalidCovarianceMatrixException {
232 setMeanAndCovariance(mu, cov, true);
233 }
234
235 /**
236 * Sets mean and covariance of this multivariate Gaussian distribution.
237 *
238 * @param mu array containing mean. Must have the same number of rows as
239 * provided covariance matrix
240 * @param cov matrix containing covariance. Must be square, symmetric
241 * and positive definite (i.e. non-singular) and must have the same number
242 * of rows as provided mean.
243 * @param validateSymmetricPositiveDefinite true if matrix must be
244 * validated to be positive definite, false to skip validation.
245 * @throws IllegalArgumentException if provided mean array has length
246 * smaller than 1 or if length of mean array is not the same as the number
247 * of rows of covariance matrix.
248 * @throws InvalidCovarianceMatrixException if provided covariance matrix is
249 * not square, symmetric and positive definite (i.e. non singular).
250 */
251 public final void setMeanAndCovariance(
252 final double[] mu, final Matrix cov, final boolean validateSymmetricPositiveDefinite)
253 throws InvalidCovarianceMatrixException {
254 if (mu.length != cov.getRows()) {
255 throw new IllegalArgumentException("mean array length must be equal to covariance number of rows");
256 }
257
258 setCovariance(cov, validateSymmetricPositiveDefinite);
259 setMean(mu);
260 }
261
262 /**
263 * Indicates whether provided matrix is a valid covariance matrix.
264 * A valid covariance matrix must be square, symmetric and positive definite
265 * (i.e. non-singular).
266 *
267 * @param cov matrix to be checked.
268 * @return true if matrix is a valid covariance matrix, false otherwise.
269 */
270 public static boolean isValidCovariance(final Matrix cov) {
271 if (cov.getRows() != cov.getColumns()) {
272 return false;
273 }
274
275 try {
276 final var decomposer = new CholeskyDecomposer(cov);
277 decomposer.decompose();
278 return decomposer.isSPD();
279 } catch (final AlgebraException e) {
280 return false;
281 }
282 }
283
284 /**
285 * Indicates whether this instance is ready for any computation, false
286 * otherwise.
287 *
288 * @return true if instance is ready, false otherwise.
289 */
290 public boolean isReady() {
291 return mu != null && cov != null &&
292 mu.length == cov.getRows();
293 }
294
295 /**
296 * Basis containing on each column the direction of each variance in the
297 * multidimensional Gaussian distribution, which is obtained from provided
298 * covariance matrix.
299 * This value is available only after the p.d.f. has been evaluated.
300 *
301 * @return basis containing on each column the direction of each variance
302 * in the multidimensional Gaussian distribution.
303 */
304 public Matrix getCovarianceBasis() {
305 return covBasis;
306 }
307
308 /**
309 * Array containing the amount of variance on each direction of the basis
310 * of the covariance in the multidimensional Gaussian distribution.
311 * This value is available only after the p.d.f. has been evaluated.
312 *
313 * @return variance on each direction of the basis of the covariance.
314 */
315 public double[] getVariances() {
316 return variances;
317 }
318
319 /**
320 * Evaluates the probability density function (p.d.f.) of a multivariate
321 * Gaussian distribution having current mean and covariance at point x.
322 *
323 * @param x array containing coordinates where p.d.f. is evaluated.
324 * @return evaluation of p.d.f.
325 * @throws NotReadyException if this instance is not ready (mean and
326 * covariance have not been provided or are not valid).
327 * @throws IllegalArgumentException if provided point length is not valid.
328 * @throws DecomposerException happens if covariance is numerically
329 * unstable (i.e. contains NaNs or very large numbers).
330 * @throws RankDeficientMatrixException happens if covariance is singular.
331 */
332 public double p(final double[] x) throws NotReadyException, DecomposerException, RankDeficientMatrixException {
333 if (!isReady()) {
334 throw new NotReadyException();
335 }
336
337 final var k = x.length;
338 if (k != mu.length) {
339 throw new IllegalArgumentException("length of point must be equal to the length of mean");
340 }
341
342 var detCov = 0.0;
343 try {
344 detCov = Utils.det(cov);
345 } catch (final WrongSizeException ignore) {
346 // never thrown
347 }
348
349 final var factor = 1.0 / (Math.sqrt(Math.pow(2.0 * Math.PI, k) * detCov));
350 return factor * Math.exp(-0.5 * squaredMahalanobisDistance(x));
351 }
352
353 /**
354 * Evaluates the cumulative distribution function (c.d.f.) of a Gaussian
355 * distribution having current mean and covariance values.
356 * The c.d.f. is equivalent to the joint probability of the multivariate
357 * Gaussian distribution of having a value less than x on each direction
358 * of the basis of independent variances obtained from covariance matrix.
359 * Because the c.d.f is a probability, it always returns values between 0.0
360 * and 1.0.
361 * NOTE: this method will resize provided basis instance if needed.
362 *
363 * @param x point where c.d.f. is evaluated.
364 * @param basis instance where is stored the basis of each direction of
365 * independent covariances, if provided.
366 * @return evaluation of c.d.f.
367 * @throws IllegalArgumentException if length of provided point is not equal
368 * to length of current mean.
369 * @throws NotReadyException if this instance is not ready (mean and
370 * covariance have not been provided or are not valid).
371 * @throws DecomposerException if covariance is numerically unstable (i.e.
372 * contains NaNs or very large numbers).
373 */
374 public double cdf(final double[] x, final Matrix basis) throws NotReadyException, DecomposerException {
375 if (!isReady()) {
376 throw new NotReadyException();
377 }
378
379 final var k = x.length;
380 if (k != mu.length) {
381 throw new IllegalArgumentException("length of point must be equal to the length of mean");
382 }
383
384 var p = 1.0;
385 try {
386 processCovariance();
387
388 if (basis != null) {
389 basis.copyFrom(covBasis);
390 }
391
392 for (int i = 0; i < k; i++) {
393 final var singleBasis = covBasis.getSubmatrixAsArray(0, i, k - 1, i);
394 final var coordX = ArrayUtils.dotProduct(x, singleBasis);
395 final var coordMu = ArrayUtils.dotProduct(mu, singleBasis);
396 p *= NormalDist.cdf(coordX, coordMu, Math.sqrt(variances[i]));
397 }
398
399 } catch (final DecomposerException e) {
400 throw e;
401 } catch (final AlgebraException ignore) {
402 // never thrown
403 }
404
405 return p;
406 }
407
408 /**
409 * Evaluates the cumulative distribution function (c.d.f.) of a Gaussian
410 * distribution having current mean and covariance values.
411 * The c.d.f. is equivalent to the joint probability of the multivariate
412 * Gaussian distribution of having a value less than x on each direction
413 * of the basis of independent variances obtained from covariance matrix.
414 * Because the c.d.f is a probability, it always returns values between 0.0
415 * and 1.0.
416 *
417 * @param x point where c.d.f. is evaluated.
418 * @return evaluation of c.d.f.
419 * @throws IllegalArgumentException if length of provided point is not equal
420 * to length of current mean.
421 * @throws NotReadyException if this instance is not ready (mean and
422 * covariance have not been provided or are not valid).
423 * @throws DecomposerException if covariance is numerically unstable (i.e.
424 * contains NaNs or very large numbers).
425 */
426 public double cdf(double[] x) throws NotReadyException, DecomposerException {
427 return cdf(x, null);
428 }
429
430 /**
431 * Computes the joint probability of all probabilities provided in the
432 * array. The joint probability is computed by multiplying all components of
433 * the array, assuming that all probabilities are independent.
434 *
435 * @param p array containing probabilities for each independent variance
436 * direction that can be obtained from provided covariance matrix.
437 * @return joint probability.
438 */
439 public static double jointProbability(final double[] p) {
440 var jointP = 1.0;
441 for (final var aP : p) {
442 jointP *= aP;
443 }
444 return jointP;
445 }
446
447 /**
448 * Evaluates the inverse cumulative distribution function of a multivariate
449 * Gaussian distribution for current mean and covariance values and provided
450 * probability values for each dimension of the multivariate Gaussian
451 * distribution.
452 * NOTE: this method will resize provided basis instance if needed.
453 *
454 * @param p array containing probability values to evaluate the inverse
455 * c.d.f. on each dimension. Values in the array must be between 0.0 and
456 * 1.0.
457 * @param result coordinates of the value x for which the c.d.f. has values
458 * p.
459 * @param basis instance where is stored the basis of each direction of
460 * independent covariances, if provided.
461 * @throws IllegalArgumentException if length of probabilities is not equal
462 * to mean length, or if result and length of probabilities are not equal,
463 * or if provided probabilities are not between 0.0 and 1.0.
464 * @throws NotReadyException if this instance is not ready (mean and
465 * covariance have not been provided or are not valid).
466 * @throws DecomposerException if covariance is numerically unstable (i.e.
467 * contains NaNs or very large numbers).
468 */
469 public void invcdf(final double[] p, final double[] result, final Matrix basis) throws NotReadyException,
470 DecomposerException {
471 if (!isReady()) {
472 throw new NotReadyException("mean and covariance not provided or invalid");
473 }
474
475 final var k = p.length;
476 if (k != mu.length) {
477 throw new IllegalArgumentException(
478 "length of probabilities must be equal to the length of mean");
479 }
480 if (k != result.length) {
481 throw new IllegalArgumentException("length of result must be equal to the length of mean");
482 }
483
484 try {
485 processCovariance();
486
487 if (basis != null) {
488 basis.copyFrom(covBasis);
489 }
490
491 // initialize to mean
492 System.arraycopy(mu, 0, result, 0, k);
493 for (var i = 0; i < k; i++) {
494 final var singleBasis = covBasis.getSubmatrixAsArray(0, i, k - 1, i);
495 final double coord = NormalDist.invcdf(p[i], mu[i], Math.sqrt(variances[i])) - mu[i];
496 // coord*singleBasis
497 ArrayUtils.multiplyByScalar(singleBasis, coord, singleBasis);
498
499 // result = mean + coord*singleBasis
500 ArrayUtils.sum(result, singleBasis, result);
501 }
502 } catch (final DecomposerException e) {
503 throw e;
504 } catch (final AlgebraException ignore) {
505 // never thrown
506 }
507 }
508
509 /**
510 * Evaluates the inverse cumulative distribution function of a multivariate
511 * Gaussian distribution for current mean and covariance values and provided
512 * probability values for each dimension of the multivariate Gaussian
513 * distribution.
514 * NOTE: this method will resize provided basis instance if needed.
515 *
516 * @param p array containing probability values to evaluate the inverse
517 * c.d.f. on each dimension. Values in the array must be between 0.0 and
518 * 1.0.
519 * @param basis instance where is stored the basis of each direction of
520 * independent covariances, if provided.
521 * @return a new array containing coordinates of the value x for which the
522 * c.d.f. has values p.
523 * @throws IllegalArgumentException if length of probabilities is not equal
524 * to mean length, or if provided probabilities are not between 0.0 and 1.0.
525 * @throws NotReadyException if this instance is not ready (mean and
526 * covariance have not been provided or are not valid).
527 * @throws DecomposerException if covariance is numerically unstable (i.e.
528 * contains NaNs or very large numbers).
529 */
530 public double[] invcdf(final double[] p, final Matrix basis) throws NotReadyException, DecomposerException {
531 if (mu == null) {
532 throw new NotReadyException("mean not defined");
533 }
534
535 final var result = new double[mu.length];
536 invcdf(p, result, basis);
537 return result;
538 }
539
540 /**
541 * Evaluates the inverse cumulative distribution function of a multivariate
542 * Gaussian distribution for current mean and covariance values and provided
543 * probability values for each dimension of the multivariate Gaussian
544 * distribution.
545 *
546 * @param p array containing probability values to evaluate the inverse
547 * c.d.f. on each dimension. Values in the array must be between 0.0 and
548 * 1.0.
549 * @param result coordinates of the value x for which the c.d.f. has values
550 * p.
551 * @throws IllegalArgumentException if length of probabilities is not equal
552 * to mean length, or if result and length of probabilities are not equal,
553 * or if provided probabilities are not between 0.0 and 1.0.
554 * @throws NotReadyException if this instance is not ready (mean and
555 * covariance have not been provided or are not valid).
556 * @throws DecomposerException if covariance is numerically unstable (i.e.
557 * contains NaNs or very large numbers).
558 */
559 public void invcdf(final double[] p, final double[] result) throws NotReadyException, DecomposerException {
560 invcdf(p, result, null);
561 }
562
563 /**
564 * Evaluates the inverse cumulative distribution function of a multivariate
565 * Gaussian distribution for current mean and covariance values and provided
566 * probability values for each dimension of the multivariate Gaussian
567 * distribution.
568 *
569 * @param p array containing probability values to evaluate the inverse
570 * c.d.f. on each dimension. Values in the array must be between 0.0 and
571 * 1.0.
572 * @return coordinates of the value x for which the c.d.f. has values
573 * * p.
574 * @throws IllegalArgumentException if length of probabilities is not equal
575 * to mean length, or if result and length of probabilities are not equal,
576 * or if provided probabilities are not between 0.0 and 1.0.
577 * @throws NotReadyException if this instance is not ready (mean and
578 * covariance have not been provided or are not valid).
579 * @throws DecomposerException if covariance is numerically unstable (i.e.
580 * contains NaNs or very large numbers).
581 */
582 public double[] invcdf(final double[] p) throws NotReadyException, DecomposerException {
583 return invcdf(p, (Matrix) null);
584 }
585
586 /**
587 * Evaluates the inverse cumulative distribution function of a multivariate
588 * Gaussian distribution for current mean and covariance values and provided
589 * probability value.
590 * Obtained result coordinates are computed taking into account the basis
591 * of independent variances computed from current covariance matrix.
592 * NOTE: notice that the inverse cdf of a mutivariate Gaussian distribution
593 * does not have a unique solution. This method simply returns one of the
594 * possible solutions by assuming equal probabilities on each dimension.
595 * NOTE: this method will resize provided basis instance if needed.
596 *
597 * @param p probability value to evaluate the inverse c.d.f. at. This value
598 * must be between 0.0 and 1.0
599 * @param result coordinates of the value x for which the c.d.f. has value
600 * p.
601 * @param basis instance where is stored the basis of each direction of
602 * independent covariances, if provided.
603 * @throws IllegalArgumentException if provided probability value is not
604 * between 0.0 and 1.0, if length of provided result array is not equal
605 * to length of current mean.
606 * @throws NotReadyException if this instance is not ready (mean and
607 * covariance have not been provided or are not valid).
608 * @throws DecomposerException if covariance is numerically unstable (i.e.
609 * contains NaNs or very large numbers).
610 */
611 public void invcdf(final double p, final double[] result, final Matrix basis)
612 throws NotReadyException, DecomposerException {
613 if (p <= 0.0 || p >= 1.0) {
614 throw new IllegalArgumentException("probability value must be between 0.0 and 1.0");
615 }
616
617 if (!isReady()) {
618 throw new NotReadyException("mean and covariance not provided or invalid");
619 }
620
621 final var k = result.length;
622 if (k != mu.length) {
623 throw new IllegalArgumentException("length of result must be equal to mean length");
624 }
625
626 final var probs = new double[k];
627 Arrays.fill(probs, Math.pow(p, 1.0 / k));
628 invcdf(probs, result, basis);
629 }
630
631 /**
632 * Evaluates the inverse cumulative distribution function of a multivariate
633 * Gaussian distribution for current mean and covariance values and provided
634 * probability value.
635 * Obtained result coordinates are computed taking into account the basis
636 * of independent variances computed from current covariance matrix.
637 * NOTE: notice that the inverse cdf of a mutivariate Gaussian distribution
638 * does not have a unique solution. This method simply returns one of the
639 * possible solutions by assuming equal probabilities on each dimension.
640 * NOTE: this method will resize provided basis instance if needed.
641 *
642 * @param p probability value to evaluate the inverse c.d.f. at. This value
643 * must be between 0.0 and 1.0
644 * @param basis instance where is stored the basis of each direction of
645 * independent covariances, if provided.
646 * @return a new array containing the coordinates of the value x for which
647 * the c.d.f. has value p.
648 * @throws IllegalArgumentException if provided probability value is not
649 * between 0.0 and 1.0.
650 * @throws NotReadyException if this instance is not ready (mean and
651 * covariance have not been provided or are not valid).
652 * @throws DecomposerException f covariance is numerically unstable (i.e.
653 * contains NaNs or very large numbers).
654 */
655 public double[] invcdf(final double p, final Matrix basis) throws NotReadyException, DecomposerException {
656 if (mu == null) {
657 throw new NotReadyException("mean not defined");
658 }
659
660 final var result = new double[mu.length];
661 invcdf(p, result, basis);
662 return result;
663 }
664
665 /**
666 * Evaluates the inverse cumulative distribution function of a multivariate
667 * Gaussian distribution for current mean and covariance values and provided
668 * probability value.
669 * Obtained result coordinates are computed taking into account the basis
670 * of independent variances computed from current covariance matrix.
671 * NOTE: notice that the inverse cdf of a mutivariate Gaussian distribution
672 * does not have a unique solution. This method simply returns one of the
673 * possible solutions by assuming equal probabilities on each dimension.
674 *
675 * @param p probability value to evaluate the inverse c.d.f. at. This value
676 * must be between 0.0 and 1.0
677 * @param result coordinates of the value x for which the c.d.f. has value
678 * p.
679 * @throws IllegalArgumentException if provided probability value is not
680 * between 0.0 and 1.0, if length of provided result array is not equal
681 * to length of current mean.
682 * @throws NotReadyException if this instance is not ready (mean and
683 * covariance have not been provided or are not valid).
684 * @throws DecomposerException f covariance is numerically unstable (i.e.
685 * contains NaNs or very large numbers).
686 */
687 public void invcdf(final double p, final double[] result) throws NotReadyException, DecomposerException {
688 invcdf(p, result, null);
689 }
690
691 /**
692 * Evaluates the inverse cumulative distribution function of a multivariate
693 * Gaussian distribution for current mean and covariance values and provided
694 * probability value.
695 * Obtained result coordinates are computed taking into account the basis
696 * of independent variances computed from current covariance matrix.
697 * NOTE: notice that the inverse cdf of a mutivariate Gaussian distribution
698 * does not have a unique solution. This method simply returns one of the
699 * possible solutions by assuming equal probabilities on each dimension.
700 *
701 * @param p probability value to evaluate the inverse c.d.f. at. This value
702 * must be between 0.0 and 1.0
703 * @return a new array containing the coordinates of the value x for which
704 * the c.d.f. has value p.
705 * @throws IllegalArgumentException if provided probability value is not
706 * between 0.0 and 1.0.
707 * @throws NotReadyException if this instance is not ready (mean and
708 * covariance have not been provided or are not valid).
709 * @throws DecomposerException f covariance is numerically unstable (i.e.
710 * contains NaNs or very large numbers).
711 */
712 public double[] invcdf(final double p) throws NotReadyException, DecomposerException {
713 return invcdf(p, (Matrix) null);
714 }
715
716 /**
717 * Computes the Mahalanobis distance of provided multivariate pot x for
718 * current mean and covariance values.
719 *
720 * @param x point where Mahalanobis distance is evaluated.
721 * @return Mahalanobis distance of provided point respect to mean.
722 * @throws DecomposerException happens if covariance is numerically
723 * unstable (i.e. contains NaNs or very large numbers).
724 * @throws RankDeficientMatrixException happens if covariance is singular.
725 */
726 public double mahalanobisDistance(final double[] x) throws DecomposerException, RankDeficientMatrixException {
727 return Math.sqrt(squaredMahalanobisDistance(x));
728 }
729
730 /**
731 * Computes the squared Mahalanobis distance of provided multivariate pot x
732 * for current mean and covariance values.
733 *
734 * @param x point where Mahalanobis distance is evaluated.
735 * @return Mahalanobis distance of provided point respect to mean.
736 * @throws DecomposerException happens if covariance is numerically
737 * unstable (i.e. contains NaNs or very large numbers).
738 * @throws RankDeficientMatrixException happens if covariance is singular.
739 */
740 public double squaredMahalanobisDistance(final double[] x) throws DecomposerException,
741 RankDeficientMatrixException {
742 final var diff = ArrayUtils.subtractAndReturnNew(x, mu);
743 final var diffMatrix = Matrix.newFromArray(diff, true);
744 final var transDiffMatrix = diffMatrix.transposeAndReturnNew();
745
746 try {
747 final var invCov = Utils.inverse(cov);
748 transDiffMatrix.multiply(invCov);
749 transDiffMatrix.multiply(diffMatrix);
750
751 } catch (final WrongSizeException ignore) {
752 // never thrown
753 }
754
755 return transDiffMatrix.getElementAtIndex(0);
756 }
757
758 /**
759 * Processes current covariance by decomposing it into a basis and its
760 * corresponding variances if needed.
761 *
762 * @throws DecomposerException happens if covariance is numerically
763 * unstable (i.e. contains NaNs or very large numbers).
764 * @throws NotReadyException never thrown because decomposer will always be
765 * ready.
766 * @throws LockedException never thrown because decomposer will never be
767 * locked.
768 * @throws NotAvailableException never thrown because first a
769 * DecomposerException will be thrown before attempting to get V or
770 * singular values.
771 */
772 public void processCovariance() throws DecomposerException, NotReadyException, LockedException,
773 NotAvailableException {
774 if (cov == null) {
775 throw new NotReadyException("covariance must be defined");
776 }
777
778 if (covBasis == null || variances == null) {
779 final var decomposer = new SingularValueDecomposer(cov);
780 decomposer.decompose();
781
782 // because matrix is symmetric positive definite:
783 // And matrices U and V are orthonormal
784 // Cov = A'*A = (U*S*V')'*(U*S*V')=V*S*U'*U*S*V' = V*S^2*V',
785
786 // where matrix S is diagonal, and contains the standard deviations
787 // on each direction of the basis V, and hence S^2 is also diagonal but
788 // containing variances on each direction.
789 // The values of S^2 are the eigenvalues of Cov, and V are the
790 // eigenvectors of Cov, hence covariance can be expressed as variances
791 // on each direction of the basis V.
792
793 // matrix containing eigenvectors (basis of directions)
794 covBasis = decomposer.getV();
795
796 // array containing the eigenvalues (variances on each direction)
797 variances = decomposer.getSingularValues();
798 }
799 }
800
801 /**
802 * Evaluates the Jacobian and a multivariate function at a certain mean
803 * point and computes the non-linear propagation of Gaussian uncertainty
804 * through such function at such point.
805 *
806 * @param evaluator interface to evaluate a multivariate function and its
807 * Jacobian at a certain point.
808 * @param mean mean of original multivariate Gaussian distribution to be
809 * propagated. Must have the length of the number of input variables of the
810 * multivariate function to be evaluated.
811 * @param covariance covariance of original Gaussian distribution to be
812 * propagated. Must be symmetric positive definite having size NxN where N
813 * is the length of provided mean.
814 * @param result instance where propagated multiavariate Gaussian
815 * distribution will be stored.
816 * @throws WrongSizeException if evaluator returns an invalid number of
817 * variables (i.e. negative or zero).
818 * @throws InvalidCovarianceMatrixException if provided covariance matrix is
819 * not valid (i.e. is not symmetric positive definite).
820 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
821 */
822 public static void propagate(
823 final JacobianEvaluator evaluator, final double[] mean, final Matrix covariance,
824 final MultivariateNormalDist result) throws WrongSizeException, InvalidCovarianceMatrixException {
825
826 final var ndims = mean.length;
827 final var nvars = evaluator.getNumberOfVariables();
828 final var evaluation = new double[nvars];
829 final var jacobian = new Matrix(nvars, ndims);
830 evaluator.evaluate(mean, evaluation, jacobian);
831
832 // [y, Y_x] = f(x)
833 // Y = Y_x * X * Y_x'
834 final var jacobianTrans = jacobian.transposeAndReturnNew();
835 jacobian.multiply(covariance);
836 jacobian.multiply(jacobianTrans);
837
838 // ensure that new covariance is symmetric positive definite
839 jacobian.symmetrize();
840
841 result.setMean(evaluation);
842 result.setCovariance(jacobian, false);
843 }
844
845 /**
846 * Evaluates the Jacobian and a multivariate function at a certain mean
847 * point and computes the non-linear propagation of Gaussian uncertainty
848 * through such function at such point.
849 *
850 * @param evaluator interface to evaluate a multivariate function and its
851 * Jacobian at a certain point.
852 * @param mean mean of original multivariate Gaussian distribution to be
853 * propagated. Must have the length of the number of input variables of the
854 * multivariate function to be evaluated.
855 * @param covariance covariance of original Gaussian distribution to be
856 * propagated. Must be symmetric positive definite having size NxN where N
857 * is the length of provided mean.
858 * @return a new propagated multivariate Gaussian distribution.
859 * @throws WrongSizeException if evaluator returns an invalid number of
860 * variables (i.e. negative or zero).
861 * @throws InvalidCovarianceMatrixException if provided covariance matrix is
862 * not valid (i.e. is not symmetric positive definite).
863 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
864 */
865 public static MultivariateNormalDist propagate(
866 final JacobianEvaluator evaluator, final double[] mean, final Matrix covariance)
867 throws WrongSizeException, InvalidCovarianceMatrixException {
868 final var result = new MultivariateNormalDist();
869 propagate(evaluator, mean, covariance, result);
870 return result;
871 }
872
873 /**
874 * Evaluates the Jacobian and a multivariate function at a certain mean
875 * point and computes the non-linear propagation of Gaussian uncertainty
876 * through such function at such point.
877 *
878 * @param evaluator interface to evaluate a multivariate function and its
879 * Jacobian at a certain point.
880 * @param dist multivariate Gaussian distribution to be propagated.
881 * @param result instance where propagated multivariate Gaussian
882 * distribution will be stored.
883 * @throws WrongSizeException if evaluator returns an invalid number of
884 * variables (i.e. negative or zero).
885 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
886 */
887 public static void propagate(final JacobianEvaluator evaluator, final MultivariateNormalDist dist,
888 final MultivariateNormalDist result) throws WrongSizeException {
889 try {
890 propagate(evaluator, dist.getMean(), dist.getCovariance(), result);
891 } catch (final InvalidCovarianceMatrixException ignore) {
892 // never thrown
893 }
894 }
895
896 /**
897 * Evaluates the Jacobian and a multivariate function at a certain mean
898 * point and computes the non-linear propagation of Gaussian uncertainty
899 * through such function at such point.
900 *
901 * @param evaluator interface to evaluate a multivariate function and its
902 * Jacobian at a certain point.
903 * @param dist multivariate Gaussian distribution to be propagated.
904 * @return a new propagated multivariate Gaussian distribution.
905 * @throws WrongSizeException if evaluator returns an invalid number of
906 * variables (i.e. negative or zero).
907 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
908 */
909 public static MultivariateNormalDist propagate(
910 final JacobianEvaluator evaluator, final MultivariateNormalDist dist) throws WrongSizeException {
911 final var result = new MultivariateNormalDist();
912 propagate(evaluator, dist, result);
913 return result;
914 }
915
916 /**
917 * Evaluates the Jacobian and a multivariate function at the mean point of
918 * this distribution and computes the non-linear propagation of Gaussian
919 * uncertainty through such function at such point.
920 *
921 * @param evaluator interface to evaluate a multivariate function and its
922 * Jacobian at a certain point.
923 * @param result instance where propagated multivariate Gaussian
924 * distribution will be stored.
925 * @throws WrongSizeException if evaluator returns an invalid number of
926 * variables (i.e. negative or zero).
927 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
928 */
929 public void propagateThisDistribution(final JacobianEvaluator evaluator, final MultivariateNormalDist result)
930 throws WrongSizeException {
931 propagate(evaluator, this, result);
932 }
933
934 /**
935 * Evaluates the Jacobian and a multivariate function at the mean point of
936 * this distribution and computes the non-linear propagation of Gaussian
937 * uncertainty through such function at such point.
938 *
939 * @param evaluator interface to evaluate a multivariate function and its
940 * Jacobian at a certain point.
941 * @return a new propagated multivariate Gaussian distribution.
942 * @throws WrongSizeException if evaluator returns an invalid number of
943 * variables (i.e. negative or zero).
944 * @see <a href="https://github.com/joansola/slamtb">propagateUncertainty.m at https://github.com/joansola/slamtb</a>
945 */
946 public MultivariateNormalDist propagateThisDistribution(final JacobianEvaluator evaluator)
947 throws WrongSizeException {
948 final var result = new MultivariateNormalDist();
949 propagateThisDistribution(evaluator, result);
950 return result;
951 }
952
953 /**
954 * Interface to evaluate a multivariate function at multivariate point x to
955 * obtain multivariate result y and its corresponding jacobian at point x.
956 */
957 public interface JacobianEvaluator {
958 /**
959 * Evaluates multivariate point
960 *
961 * @param x array containing multivariate point where function is
962 * evaluated.
963 * @param y result of evaluating multivariate point.
964 * @param jacobian jacobian of multivariate function at point x.
965 */
966 void evaluate(final double[] x, final double[] y, final Matrix jacobian);
967
968 /**
969 * Number of variables in output of evaluated function. This is equal
970 * to the length of the array y obtained as function evaluations.
971 *
972 * @return number of variables of the function.
973 */
974 int getNumberOfVariables();
975 }
976 }