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: 811827 $ $Date: 2009-09-06 11:32:50 -0400 (Sun, 06 Sep 2009) $
059     * @since 1.2
060     */
061    
062    public class CorrelatedRandomVectorGenerator
063        implements RandomVectorGenerator {
064    
065        /** Mean vector. */
066        private final double[] mean;
067    
068        /** Underlying generator. */
069        private final NormalizedRandomGenerator generator;
070    
071        /** Storage for the normalized vector. */
072        private final double[] normalized;
073    
074        /** Permutated Cholesky root of the covariance matrix. */
075        private RealMatrix root;
076    
077        /** Rank of the covariance matrix. */
078        private int rank;
079    
080        /** Simple constructor.
081         * <p>Build a correlated random vector generator from its mean
082         * vector and covariance matrix.</p>
083         * @param mean expected mean values for all components
084         * @param covariance covariance matrix
085         * @param small diagonal elements threshold under which  column are
086         * considered to be dependent on previous ones and are discarded
087         * @param generator underlying generator for uncorrelated normalized
088         * components
089         * @exception IllegalArgumentException if there is a dimension
090         * mismatch between the mean vector and the covariance matrix
091         * @exception NotPositiveDefiniteMatrixException if the
092         * covariance matrix is not strictly positive definite
093         * @exception DimensionMismatchException if the mean and covariance
094         * arrays dimensions don't match
095         */
096        public CorrelatedRandomVectorGenerator(double[] mean,
097                                               RealMatrix covariance, double small,
098                                               NormalizedRandomGenerator generator)
099        throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
100    
101            int order = covariance.getRowDimension();
102            if (mean.length != order) {
103                throw new DimensionMismatchException(mean.length, order);
104            }
105            this.mean = mean.clone();
106    
107            decompose(covariance, small);
108    
109            this.generator = generator;
110            normalized = new double[rank];
111    
112        }
113    
114        /** Simple constructor.
115         * <p>Build a null mean random correlated vector generator from its
116         * covariance matrix.</p>
117         * @param covariance covariance matrix
118         * @param small diagonal elements threshold under which  column are
119         * considered to be dependent on previous ones and are discarded
120         * @param generator underlying generator for uncorrelated normalized
121         * components
122         * @exception NotPositiveDefiniteMatrixException if the
123         * covariance matrix is not strictly positive definite
124         */
125        public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
126                                               NormalizedRandomGenerator generator)
127        throws NotPositiveDefiniteMatrixException {
128    
129            int order = covariance.getRowDimension();
130            mean = new double[order];
131            for (int i = 0; i < order; ++i) {
132                mean[i] = 0;
133            }
134    
135            decompose(covariance, small);
136    
137            this.generator = generator;
138            normalized = new double[rank];
139    
140        }
141    
142        /** Get the underlying normalized components generator.
143         * @return underlying uncorrelated components generator
144         */
145        public NormalizedRandomGenerator getGenerator() {
146            return generator;
147        }
148    
149        /** Get the root of the covariance matrix.
150         * The root is the rectangular matrix <code>B</code> such that
151         * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
152         * @return root of the square matrix
153         * @see #getRank()
154         */
155        public RealMatrix getRootMatrix() {
156            return root;
157        }
158    
159        /** Get the rank of the covariance matrix.
160         * The rank is the number of independent rows in the covariance
161         * matrix, it is also the number of columns of the rectangular
162         * matrix of the decomposition.
163         * @return rank of the square matrix.
164         * @see #getRootMatrix()
165         */
166        public int getRank() {
167            return rank;
168        }
169    
170        /** Decompose the original square matrix.
171         * <p>The decomposition is based on a Choleski decomposition
172         * where additional transforms are performed:
173         * <ul>
174         *   <li>the rows of the decomposed matrix are permuted</li>
175         *   <li>columns with the too small diagonal element are discarded</li>
176         *   <li>the matrix is permuted</li>
177         * </ul>
178         * This means that rather than computing M = U<sup>T</sup>.U where U
179         * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
180         * where B is a rectangular matrix.
181         * @param covariance covariance matrix
182         * @param small diagonal elements threshold under which  column are
183         * considered to be dependent on previous ones and are discarded
184         * @exception NotPositiveDefiniteMatrixException if the
185         * covariance matrix is not strictly positive definite
186         */
187        private void decompose(RealMatrix covariance, double small)
188        throws NotPositiveDefiniteMatrixException {
189    
190            int order = covariance.getRowDimension();
191            double[][] c = covariance.getData();
192            double[][] b = new double[order][order];
193    
194            int[] swap  = new int[order];
195            int[] index = new int[order];
196            for (int i = 0; i < order; ++i) {
197                index[i] = i;
198            }
199    
200            rank = 0;
201            for (boolean loop = true; loop;) {
202    
203                // find maximal diagonal element
204                swap[rank] = rank;
205                for (int i = rank + 1; i < order; ++i) {
206                    int ii  = index[i];
207                    int isi = index[swap[i]];
208                    if (c[ii][ii] > c[isi][isi]) {
209                        swap[rank] = i;
210                    }
211                }
212    
213    
214                // swap elements
215                if (swap[rank] != rank) {
216                    int tmp = index[rank];
217                    index[rank] = index[swap[rank]];
218                    index[swap[rank]] = tmp;
219                }
220    
221                // check diagonal element
222                int ir = index[rank];
223                if (c[ir][ir] < small) {
224    
225                    if (rank == 0) {
226                        throw new NotPositiveDefiniteMatrixException();
227                    }
228    
229                    // check remaining diagonal elements
230                    for (int i = rank; i < order; ++i) {
231                        if (c[index[i]][index[i]] < -small) {
232                            // there is at least one sufficiently negative diagonal element,
233                            // the covariance matrix is wrong
234                            throw new NotPositiveDefiniteMatrixException();
235                        }
236                    }
237    
238                    // all remaining diagonal elements are close to zero,
239                    // we consider we have found the rank of the covariance matrix
240                    ++rank;
241                    loop = false;
242    
243                } else {
244    
245                    // transform the matrix
246                    double sqrt = Math.sqrt(c[ir][ir]);
247                    b[rank][rank] = sqrt;
248                    double inverse = 1 / sqrt;
249                    for (int i = rank + 1; i < order; ++i) {
250                        int ii = index[i];
251                        double e = inverse * c[ii][ir];
252                        b[i][rank] = e;
253                        c[ii][ii] -= e * e;
254                        for (int j = rank + 1; j < i; ++j) {
255                            int ij = index[j];
256                            double f = c[ii][ij] - e * b[j][rank];
257                            c[ii][ij] = f;
258                            c[ij][ii] = f;
259                        }
260                    }
261    
262                    // prepare next iteration
263                    loop = ++rank < order;
264    
265                }
266    
267            }
268    
269            // build the root matrix
270            root = MatrixUtils.createRealMatrix(order, rank);
271            for (int i = 0; i < order; ++i) {
272                for (int j = 0; j < rank; ++j) {
273                    root.setEntry(index[i], j, b[i][j]);
274                }
275            }
276    
277        }
278    
279        /** Generate a correlated random vector.
280         * @return a random vector as an array of double. The returned array
281         * is created at each call, the caller can do what it wants with it.
282         */
283        public double[] nextVector() {
284    
285            // generate uncorrelated vector
286            for (int i = 0; i < rank; ++i) {
287                normalized[i] = generator.nextNormalizedDouble();
288            }
289    
290            // compute correlated vector
291            double[] correlated = new double[mean.length];
292            for (int i = 0; i < correlated.length; ++i) {
293                correlated[i] = mean[i];
294                for (int j = 0; j < rank; ++j) {
295                    correlated[i] += root.getEntry(i, j) * normalized[j];
296                }
297            }
298    
299            return correlated;
300    
301        }
302    
303    }