    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
    018    package org.apache.commons.math3.random;
    019
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.linear.RealMatrix;
    022    import org.apache.commons.math3.linear.RectangularCholeskyDecomposition;
    023
    024    /**
    025     * A {@link RandomVectorGenerator} that generates vectors with with
    026     * correlated components.
    027     * <p>Random vectors with correlated components are built by combining
    028     * the uncorrelated components of another random vector in such a way that
    029     * the resulting correlations are the ones specified by a positive
    030     * definite covariance matrix.</p>
    031     * <p>The main use for correlated random vector generation is for Monte-Carlo
    032     * simulation of physical problems with several variables, for example to
    033     * generate error vectors to be added to a nominal vector. A particularly
    034     * interesting case is when the generated vector should be drawn from a <a
    035     * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
    036     * Multivariate Normal Distribution</a>. The approach using a Cholesky
    037     * decomposition is quite usual in this case. However, it can be extended
    038     * to other cases as long as the underlying random generator provides
    039     * {@link NormalizedRandomGenerator normalized values} like {@link
    040     * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
    041     * <p>Sometimes, the covariance matrix for a given simulation is not
    042     * strictly positive definite. This means that the correlations are
    043     * not all independent from each other. In this case, however, the non
    044     * strictly positive elements found during the Cholesky decomposition
    045     * of the covariance matrix should not be negative either, they
    046     * should be null. Another non-conventional extension handling this case
    047     * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
    048     * where <code>C</code> is the covariance matrix and <code>U</code>
    049     * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code>
    050     * where <code>B</code> is a rectangular matrix having
    051     * more rows than columns. The number of columns of <code>B</code> is
    052     * the rank of the covariance matrix, and it is the dimension of the
    053     * uncorrelated random vector that is needed to compute the component
    054     * of the correlated vector. This class handles this situation
    055     * automatically.</p>
    056     *
    057     * @version $Id: CorrelatedRandomVectorGenerator.java 1416643 2012-12-03 19:37:14Z tn $
    058     * @since 1.2
    059     */
    060
    061    public class CorrelatedRandomVectorGenerator
    062        implements RandomVectorGenerator {
    063        /** Mean vector. */
    064        private final double[] mean;
    065        /** Underlying generator. */
    066        private final NormalizedRandomGenerator generator;
    067        /** Storage for the normalized vector. */
    068        private final double[] normalized;
    069        /** Root of the covariance matrix. */
    070        private final RealMatrix root;
    071
    072        /**
    073         * Builds a correlated random vector generator from its mean
    074         * vector and covariance matrix.
    075         *
    076         * @param mean Expected mean values for all components.
    077         * @param covariance Covariance matrix.
    078         * @param small Diagonal elements threshold under which  column are
    079         * considered to be dependent on previous ones and are discarded
    080         * @param generator underlying generator for uncorrelated normalized
    081         * components.
    082         * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException
    083         * if the covariance matrix is not strictly positive definite.
    084         * @throws DimensionMismatchException if the mean and covariance
    085         * arrays dimensions do not match.
    086         */
    087        public CorrelatedRandomVectorGenerator(double[] mean,
    088                                               RealMatrix covariance, double small,
    089                                               NormalizedRandomGenerator generator) {
    090            int order = covariance.getRowDimension();
    091            if (mean.length != order) {
    092                throw new DimensionMismatchException(mean.length, order);
    093            }
    094            this.mean = mean.clone();
    095
    096            final RectangularCholeskyDecomposition decomposition =
    097                new RectangularCholeskyDecomposition(covariance, small);
    098            root = decomposition.getRootMatrix();
    099
    100            this.generator = generator;
    101            normalized = new double[decomposition.getRank()];
    102
    103        }
    104
    105        /**
    106         * Builds a null mean random correlated vector generator from its
    107         * covariance matrix.
    108         *
    109         * @param covariance Covariance matrix.
    110         * @param small Diagonal elements threshold under which  column are
    111         * considered to be dependent on previous ones and are discarded.
    112         * @param generator Underlying generator for uncorrelated normalized
    113         * components.
    114         * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException
    115         * if the covariance matrix is not strictly positive definite.
    116         */
    117        public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
    118                                               NormalizedRandomGenerator generator) {
    119            int order = covariance.getRowDimension();
    120            mean = new double[order];
    121            for (int i = 0; i < order; ++i) {
    122                mean[i] = 0;
    123            }
    124
    125            final RectangularCholeskyDecomposition decomposition =
    126                new RectangularCholeskyDecomposition(covariance, small);
    127            root = decomposition.getRootMatrix();
    128
    129            this.generator = generator;
    130            normalized = new double[decomposition.getRank()];
    131
    132        }
    133
    134        /** Get the underlying normalized components generator.
    135         * @return underlying uncorrelated components generator
    136         */
    137        public NormalizedRandomGenerator getGenerator() {
    138            return generator;
    139        }
    140
    141        /** Get the rank of the covariance matrix.
    142         * The rank is the number of independent rows in the covariance
    143         * matrix, it is also the number of columns of the root matrix.
    144         * @return rank of the square matrix.
    145         * @see #getRootMatrix()
    146         */
    147        public int getRank() {
    148            return normalized.length;
    149        }
    150
    151        /** Get the root of the covariance matrix.
    152         * The root is the rectangular matrix <code>B</code> such that
    153         * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
    154         * @return root of the square matrix
    155         * @see #getRank()
    156         */
    157        public RealMatrix getRootMatrix() {
    158            return root;
    159        }
    160
    161        /** Generate a correlated random vector.
    162         * @return a random vector as an array of double. The returned array
    163         * is created at each call, the caller can do what it wants with it.
    164         */
    165        public double[] nextVector() {
    166
    167            // generate uncorrelated vector
    168            for (int i = 0; i < normalized.length; ++i) {
    169                normalized[i] = generator.nextNormalizedDouble();
    170            }
    171
    172            // compute correlated vector
    173            double[] correlated = new double[mean.length];
    174            for (int i = 0; i < correlated.length; ++i) {
    175                correlated[i] = mean[i];
    176                for (int j = 0; j < root.getColumnDimension(); ++j) {
    177                    correlated[i] += root.getEntry(i, j) * normalized[j];
    178                }
    179            }
    180
    181            return correlated;
    182
    183        }
    184
    185    }