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    }