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