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 */
017package org.apache.commons.math3.distribution;
018
019import org.apache.commons.math3.exception.DimensionMismatchException;
020import org.apache.commons.math3.linear.Array2DRowRealMatrix;
021import org.apache.commons.math3.linear.EigenDecomposition;
022import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
023import org.apache.commons.math3.linear.RealMatrix;
024import org.apache.commons.math3.linear.SingularMatrixException;
025import org.apache.commons.math3.random.RandomGenerator;
026import org.apache.commons.math3.random.Well19937c;
027import org.apache.commons.math3.util.FastMath;
028import 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 * @since 3.1
039 */
040public class MultivariateNormalDistribution
041    extends AbstractMultivariateRealDistribution {
042    /** Vector of means. */
043    private final double[] means;
044    /** Covariance matrix. */
045    private final RealMatrix covarianceMatrix;
046    /** The matrix inverse of the covariance matrix. */
047    private final RealMatrix covarianceMatrixInverse;
048    /** The determinant of the covariance matrix. */
049    private final double covarianceMatrixDeterminant;
050    /** Matrix used in computation of samples. */
051    private final RealMatrix samplingMatrix;
052
053    /**
054     * Creates a multivariate normal distribution with the given mean vector and
055     * covariance matrix.
056     * <br/>
057     * The number of dimensions is equal to the length of the mean vector
058     * and to the number of rows and columns of the covariance matrix.
059     * It is frequently written as "p" in formulae.
060     *
061     * @param means Vector of means.
062     * @param covariances Covariance matrix.
063     * @throws DimensionMismatchException if the arrays length are
064     * inconsistent.
065     * @throws SingularMatrixException if the eigenvalue decomposition cannot
066     * be performed on the provided covariance matrix.
067     * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is
068     * negative.
069     */
070    public MultivariateNormalDistribution(final double[] means,
071                                          final double[][] covariances)
072        throws SingularMatrixException,
073               DimensionMismatchException,
074               NonPositiveDefiniteMatrixException {
075        this(new Well19937c(), means, covariances);
076    }
077
078    /**
079     * Creates a multivariate normal distribution with the given mean vector and
080     * covariance matrix.
081     * <br/>
082     * The number of dimensions is equal to the length of the mean vector
083     * and to the number of rows and columns of the covariance matrix.
084     * It is frequently written as "p" in formulae.
085     *
086     * @param rng Random Number Generator.
087     * @param means Vector of means.
088     * @param covariances Covariance matrix.
089     * @throws DimensionMismatchException if the arrays length are
090     * inconsistent.
091     * @throws SingularMatrixException if the eigenvalue decomposition cannot
092     * be performed on the provided covariance matrix.
093     * @throws NonPositiveDefiniteMatrixException if any of the eigenvalues is
094     * negative.
095     */
096    public MultivariateNormalDistribution(RandomGenerator rng,
097                                          final double[] means,
098                                          final double[][] covariances)
099            throws SingularMatrixException,
100                   DimensionMismatchException,
101                   NonPositiveDefiniteMatrixException {
102        super(rng, means.length);
103
104        final int dim = means.length;
105
106        if (covariances.length != dim) {
107            throw new DimensionMismatchException(covariances.length, dim);
108        }
109
110        for (int i = 0; i < dim; i++) {
111            if (dim != covariances[i].length) {
112                throw new DimensionMismatchException(covariances[i].length, dim);
113            }
114        }
115
116        this.means = MathArrays.copyOf(means);
117
118        covarianceMatrix = new Array2DRowRealMatrix(covariances);
119
120        // Covariance matrix eigen decomposition.
121        final EigenDecomposition covMatDec = new EigenDecomposition(covarianceMatrix);
122
123        // Compute and store the inverse.
124        covarianceMatrixInverse = covMatDec.getSolver().getInverse();
125        // Compute and store the determinant.
126        covarianceMatrixDeterminant = covMatDec.getDeterminant();
127
128        // Eigenvalues of the covariance matrix.
129        final double[] covMatEigenvalues = covMatDec.getRealEigenvalues();
130
131        for (int i = 0; i < covMatEigenvalues.length; i++) {
132            if (covMatEigenvalues[i] < 0) {
133                throw new NonPositiveDefiniteMatrixException(covMatEigenvalues[i], i, 0);
134            }
135        }
136
137        // Matrix where each column is an eigenvector of the covariance matrix.
138        final Array2DRowRealMatrix covMatEigenvectors = new Array2DRowRealMatrix(dim, dim);
139        for (int v = 0; v < dim; v++) {
140            final double[] evec = covMatDec.getEigenvector(v).toArray();
141            covMatEigenvectors.setColumn(v, evec);
142        }
143
144        final RealMatrix tmpMatrix = covMatEigenvectors.transpose();
145
146        // Scale each eigenvector by the square root of its eigenvalue.
147        for (int row = 0; row < dim; row++) {
148            final double factor = FastMath.sqrt(covMatEigenvalues[row]);
149            for (int col = 0; col < dim; col++) {
150                tmpMatrix.multiplyEntry(row, col, factor);
151            }
152        }
153
154        samplingMatrix = covMatEigenvectors.multiply(tmpMatrix);
155    }
156
157    /**
158     * Gets the mean vector.
159     *
160     * @return the mean vector.
161     */
162    public double[] getMeans() {
163        return MathArrays.copyOf(means);
164    }
165
166    /**
167     * Gets the covariance matrix.
168     *
169     * @return the covariance matrix.
170     */
171    public RealMatrix getCovariances() {
172        return covarianceMatrix.copy();
173    }
174
175    /** {@inheritDoc} */
176    public double density(final double[] vals) throws DimensionMismatchException {
177        final int dim = getDimension();
178        if (vals.length != dim) {
179            throw new DimensionMismatchException(vals.length, dim);
180        }
181
182        return FastMath.pow(2 * FastMath.PI, -0.5 * dim) *
183            FastMath.pow(covarianceMatrixDeterminant, -0.5) *
184            getExponentTerm(vals);
185    }
186
187    /**
188     * Gets the square root of each element on the diagonal of the covariance
189     * matrix.
190     *
191     * @return the standard deviations.
192     */
193    public double[] getStandardDeviations() {
194        final int dim = getDimension();
195        final double[] std = new double[dim];
196        final double[][] s = covarianceMatrix.getData();
197        for (int i = 0; i < dim; i++) {
198            std[i] = FastMath.sqrt(s[i][i]);
199        }
200        return std;
201    }
202
203    /** {@inheritDoc} */
204    @Override
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}