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 }