/*
* 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;
}
}