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    }