    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.distribution;
    018
    019    import org.apache.commons.math3.exception.DimensionMismatchException;
    020    import org.apache.commons.math3.linear.Array2DRowRealMatrix;
    021    import org.apache.commons.math3.linear.EigenDecomposition;
    022    import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
    023    import org.apache.commons.math3.linear.RealMatrix;
    024    import org.apache.commons.math3.linear.SingularMatrixException;
    025    import org.apache.commons.math3.random.RandomGenerator;
    026    import org.apache.commons.math3.random.Well19937c;
    027    import org.apache.commons.math3.util.FastMath;
    028    import org.apache.commons.math3.util.MathArrays;
    029
    030    /**
    031     * Implementation of the multivariate normal (Gaussian) distribution.
    032     *
    033     * @see <a href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
    034     * Multivariate normal distribution (Wikipedia)</a>
    035     * @see <a href="http://mathworld.wolfram.com/MultivariateNormalDistribution.html">
    036     * Multivariate normal distribution (MathWorld)</a>
    037     *
    038     * @version $Id: MultivariateNormalDistribution.java 1416643 2012-12-03 19:37:14Z tn $
    039     * @since 3.1
    040     */
    041    public class MultivariateNormalDistribution
    042        extends AbstractMultivariateRealDistribution {
    043        /** Vector of means. */
    044        private final double[] means;
    045        /** Covariance matrix. */
    046        private final RealMatrix covarianceMatrix;
    047        /** The matrix inverse of the covariance matrix. */
    048        private final RealMatrix covarianceMatrixInverse;
    049        /** The determinant of the covariance matrix. */
    050        private final double covarianceMatrixDeterminant;
    051        /** Matrix used in computation of samples. */
    052        private final RealMatrix samplingMatrix;
    053
    054        /**
    055         * Creates a multivariate normal distribution with the given mean vector and
    056         * covariance matrix.
    057         * <br/>
    058         * The number of dimensions is equal to the length of the mean vector
    059         * and to the number of rows and columns of the covariance matrix.
    060         * It is frequently written as "p" in formulae.
    061         *
    062         * @param means Vector of means.
    063         * @param covariances Covariance matrix.
    064         * @throws DimensionMismatchException if the arrays length are
    065         * inconsistent.
    066         * @throws SingularMatrixException if the eigenvalue decomposition cannot
    067         * be performed on the provided covariance matrix.
    068         * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is
    069         * negative.
    070         */
    071        public MultivariateNormalDistribution(final double[] means,
    072                                              final double[][] covariances)
    073            throws SingularMatrixException,
    074                   DimensionMismatchException,
    075                   NonPositiveDefiniteMatrixException {
    076            this(new Well19937c(), means, covariances);
    077        }
    078
    079        /**
    080         * Creates a multivariate normal distribution with the given mean vector and
    081         * covariance matrix.
    082         * <br/>
    083         * The number of dimensions is equal to the length of the mean vector
    084         * and to the number of rows and columns of the covariance matrix.
    085         * It is frequently written as "p" in formulae.
    086         *
    087         * @param rng Random Number Generator.
    088         * @param means Vector of means.
    089         * @param covariances Covariance matrix.
    090         * @throws DimensionMismatchException if the arrays length are
    091         * inconsistent.
    092         * @throws SingularMatrixException if the eigenvalue decomposition cannot
    093         * be performed on the provided covariance matrix.
    094         * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is
    095         * negative.
    096         */
    097        public MultivariateNormalDistribution(RandomGenerator rng,
    098                                              final double[] means,
    099                                              final double[][] covariances)
    100                throws SingularMatrixException,
    101                       DimensionMismatchException,
    102                       NonPositiveDefiniteMatrixException {
    103            super(rng, means.length);
    104
    105            final int dim = means.length;
    106
    107            if (covariances.length != dim) {
    108                throw new DimensionMismatchException(covariances.length, dim);
    109            }
    110
    111            for (int i = 0; i < dim; i++) {
    112                if (dim != covariances[i].length) {
    113                    throw new DimensionMismatchException(covariances[i].length, dim);
    114                }
    115            }
    116
    117            this.means = MathArrays.copyOf(means);
    118
    119            covarianceMatrix = new Array2DRowRealMatrix(covariances);
    120
    121            // Covariance matrix eigen decomposition.
    122            final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix);
    123
    124            // Compute and store the inverse.
    125            covarianceMatrixInverse = covMatDec.getSolver().getInverse();
    126            // Compute and store the determinant.
    127            covarianceMatrixDeterminant = covMatDec.getDeterminant();
    128
    129            // Eigenvalues of the covariance matrix.
    130            final double[] covMatEigenvalues = covMatDec.getRealEigenvalues();
    131
    132            for (int i = 0; i < covMatEigenvalues.length; i++) {
    133                if (covMatEigenvalues[i] < 0) {
    134                    throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0);
    135                }
    136            }
    137
    138            // Matrix where each column is an eigenvector of the covariance matrix.
    139            final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim);
    140            for (int v = 0; v < dim; v++) {
    141                final double[] evec = covMatDec.getEigenvector(v).toArray();
    142                covMatEigenvectors.setColumn(v, evec);
    143            }
    144
    145            final RealMatrix tmpMatrix = covMatEigenvectors.transpose();
    146
    147            // Scale each eigenvector by the square root of its eigenvalue.
    148            for (int row = 0; row < dim; row++) {
    149                final double factor = FastMath.sqrt(covMatEigenvalues[row]);
    150                for (int col = 0; col < dim; col++) {
    151                    tmpMatrix.multiplyEntry(row, col, factor);
    152                }
    153            }
    154
    155            samplingMatrix = covMatEigenvectors.multiply(tmpMatrix);
    156        }
    157
    158        /**
    159         * Gets the mean vector.
    160         *
    161         * @return the mean vector.
    162         */
    163        public double[] getMeans() {
    164            return MathArrays.copyOf(means);
    165        }
    166
    167        /**
    168         * Gets the covariance matrix.
    169         *
    170         * @return the covariance matrix.
    171         */
    172        public RealMatrix getCovariances() {
    173            return covarianceMatrix.copy();
    174        }
    175
    176        /** {@inheritDoc} */
    177        public double density(final double[] vals) throws DimensionMismatchException {
    178            final int dim = getDimension();
    179            if (vals.length != dim) {
    180                throw new DimensionMismatchException(vals.length, dim);
    181            }
    182
    183            return FastMath.pow(2 * FastMath.PI, -dim / 2) *
    184                FastMath.pow(covarianceMatrixDeterminant, -0.5) *
    185                getExponentTerm(vals);
    186        }
    187
    188        /**
    189         * Gets the square root of each element on the diagonal of the covariance
    190         * matrix.
    191         *
    192         * @return the standard deviations.
    193         */
    194        public double[] getStandardDeviations() {
    195            final int dim = getDimension();
    196            final double[] std = new double[dim];
    197            final double[][] s = covarianceMatrix.getData();
    198            for (int i = 0; i < dim; i++) {
    199                std[i] = FastMath.sqrt(s[i][i]);
    200            }
    201            return std;
    202        }
    203
    204        /** {@inheritDoc} */
    205        public double[] sample() {
    206            final int dim = getDimension();
    207            final double[] normalVals = new double[dim];
    208
    209            for (int i = 0; i < dim; i++) {
    210                normalVals[i] = random.nextGaussian();
    211            }
    212
    213            final double[] vals = samplingMatrix.operate(normalVals);
    214
    215            for (int i = 0; i < dim; i++) {
    216                vals[i] += means[i];
    217            }
    218
    219            return vals;
    220        }
    221
    222        /**
    223         * Computes the term used in the exponent (see definition of the distribution).
    224         *
    225         * @param values Values at which to compute density.
    226         * @return the multiplication factor of density calculations.
    227         */
    228        private double getExponentTerm(final double[] values) {
    229            final double[] centered = new double[values.length];
    230            for (int i = 0; i < centered.length; i++) {
    231                centered[i] = values[i] - getMeans()[i];
    232            }
    233            final double[] preMultiplied = covarianceMatrixInverse.preMultiply(centered);
    234            double sum = 0;
    235            for (int i = 0; i < preMultiplied.length; i++) {
    236                sum += preMultiplied[i] * centered[i];
    237            }
    238            return FastMath.exp(-0.5 * sum);
    239        }
    240    }