    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.math.random;
    019
    020    import org.apache.commons.math.DimensionMismatchException;
    021    import org.apache.commons.math.linear.MatrixUtils;
    022    import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
    023    import org.apache.commons.math.linear.RealMatrix;
    024    import org.apache.commons.math.util.FastMath;
    025
    026    /**
    027     * A {@link RandomVectorGenerator} that generates vectors with with
    028     * correlated components.
    029     * <p>Random vectors with correlated components are built by combining
    030     * the uncorrelated components of another random vector in such a way that
    031     * the resulting correlations are the ones specified by a positive
    032     * definite covariance matrix.</p>
    033     * <p>The main use for correlated random vector generation is for Monte-Carlo
    034     * simulation of physical problems with several variables, for example to
    035     * generate error vectors to be added to a nominal vector. A particularly
    036     * interesting case is when the generated vector should be drawn from a <a
    037     * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
    038     * Multivariate Normal Distribution</a>. The approach using a Cholesky
    039     * decomposition is quite usual in this case. However, it can be extended
    040     * to other cases as long as the underlying random generator provides
    041     * {@link NormalizedRandomGenerator normalized values} like {@link
    042     * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
    043     * <p>Sometimes, the covariance matrix for a given simulation is not
    044     * strictly positive definite. This means that the correlations are
    045     * not all independent from each other. In this case, however, the non
    046     * strictly positive elements found during the Cholesky decomposition
    047     * of the covariance matrix should not be negative either, they
    048     * should be null. Another non-conventional extension handling this case
    049     * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
    050     * where <code>C</code> is the covariance matrix and <code>U</code>
    051     * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code>
    052     * where <code>B</code> is a rectangular matrix having
    053     * more rows than columns. The number of columns of <code>B</code> is
    054     * the rank of the covariance matrix, and it is the dimension of the
    055     * uncorrelated random vector that is needed to compute the component
    056     * of the correlated vector. This class handles this situation
    057     * automatically.</p>
    058     *
    059     * @version $Revision: 1043908 $ $Date: 2010-12-09 12:53:14 +0100 (jeu. 09 d??c. 2010) $
    060     * @since 1.2
    061     */
    062
    063    public class CorrelatedRandomVectorGenerator
    064        implements RandomVectorGenerator {
    065
    066        /** Mean vector. */
    067        private final double[] mean;
    068
    069        /** Underlying generator. */
    070        private final NormalizedRandomGenerator generator;
    071
    072        /** Storage for the normalized vector. */
    073        private final double[] normalized;
    074
    075        /** Permutated Cholesky root of the covariance matrix. */
    076        private RealMatrix root;
    077
    078        /** Rank of the covariance matrix. */
    079        private int rank;
    080
    081        /** Simple constructor.
    082         * <p>Build a correlated random vector generator from its mean
    083         * vector and covariance matrix.</p>
    084         * @param mean expected mean values for all components
    085         * @param covariance covariance matrix
    086         * @param small diagonal elements threshold under which  column are
    087         * considered to be dependent on previous ones and are discarded
    088         * @param generator underlying generator for uncorrelated normalized
    089         * components
    090         * @exception IllegalArgumentException if there is a dimension
    091         * mismatch between the mean vector and the covariance matrix
    092         * @exception NotPositiveDefiniteMatrixException if the
    093         * covariance matrix is not strictly positive definite
    094         * @exception DimensionMismatchException if the mean and covariance
    095         * arrays dimensions don't match
    096         */
    097        public CorrelatedRandomVectorGenerator(double[] mean,
    098                                               RealMatrix covariance, double small,
    099                                               NormalizedRandomGenerator generator)
    100        throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
    101
    102            int order = covariance.getRowDimension();
    103            if (mean.length != order) {
    104                throw new DimensionMismatchException(mean.length, order);
    105            }
    106            this.mean = mean.clone();
    107
    108            decompose(covariance, small);
    109
    110            this.generator = generator;
    111            normalized = new double[rank];
    112
    113        }
    114
    115        /** Simple constructor.
    116         * <p>Build a null mean random correlated vector generator from its
    117         * covariance matrix.</p>
    118         * @param covariance covariance matrix
    119         * @param small diagonal elements threshold under which  column are
    120         * considered to be dependent on previous ones and are discarded
    121         * @param generator underlying generator for uncorrelated normalized
    122         * components
    123         * @exception NotPositiveDefiniteMatrixException if the
    124         * covariance matrix is not strictly positive definite
    125         */
    126        public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
    127                                               NormalizedRandomGenerator generator)
    128        throws NotPositiveDefiniteMatrixException {
    129
    130            int order = covariance.getRowDimension();
    131            mean = new double[order];
    132            for (int i = 0; i < order; ++i) {
    133                mean[i] = 0;
    134            }
    135
    136            decompose(covariance, small);
    137
    138            this.generator = generator;
    139            normalized = new double[rank];
    140
    141        }
    142
    143        /** Get the underlying normalized components generator.
    144         * @return underlying uncorrelated components generator
    145         */
    146        public NormalizedRandomGenerator getGenerator() {
    147            return generator;
    148        }
    149
    150        /** Get the root of the covariance matrix.
    151         * The root is the rectangular matrix <code>B</code> such that
    152         * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
    153         * @return root of the square matrix
    154         * @see #getRank()
    155         */
    156        public RealMatrix getRootMatrix() {
    157            return root;
    158        }
    159
    160        /** Get the rank of the covariance matrix.
    161         * The rank is the number of independent rows in the covariance
    162         * matrix, it is also the number of columns of the rectangular
    163         * matrix of the decomposition.
    164         * @return rank of the square matrix.
    165         * @see #getRootMatrix()
    166         */
    167        public int getRank() {
    168            return rank;
    169        }
    170
    171        /** Decompose the original square matrix.
    172         * <p>The decomposition is based on a Choleski decomposition
    173         * where additional transforms are performed:
    174         * <ul>
    175         *   <li>the rows of the decomposed matrix are permuted</li>
    176         *   <li>columns with the too small diagonal element are discarded</li>
    177         *   <li>the matrix is permuted</li>
    178         * </ul>
    179         * This means that rather than computing M = U<sup>T</sup>.U where U
    180         * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
    181         * where B is a rectangular matrix.
    182         * @param covariance covariance matrix
    183         * @param small diagonal elements threshold under which  column are
    184         * considered to be dependent on previous ones and are discarded
    185         * @exception NotPositiveDefiniteMatrixException if the
    186         * covariance matrix is not strictly positive definite
    187         */
    188        private void decompose(RealMatrix covariance, double small)
    189        throws NotPositiveDefiniteMatrixException {
    190
    191            int order = covariance.getRowDimension();
    192            double[][] c = covariance.getData();
    193            double[][] b = new double[order][order];
    194
    195            int[] swap  = new int[order];
    196            int[] index = new int[order];
    197            for (int i = 0; i < order; ++i) {
    198                index[i] = i;
    199            }
    200
    201            rank = 0;
    202            for (boolean loop = true; loop;) {
    203
    204                // find maximal diagonal element
    205                swap[rank] = rank;
    206                for (int i = rank + 1; i < order; ++i) {
    207                    int ii  = index[i];
    208                    int isi = index[swap[i]];
    209                    if (c[ii][ii] > c[isi][isi]) {
    210                        swap[rank] = i;
    211                    }
    212                }
    213
    214
    215                // swap elements
    216                if (swap[rank] != rank) {
    217                    int tmp = index[rank];
    218                    index[rank] = index[swap[rank]];
    219                    index[swap[rank]] = tmp;
    220                }
    221
    222                // check diagonal element
    223                int ir = index[rank];
    224                if (c[ir][ir] < small) {
    225
    226                    if (rank == 0) {
    227                        throw new NotPositiveDefiniteMatrixException();
    228                    }
    229
    230                    // check remaining diagonal elements
    231                    for (int i = rank; i < order; ++i) {
    232                        if (c[index[i]][index[i]] < -small) {
    233                            // there is at least one sufficiently negative diagonal element,
    234                            // the covariance matrix is wrong
    235                            throw new NotPositiveDefiniteMatrixException();
    236                        }
    237                    }
    238
    239                    // all remaining diagonal elements are close to zero,
    240                    // we consider we have found the rank of the covariance matrix
    241                    ++rank;
    242                    loop = false;
    243
    244                } else {
    245
    246                    // transform the matrix
    247                    double sqrt = FastMath.sqrt(c[ir][ir]);
    248                    b[rank][rank] = sqrt;
    249                    double inverse = 1 / sqrt;
    250                    for (int i = rank + 1; i < order; ++i) {
    251                        int ii = index[i];
    252                        double e = inverse * c[ii][ir];
    253                        b[i][rank] = e;
    254                        c[ii][ii] -= e * e;
    255                        for (int j = rank + 1; j < i; ++j) {
    256                            int ij = index[j];
    257                            double f = c[ii][ij] - e * b[j][rank];
    258                            c[ii][ij] = f;
    259                            c[ij][ii] = f;
    260                        }
    261                    }
    262
    263                    // prepare next iteration
    264                    loop = ++rank < order;
    265
    266                }
    267
    268            }
    269
    270            // build the root matrix
    271            root = MatrixUtils.createRealMatrix(order, rank);
    272            for (int i = 0; i < order; ++i) {
    273                for (int j = 0; j < rank; ++j) {
    274                    root.setEntry(index[i], j, b[i][j]);
    275                }
    276            }
    277
    278        }
    279
    280        /** Generate a correlated random vector.
    281         * @return a random vector as an array of double. The returned array
    282         * is created at each call, the caller can do what it wants with it.
    283         */
    284        public double[] nextVector() {
    285
    286            // generate uncorrelated vector
    287            for (int i = 0; i < rank; ++i) {
    288                normalized[i] = generator.nextNormalizedDouble();
    289            }
    290
    291            // compute correlated vector
    292            double[] correlated = new double[mean.length];
    293            for (int i = 0; i < correlated.length; ++i) {
    294                correlated[i] = mean[i];
    295                for (int j = 0; j < rank; ++j) {
    296                    correlated[i] += root.getEntry(i, j) * normalized[j];
    297                }
    298            }
    299
    300            return correlated;
    301
    302        }
    303
    304    }