Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 déc. 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 | } |