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 }