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 }