| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| package org.apache.commons.math.random; |
| |
| import org.apache.commons.math.DimensionMismatchException; |
| import org.apache.commons.math.linear.MatrixUtils; |
| import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException; |
| import org.apache.commons.math.linear.RealMatrix; |
| import org.apache.commons.math.util.FastMath; |
| |
| /** |
| * A {@link RandomVectorGenerator} that generates vectors with with |
| * correlated components. |
| * <p>Random vectors with correlated components are built by combining |
| * the uncorrelated components of another random vector in such a way that |
| * the resulting correlations are the ones specified by a positive |
| * definite covariance matrix.</p> |
| * <p>The main use for correlated random vector generation is for Monte-Carlo |
| * simulation of physical problems with several variables, for example to |
| * generate error vectors to be added to a nominal vector. A particularly |
| * interesting case is when the generated vector should be drawn from a <a |
| * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution"> |
| * Multivariate Normal Distribution</a>. The approach using a Cholesky |
| * decomposition is quite usual in this case. However, it can be extended |
| * to other cases as long as the underlying random generator provides |
| * {@link NormalizedRandomGenerator normalized values} like {@link |
| * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p> |
| * <p>Sometimes, the covariance matrix for a given simulation is not |
| * strictly positive definite. This means that the correlations are |
| * not all independent from each other. In this case, however, the non |
| * strictly positive elements found during the Cholesky decomposition |
| * of the covariance matrix should not be negative either, they |
| * should be null. Another non-conventional extension handling this case |
| * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code> |
| * where <code>C</code> is the covariance matrix and <code>U</code> |
| * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code> |
| * where <code>B</code> is a rectangular matrix having |
| * more rows than columns. The number of columns of <code>B</code> is |
| * the rank of the covariance matrix, and it is the dimension of the |
| * uncorrelated random vector that is needed to compute the component |
| * of the correlated vector. This class handles this situation |
| * automatically.</p> |
| * |
| * @version $Revision: 1043908 $ $Date: 2010-12-09 12:53:14 +0100 (jeu. 09 déc. 2010) $ |
| * @since 1.2 |
| */ |
| |
| public class CorrelatedRandomVectorGenerator |
| implements RandomVectorGenerator { |
| |
| /** Mean vector. */ |
| private final double[] mean; |
| |
| /** Underlying generator. */ |
| private final NormalizedRandomGenerator generator; |
| |
| /** Storage for the normalized vector. */ |
| private final double[] normalized; |
| |
| /** Permutated Cholesky root of the covariance matrix. */ |
| private RealMatrix root; |
| |
| /** Rank of the covariance matrix. */ |
| private int rank; |
| |
| /** Simple constructor. |
| * <p>Build a correlated random vector generator from its mean |
| * vector and covariance matrix.</p> |
| * @param mean expected mean values for all components |
| * @param covariance covariance matrix |
| * @param small diagonal elements threshold under which column are |
| * considered to be dependent on previous ones and are discarded |
| * @param generator underlying generator for uncorrelated normalized |
| * components |
| * @exception IllegalArgumentException if there is a dimension |
| * mismatch between the mean vector and the covariance matrix |
| * @exception NotPositiveDefiniteMatrixException if the |
| * covariance matrix is not strictly positive definite |
| * @exception DimensionMismatchException if the mean and covariance |
| * arrays dimensions don't match |
| */ |
| public CorrelatedRandomVectorGenerator(double[] mean, |
| RealMatrix covariance, double small, |
| NormalizedRandomGenerator generator) |
| throws NotPositiveDefiniteMatrixException, DimensionMismatchException { |
| |
| int order = covariance.getRowDimension(); |
| if (mean.length != order) { |
| throw new DimensionMismatchException(mean.length, order); |
| } |
| this.mean = mean.clone(); |
| |
| decompose(covariance, small); |
| |
| this.generator = generator; |
| normalized = new double[rank]; |
| |
| } |
| |
| /** Simple constructor. |
| * <p>Build a null mean random correlated vector generator from its |
| * covariance matrix.</p> |
| * @param covariance covariance matrix |
| * @param small diagonal elements threshold under which column are |
| * considered to be dependent on previous ones and are discarded |
| * @param generator underlying generator for uncorrelated normalized |
| * components |
| * @exception NotPositiveDefiniteMatrixException if the |
| * covariance matrix is not strictly positive definite |
| */ |
| public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small, |
| NormalizedRandomGenerator generator) |
| throws NotPositiveDefiniteMatrixException { |
| |
| int order = covariance.getRowDimension(); |
| mean = new double[order]; |
| for (int i = 0; i < order; ++i) { |
| mean[i] = 0; |
| } |
| |
| decompose(covariance, small); |
| |
| this.generator = generator; |
| normalized = new double[rank]; |
| |
| } |
| |
| /** Get the underlying normalized components generator. |
| * @return underlying uncorrelated components generator |
| */ |
| public NormalizedRandomGenerator getGenerator() { |
| return generator; |
| } |
| |
| /** Get the root of the covariance matrix. |
| * The root is the rectangular matrix <code>B</code> such that |
| * the covariance matrix is equal to <code>B.B<sup>T</sup></code> |
| * @return root of the square matrix |
| * @see #getRank() |
| */ |
| public RealMatrix getRootMatrix() { |
| return root; |
| } |
| |
| /** Get the rank of the covariance matrix. |
| * The rank is the number of independent rows in the covariance |
| * matrix, it is also the number of columns of the rectangular |
| * matrix of the decomposition. |
| * @return rank of the square matrix. |
| * @see #getRootMatrix() |
| */ |
| public int getRank() { |
| return rank; |
| } |
| |
| /** Decompose the original square matrix. |
| * <p>The decomposition is based on a Choleski decomposition |
| * where additional transforms are performed: |
| * <ul> |
| * <li>the rows of the decomposed matrix are permuted</li> |
| * <li>columns with the too small diagonal element are discarded</li> |
| * <li>the matrix is permuted</li> |
| * </ul> |
| * This means that rather than computing M = U<sup>T</sup>.U where U |
| * is an upper triangular matrix, this method computed M=B.B<sup>T</sup> |
| * where B is a rectangular matrix. |
| * @param covariance covariance matrix |
| * @param small diagonal elements threshold under which column are |
| * considered to be dependent on previous ones and are discarded |
| * @exception NotPositiveDefiniteMatrixException if the |
| * covariance matrix is not strictly positive definite |
| */ |
| private void decompose(RealMatrix covariance, double small) |
| throws NotPositiveDefiniteMatrixException { |
| |
| int order = covariance.getRowDimension(); |
| double[][] c = covariance.getData(); |
| double[][] b = new double[order][order]; |
| |
| int[] swap = new int[order]; |
| int[] index = new int[order]; |
| for (int i = 0; i < order; ++i) { |
| index[i] = i; |
| } |
| |
| rank = 0; |
| for (boolean loop = true; loop;) { |
| |
| // find maximal diagonal element |
| swap[rank] = rank; |
| for (int i = rank + 1; i < order; ++i) { |
| int ii = index[i]; |
| int isi = index[swap[i]]; |
| if (c[ii][ii] > c[isi][isi]) { |
| swap[rank] = i; |
| } |
| } |
| |
| |
| // swap elements |
| if (swap[rank] != rank) { |
| int tmp = index[rank]; |
| index[rank] = index[swap[rank]]; |
| index[swap[rank]] = tmp; |
| } |
| |
| // check diagonal element |
| int ir = index[rank]; |
| if (c[ir][ir] < small) { |
| |
| if (rank == 0) { |
| throw new NotPositiveDefiniteMatrixException(); |
| } |
| |
| // check remaining diagonal elements |
| for (int i = rank; i < order; ++i) { |
| if (c[index[i]][index[i]] < -small) { |
| // there is at least one sufficiently negative diagonal element, |
| // the covariance matrix is wrong |
| throw new NotPositiveDefiniteMatrixException(); |
| } |
| } |
| |
| // all remaining diagonal elements are close to zero, |
| // we consider we have found the rank of the covariance matrix |
| ++rank; |
| loop = false; |
| |
| } else { |
| |
| // transform the matrix |
| double sqrt = FastMath.sqrt(c[ir][ir]); |
| b[rank][rank] = sqrt; |
| double inverse = 1 / sqrt; |
| for (int i = rank + 1; i < order; ++i) { |
| int ii = index[i]; |
| double e = inverse * c[ii][ir]; |
| b[i][rank] = e; |
| c[ii][ii] -= e * e; |
| for (int j = rank + 1; j < i; ++j) { |
| int ij = index[j]; |
| double f = c[ii][ij] - e * b[j][rank]; |
| c[ii][ij] = f; |
| c[ij][ii] = f; |
| } |
| } |
| |
| // prepare next iteration |
| loop = ++rank < order; |
| |
| } |
| |
| } |
| |
| // build the root matrix |
| root = MatrixUtils.createRealMatrix(order, rank); |
| for (int i = 0; i < order; ++i) { |
| for (int j = 0; j < rank; ++j) { |
| root.setEntry(index[i], j, b[i][j]); |
| } |
| } |
| |
| } |
| |
| /** Generate a correlated random vector. |
| * @return a random vector as an array of double. The returned array |
| * is created at each call, the caller can do what it wants with it. |
| */ |
| public double[] nextVector() { |
| |
| // generate uncorrelated vector |
| for (int i = 0; i < rank; ++i) { |
| normalized[i] = generator.nextNormalizedDouble(); |
| } |
| |
| // compute correlated vector |
| double[] correlated = new double[mean.length]; |
| for (int i = 0; i < correlated.length; ++i) { |
| correlated[i] = mean[i]; |
| for (int j = 0; j < rank; ++j) { |
| correlated[i] += root.getEntry(i, j) * normalized[j]; |
| } |
| } |
| |
| return correlated; |
| |
| } |
| |
| } |