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.MatrixUtils;
22  import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
23  import org.apache.commons.math.linear.RealMatrix;
24  
25  /** 
26   * A {@link RandomVectorGenerator} that generates vectors with with 
27   * correlated components.
28   * <p>Random vectors with correlated components are built by combining
29   * the uncorrelated components of another random vector in such a way that
30   * the resulting correlations are the ones specified by a positive
31   * definite covariance matrix.</p>
32   * <p>The main use for correlated random vector generation is for Monte-Carlo
33   * simulation of physical problems with several variables, for example to
34   * generate error vectors to be added to a nominal vector. A particularly
35   * interesting case is when the generated vector should be drawn from a <a
36   * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
37   * Multivariate Normal Distribution</a>. The approach using a Cholesky
38   * decomposition is quite usual in this case. However, it cas be extended
39   * to other cases as long as the underlying random generator provides
40   * {@link NormalizedRandomGenerator normalized values} like {@link
41   * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
42   * <p>Sometimes, the covariance matrix for a given simulation is not
43   * strictly positive definite. This means that the correlations are
44   * not all independent from each other. In this case, however, the non
45   * strictly positive elements found during the Cholesky decomposition
46   * of the covariance matrix should not be negative either, they
47   * should be null. Another non-conventional extension handling this case
48   * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
49   * where <code>C</code> is the covariance matrix and <code>U</code>
50   * is an uppertriangular matrix, we compute <code>C = B.B<sup>T</sup></code>
51   * where <code>B</code> is a rectangular matrix having
52   * more rows than columns. The number of columns of <code>B</code> is
53   * the rank of the covariance matrix, and it is the dimension of the
54   * uncorrelated random vector that is needed to compute the component
55   * of the correlated vector. This class handles this situation
56   * automatically.</p>
57   *
58   * @version $Revision: 781122 $ $Date: 2009-06-02 14:53:23 -0400 (Tue, 02 Jun 2009) $
59   * @since 1.2
60   */
61  
62  public class CorrelatedRandomVectorGenerator
63      implements RandomVectorGenerator {
64  
65      /** Simple constructor.
66       * <p>Build a correlated random vector generator from its mean
67       * vector and covariance matrix.</p>
68       * @param mean expected mean values for all components
69       * @param covariance covariance matrix
70       * @param small diagonal elements threshold under which  column are
71       * considered to be dependent on previous ones and are discarded
72       * @param generator underlying generator for uncorrelated normalized
73       * components
74       * @exception IllegalArgumentException if there is a dimension
75       * mismatch between the mean vector and the covariance matrix
76       * @exception NotPositiveDefiniteMatrixException if the
77       * covariance matrix is not strictly positive definite
78       * @exception DimensionMismatchException if the mean and covariance
79       * arrays dimensions don't match
80       */
81      public CorrelatedRandomVectorGenerator(double[] mean,
82                                             RealMatrix covariance, double small,
83                                             NormalizedRandomGenerator generator)
84      throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
85  
86          int order = covariance.getRowDimension();
87          if (mean.length != order) {
88              throw new DimensionMismatchException(mean.length, order);
89          }
90          this.mean = mean.clone();
91  
92          decompose(covariance, small);
93  
94          this.generator = generator;
95          normalized = new double[rank];
96  
97      }
98  
99      /** 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 }