View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.random;
19  
20  import org.apache.commons.math.DimensionMismatchException;
21  import org.apache.commons.math.linear.RealMatrix;
22  import org.apache.commons.math.linear.RealMatrixImpl;
23  
24  /** 
25   * A {@link RandomVectorGenerator} that generates vectors with with 
26   * correlated components.
27   * <p>Random vectors with correlated components are built by combining
28   * the uncorrelated components of another random vector in such a way that
29   * the resulting correlations are the ones specified by a positive
30   * definite covariance matrix.</p>
31   * <p>The main use for correlated random vector generation is for Monte-Carlo
32   * simulation of physical problems with several variables, for example to
33   * generate error vectors to be added to a nominal vector. A particularly
34   * interesting case is when the generated vector should be drawn from a <a
35   * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
36   * Multivariate Normal Distribution</a>. The approach using a Cholesky
37   * decomposition is quite usual in this case. However, it cas be extended
38   * to other cases as long as the underlying random generator provides
39   * {@link NormalizedRandomGenerator normalized values} like {@link
40   * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
41   * <p>Sometimes, the covariance matrix for a given simulation is not
42   * strictly positive definite. This means that the correlations are
43   * not all independent from each other. In this case, however, the non
44   * strictly positive elements found during the Cholesky decomposition
45   * of the covariance matrix should not be negative either, they
46   * should be null. Another non-conventional extension handling this case
47   * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
48   * where <code>C</code> is the covariance matrix and <code>U</code>
49   * is an uppertriangular matrix, we compute <code>C = B.B<sup>T</sup></code>
50   * where <code>B</code> is a rectangular matrix having
51   * more rows than columns. The number of columns of <code>B</code> is
52   * the rank of the covariance matrix, and it is the dimension of the
53   * uncorrelated random vector that is needed to compute the component
54   * of the correlated vector. This class handles this situation
55   * automatically.</p>
56   *
57   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
58   * @since 1.2
59   */
60  
61  public class CorrelatedRandomVectorGenerator
62  implements RandomVectorGenerator {
63  
64      /** Simple constructor.
65       * <p>Build a correlated random vector generator from its mean
66       * vector and covariance matrix.</p>
67       * @param mean expected mean values for all components
68       * @param covariance covariance matrix
69       * @param small diagonal elements threshold under which  column are
70       * considered to be dependent on previous ones and are discarded
71       * @param generator underlying generator for uncorrelated normalized
72       * components
73       * @exception IllegalArgumentException if there is a dimension
74       * mismatch between the mean vector and the covariance matrix
75       * @exception NotPositiveDefiniteMatrixException if the
76       * covariance matrix is not strictly positive definite
77       * @exception DimensionMismatchException if the mean and covariance
78       * arrays dimensions don't match
79       */
80      public CorrelatedRandomVectorGenerator(double[] mean,
81                                             RealMatrix covariance, double small,
82                                             NormalizedRandomGenerator generator)
83      throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
84  
85          int order = covariance.getRowDimension();
86          if (mean.length != order) {
87              throw new DimensionMismatchException(mean.length, order);
88          }
89          this.mean = (double[]) mean.clone();
90  
91          decompose(covariance, small);
92  
93          this.generator = generator;
94          normalized = new double[rank];
95  
96      }
97  
98      /** Simple constructor.
99       * <p>Build a null mean random correlated vector generator from its
100      * covariance matrix.</p>
101      * @param covariance covariance matrix
102      * @param small diagonal elements threshold under which  column are
103      * considered to be dependent on previous ones and are discarded
104      * @param generator underlying generator for uncorrelated normalized
105      * components
106      * @exception NotPositiveDefiniteMatrixException if the
107      * covariance matrix is not strictly positive definite
108      */
109     public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
110                                            NormalizedRandomGenerator generator)
111     throws NotPositiveDefiniteMatrixException {
112 
113         int order = covariance.getRowDimension();
114         mean = new double[order];
115         for (int i = 0; i < order; ++i) {
116             mean[i] = 0;
117         }
118 
119         decompose(covariance, small);
120 
121         this.generator = generator;
122         normalized = new double[rank];
123 
124     }
125 
126     /** Get the underlying normalized components generator.
127      * @return underlying uncorrelated components generator
128      */
129     public NormalizedRandomGenerator getGenerator() {
130         return generator;
131     }
132 
133     /** Get the root of the covariance matrix.
134      * The root is the rectangular matrix <code>B</code> such that
135      * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
136      * @return root of the square matrix
137      * @see #getRank()
138      */
139     public RealMatrix getRootMatrix() {
140         return root;
141     }
142 
143     /** Get the rank of the covariance matrix.
144      * The rank is the number of independent rows in the covariance
145      * matrix, it is also the number of columns of the rectangular
146      * matrix of the decomposition.
147      * @return rank of the square matrix.
148      * @see #getRootMatrix()
149      */
150     public int getRank() {
151         return rank;
152     }
153 
154     /** Decompose the original square matrix.
155      * <p>The decomposition is based on a Choleski decomposition
156      * where additional transforms are performed:
157      * <ul>
158      *   <li>the rows of the decomposed matrix are permuted</li>
159      *   <li>columns with the too small diagonal element are discarded</li>
160      *   <li>the matrix is permuted</li>
161      * </ul>
162      * This means that rather than computing M = U<sup>T</sup>.U where U
163      * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
164      * where B is a rectangular matrix.
165      * @param covariance covariance matrix
166      * @param small diagonal elements threshold under which  column are
167      * considered to be dependent on previous ones and are discarded
168      * @exception NotPositiveDefiniteMatrixException if the
169      * covariance matrix is not strictly positive definite
170      */
171     private void decompose(RealMatrix covariance, double small)
172     throws NotPositiveDefiniteMatrixException {
173 
174         int order = covariance.getRowDimension();
175         double[][] c = covariance.getData();
176         double[][] b = new double[order][order];
177 
178         int[] swap  = new int[order];
179         int[] index = new int[order];
180         for (int i = 0; i < order; ++i) {
181             index[i] = i;
182         }
183 
184         rank = 0;
185         for (boolean loop = true; loop;) {
186 
187             // find maximal diagonal element
188             swap[rank] = rank;
189             for (int i = rank + 1; i < order; ++i) {
190                 int ii  = index[i];
191                 int isi = index[swap[i]];
192                 if (c[ii][ii] > c[isi][isi]) {
193                     swap[rank] = i;
194                 }
195             }
196 
197 
198             // swap elements
199             if (swap[rank] != rank) {
200                 int tmp = index[rank];
201                 index[rank] = index[swap[rank]];
202                 index[swap[rank]] = tmp;
203             }
204 
205             // check diagonal element
206             int ir = index[rank];
207             if (c[ir][ir] < small) {
208 
209                 if (rank == 0) {
210                     throw new NotPositiveDefiniteMatrixException();
211                 }
212 
213                 // check remaining diagonal elements
214                 for (int i = rank; i < order; ++i) {
215                     if (c[index[i]][index[i]] < -small) {
216                         // there is at least one sufficiently negative diagonal element,
217                         // the covariance matrix is wrong
218                         throw new NotPositiveDefiniteMatrixException();
219                     }
220                 }
221 
222                 // all remaining diagonal elements are close to zero,
223                 // we consider we have found the rank of the covariance matrix
224                 ++rank;
225                 loop = false;
226 
227             } else {
228 
229                 // transform the matrix
230                 double sqrt = Math.sqrt(c[ir][ir]);
231                 b[rank][rank] = sqrt;
232                 double inverse = 1 / sqrt;
233                 for (int i = rank + 1; i < order; ++i) {
234                     int ii = index[i];
235                     double e = inverse * c[ii][ir];
236                     b[i][rank] = e;
237                     c[ii][ii] -= e * e;
238                     for (int j = rank + 1; j < i; ++j) {
239                         int ij = index[j];
240                         double f = c[ii][ij] - e * b[j][rank];
241                         c[ii][ij] = f;
242                         c[ij][ii] = f;
243                     }
244                 }
245 
246                 // prepare next iteration
247                 loop = ++rank < order;
248 
249             }
250 
251         }
252 
253         // build the root matrix
254         root = new RealMatrixImpl(order, rank);
255         for (int i = 0; i < order; ++i) {
256             System.arraycopy(b[i], 0, root.getDataRef()[swap[i]], 0, rank);
257         }
258 
259     }
260 
261     /** Generate a correlated random vector.
262      * @return a random vector as an array of double. The returned array
263      * is created at each call, the caller can do what it wants with it.
264      */
265     public double[] nextVector() {
266 
267         // generate uncorrelated vector
268         for (int i = 0; i < rank; ++i) {
269             normalized[i] = generator.nextNormalizedDouble();
270         }
271 
272         // compute correlated vector
273         double[] correlated = new double[mean.length];
274         for (int i = 0; i < correlated.length; ++i) {
275             correlated[i] = mean[i];
276             for (int j = 0; j < rank; ++j) {
277                 correlated[i] += root.getEntry(i, j) * normalized[j];
278             }
279         }
280 
281         return correlated;
282 
283     }
284 
285     /** Mean vector. */
286     private double[] mean;
287 
288     /** Permutated Cholesky root of the covariance matrix. */
289     private RealMatrixImpl root;
290 
291     /** Rank of the covariance matrix. */
292     private int rank;
293 
294     /** Underlying generator. */
295     private NormalizedRandomGenerator generator;
296 
297     /** Storage for the normalized vector. */
298     private double[] normalized;
299 
300 }