MultivariateNormalDistribution.java

  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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.distribution.multivariate;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.linear.Array2DRowRealMatrix;
  27. import org.hipparchus.linear.EigenDecompositionSymmetric;
  28. import org.hipparchus.linear.RealMatrix;
  29. import org.hipparchus.random.RandomGenerator;
  30. import org.hipparchus.random.Well19937c;
  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.Precision;

  34. /**
  35.  * Implementation of the multivariate normal (Gaussian) distribution.
  36.  *
  37.  * @see <a href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
  38.  * Multivariate normal distribution (Wikipedia)</a>
  39.  * @see <a href="http://mathworld.wolfram.com/MultivariateNormalDistribution.html">
  40.  * Multivariate normal distribution (MathWorld)</a>
  41.  */
  42. public class MultivariateNormalDistribution
  43.     extends AbstractMultivariateRealDistribution {
  44.     /** Default singular matrix tolerance check value **/
  45.     private static final double DEFAULT_TOLERANCE = Precision.EPSILON;

  47.     /** Vector of means. */
  48.     private final double[] means;
  49.     /** Covariance matrix. */
  50.     private final RealMatrix covarianceMatrix;
  51.     /** The matrix inverse of the covariance matrix. */
  52.     private final RealMatrix covarianceMatrixInverse;
  53.     /** The determinant of the covariance matrix. */
  54.     private final double covarianceMatrixDeterminant;
  55.     /** Matrix used in computation of samples. */
  56.     private final RealMatrix samplingMatrix;
  57.     /** Inverse singular check tolerance when testing if invertable **/
  58.     private final double singularMatrixCheckTolerance;

  60.     /**
  61.      * Creates a multivariate normal distribution with the given mean vector and
  62.      * covariance matrix.<br>
  63.      * The number of dimensions is equal to the length of the mean vector
  64.      * and to the number of rows and columns of the covariance matrix.
  65.      * It is frequently written as "p" in formulae.
  66.      * <p>
  67.      * <b>Note:</b> this constructor will implicitly create an instance of
  68.      * {@link Well19937c} as random generator to be used for sampling only (see
  69.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  70.      * needed for the created distribution, it is advised to pass {@code null}
  71.      * as random generator via the appropriate constructors to avoid the
  72.      * additional initialisation overhead.
  73.      *
  74.      * @param means Vector of means.
  75.      * @param covariances Covariance matrix.
  76.      * @throws MathIllegalArgumentException if the arrays length are
  77.      * inconsistent.
  78.      * @throws MathIllegalArgumentException if the eigenvalue decomposition cannot
  79.      * be performed on the provided covariance matrix.
  80.      * @throws MathIllegalArgumentException if any of the eigenvalues is
  81.      * negative.
  82.      */
  83.     public MultivariateNormalDistribution(final double[] means,
  84.                                           final double[][] covariances)
  85.         throws MathIllegalArgumentException {
  86.         this(means, covariances, DEFAULT_TOLERANCE);
  87.     }

  89.     /**
  90.      * Creates a multivariate normal distribution with the given mean vector and
  91.      * covariance matrix.<br>
  92.      * The number of dimensions is equal to the length of the mean vector
  93.      * and to the number of rows and columns of the covariance matrix.
  94.      * It is frequently written as "p" in formulae.
  95.      * <p>
  96.      * <b>Note:</b> this constructor will implicitly create an instance of
  97.      * {@link Well19937c} as random generator to be used for sampling only (see
  98.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  99.      * needed for the created distribution, it is advised to pass {@code null}
  100.      * as random generator via the appropriate constructors to avoid the
  101.      * additional initialisation overhead.
  102.      *
  103.      * @param means Vector of means.
  104.      * @param covariances Covariance matrix.
  105.      * @param singularMatrixCheckTolerance Tolerance used during the singular matrix check before inversion
  106.      * @throws MathIllegalArgumentException if the arrays length are
  107.      * inconsistent.
  108.      * @throws MathIllegalArgumentException if the eigenvalue decomposition cannot
  109.      * be performed on the provided covariance matrix.
  110.      * @throws MathIllegalArgumentException if any of the eigenvalues is
  111.      * negative.
  112.      */
  113.     public MultivariateNormalDistribution(final double[] means,
  114.                                           final double[][] covariances,
  115.                                           final double singularMatrixCheckTolerance)
  116.         throws MathIllegalArgumentException {
  117.         this(new Well19937c(), means, covariances, singularMatrixCheckTolerance);
  118.     }


  121.     /**
  122.      * Creates a multivariate normal distribution with the given mean vector and
  123.      * covariance matrix.
  124.      * <br>
  125.      * The number of dimensions is equal to the length of the mean vector
  126.      * and to the number of rows and columns of the covariance matrix.
  127.      * It is frequently written as "p" in formulae.
  128.      *
  129.      * @param rng Random Number Generator.
  130.      * @param means Vector of means.
  131.      * @param covariances Covariance matrix.
  132.      * @throws MathIllegalArgumentException if the arrays length are
  133.      * inconsistent.
  134.      * @throws MathIllegalArgumentException if the eigenvalue decomposition cannot
  135.      * be performed on the provided covariance matrix.
  136.      * @throws MathIllegalArgumentException if any of the eigenvalues is
  137.      * negative.
  138.      */
  139.     public MultivariateNormalDistribution(RandomGenerator rng,
  140.                                           final double[] means,
  141.                                           final double[][] covariances) {
  142.         this(rng, means, covariances, DEFAULT_TOLERANCE);
  143.     }

  145.     /**
  146.      * Creates a multivariate normal distribution with the given mean vector and
  147.      * covariance matrix.
  148.      * <br>
  149.      * The number of dimensions is equal to the length of the mean vector
  150.      * and to the number of rows and columns of the covariance matrix.
  151.      * It is frequently written as "p" in formulae.
  152.      *
  153.      * @param rng Random Number Generator.
  154.      * @param means Vector of means.
  155.      * @param covariances Covariance matrix.
  156.      * @param singularMatrixCheckTolerance Tolerance used during the singular matrix check before inversion
  157.      * @throws MathIllegalArgumentException if the arrays length are
  158.      * inconsistent.
  159.      * @throws MathIllegalArgumentException if the eigenvalue decomposition cannot
  160.      * be performed on the provided covariance matrix.
  161.      * @throws MathIllegalArgumentException if any of the eigenvalues is
  162.      * negative.
  163.      */
  164.     public MultivariateNormalDistribution(RandomGenerator rng,
  165.                                           final double[] means,
  166.                                           final double[][] covariances,
  167.                                           final double singularMatrixCheckTolerance)
  168.             throws MathIllegalArgumentException {
  169.         super(rng, means.length);

  171.         final int dim = means.length;

  173.         if (covariances.length != dim) {
  174.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  175.                                                    covariances.length, dim);
  176.         }

  178.         for (int i = 0; i < dim; i++) {
  179.             if (dim != covariances[i].length) {
  180.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  181.                                                        covariances[i].length, dim);
  182.             }
  183.         }

  185.         this.means = means.clone();
  186.         this.singularMatrixCheckTolerance = singularMatrixCheckTolerance;

  188.         covarianceMatrix = new Array2DRowRealMatrix(covariances);

  190.         // Covariance matrix eigen decomposition.
  191.         final EigenDecompositionSymmetric covMatDec =
  192.                         new EigenDecompositionSymmetric(covarianceMatrix, singularMatrixCheckTolerance, true);

  194.         // Compute and store the inverse.
  195.         covarianceMatrixInverse = covMatDec.getSolver().getInverse();
  196.         // Compute and store the determinant.
  197.         covarianceMatrixDeterminant = covMatDec.getDeterminant();

  199.         // Eigenvalues of the covariance matrix.
  200.         final double[] covMatEigenvalues = covMatDec.getEigenvalues();

  202.         for (double covMatEigenvalue : covMatEigenvalues) {
  203.             if (covMatEigenvalue < 0) {
  204.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_POSITIVE_DEFINITE_MATRIX);
  205.             }
  206.         }

  208.         // Matrix where each column is an eigenvector of the covariance matrix.
  209.         final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim);
  210.         for (int v = 0; v < dim; v++) {
  211.             final double[] evec = covMatDec.getEigenvector(v).toArray();
  212.             covMatEigenvectors.setColumn(v, evec);
  213.         }

  215.         final RealMatrix tmpMatrix = covMatEigenvectors.transpose();

  217.         // Scale each eigenvector by the square root of its eigenvalue.
  218.         for (int row = 0; row < dim; row++) {
  219.             final double factor = FastMath.sqrt(covMatEigenvalues[row]);
  220.             for (int col = 0; col < dim; col++) {
  221.                 tmpMatrix.multiplyEntry(row, col, factor);
  222.             }
  223.         }

  225.         samplingMatrix = covMatEigenvectors.multiply(tmpMatrix);
  226.     }

  228.     /**
  229.      * Gets the mean vector.
  230.      *
  231.      * @return the mean vector.
  232.      */
  233.     public double[] getMeans() {
  234.         return means.clone();
  235.     }

  237.     /**
  238.      * Gets the covariance matrix.
  239.      *
  240.      * @return the covariance matrix.
  241.      */
  242.     public RealMatrix getCovariances() {
  243.         return covarianceMatrix.copy();
  244.     }

  246.     /**
  247.      * Gets the current setting for the tolerance check used during singular checks before inversion
  248.      * @return tolerance
  249.      */
  250.     public double getSingularMatrixCheckTolerance() { return singularMatrixCheckTolerance; }

  252.     /** {@inheritDoc} */
  253.     @Override
  254.     public double density(final double[] vals) throws MathIllegalArgumentException {
  255.         final int dim = getDimension();
  256.         if (vals.length != dim) {
  257.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  258.                                                    vals.length, dim);
  259.         }

  261.         return FastMath.pow(2 * FastMath.PI, -0.5 * dim) *
  262.             FastMath.pow(covarianceMatrixDeterminant, -0.5) *
  263.             getExponentTerm(vals);
  264.     }

  266.     /**
  267.      * Gets the square root of each element on the diagonal of the covariance
  268.      * matrix.
  269.      *
  270.      * @return the standard deviations.
  271.      */
  272.     public double[] getStandardDeviations() {
  273.         final int dim = getDimension();
  274.         final double[] std = new double[dim];
  275.         final double[][] s = covarianceMatrix.getData();
  276.         for (int i = 0; i < dim; i++) {
  277.             std[i] = FastMath.sqrt(s[i][i]);
  278.         }
  279.         return std;
  280.     }

  282.     /** {@inheritDoc} */
  283.     @Override
  284.     public double[] sample() {
  285.         final int dim = getDimension();
  286.         final double[] normalVals = new double[dim];

  288.         for (int i = 0; i < dim; i++) {
  289.             normalVals[i] = random.nextGaussian();
  290.         }

  292.         final double[] vals = samplingMatrix.operate(normalVals);

  294.         for (int i = 0; i < dim; i++) {
  295.             vals[i] += means[i];
  296.         }

  298.         return vals;
  299.     }

  301.     /**
  302.      * Computes the term used in the exponent (see definition of the distribution).
  303.      *
  304.      * @param values Values at which to compute density.
  305.      * @return the multiplication factor of density calculations.
  306.      */
  307.     private double getExponentTerm(final double[] values) {
  308.         final double[] centered = new double[values.length];
  309.         for (int i = 0; i < centered.length; i++) {
  310.             centered[i] = values[i] - getMeans()[i];
  311.         }
  312.         final double[] preMultiplied = covarianceMatrixInverse.preMultiply(centered);
  313.         double sum = 0;
  314.         for (int i = 0; i < preMultiplied.length; i++) {
  315.             sum += preMultiplied[i] * centered[i];
  316.         }
  317.         return FastMath.exp(-0.5 * sum);
  318.     }
  319. }
Created with JaCoCo 0.8.12.202403310830