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