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 }