001/*002 * Licensed to the Apache Software Foundation (ASF) under one or more003 * contributor license agreements. See the NOTICE file distributed with004 * this work for additional information regarding copyright ownership.005 * The ASF licenses this file to You under the Apache License, Version 2.0006 * (the "License"); you may not use this file except in compliance with007 * the License. You may obtain a copy of the License at008 *009 * http://www.apache.org/licenses/LICENSE-2.0010 *011 * Unless required by applicable law or agreed to in writing, software012 * 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 and015 * limitations under the License.016 */017018package org.apache.commons.math3.random;019020import org.apache.commons.math3.exception.DimensionMismatchException;021import org.apache.commons.math3.linear.RealMatrix;022import org.apache.commons.math3.linear.RectangularCholeskyDecomposition;023024/**025 * A {@link RandomVectorGenerator} that generates vectors with with026 * correlated components.027 * <p>Random vectors with correlated components are built by combining028 * the uncorrelated components of another random vector in such a way that029 * the resulting correlations are the ones specified by a positive030 * definite covariance matrix.</p>031 * <p>The main use for correlated random vector generation is for Monte-Carlo032 * simulation of physical problems with several variables, for example to033 * generate error vectors to be added to a nominal vector. A particularly034 * interesting case is when the generated vector should be drawn from a <a035 * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">036 * Multivariate Normal Distribution</a>. The approach using a Cholesky037 * decomposition is quite usual in this case. However, it can be extended038 * to other cases as long as the underlying random generator provides039 * {@link NormalizedRandomGenerator normalized values} like {@link040 * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>041 * <p>Sometimes, the covariance matrix for a given simulation is not042 * strictly positive definite. This means that the correlations are043 * not all independent from each other. In this case, however, the non044 * strictly positive elements found during the Cholesky decomposition045 * of the covariance matrix should not be negative either, they046 * should be null. Another non-conventional extension handling this case047 * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>048 * where <code>C</code> is the covariance matrix and <code>U</code>049 * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code>050 * where <code>B</code> is a rectangular matrix having051 * more rows than columns. The number of columns of <code>B</code> is052 * the rank of the covariance matrix, and it is the dimension of the053 * uncorrelated random vector that is needed to compute the component054 * of the correlated vector. This class handles this situation055 * automatically.</p>056 *057 * @since 1.2058 */059060public class CorrelatedRandomVectorGenerator061 implements RandomVectorGenerator {062 /** Mean vector. */063 private final double[] mean;064 /** Underlying generator. */065 private final NormalizedRandomGenerator generator;066 /** Storage for the normalized vector. */067 private final double[] normalized;068 /** Root of the covariance matrix. */069 private final RealMatrix root;070071 /**072 * Builds a correlated random vector generator from its mean073 * vector and covariance matrix.074 *075 * @param mean Expected mean values for all components.076 * @param covariance Covariance matrix.077 * @param small Diagonal elements threshold under which column are078 * considered to be dependent on previous ones and are discarded079 * @param generator underlying generator for uncorrelated normalized080 * components.081 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException082 * if the covariance matrix is not strictly positive definite.083 * @throws DimensionMismatchException if the mean and covariance084 * arrays dimensions do not match.085 */086 public CorrelatedRandomVectorGenerator(double[] mean,087 RealMatrix covariance, double small,088 NormalizedRandomGenerator generator) {089 int order = covariance.getRowDimension();090 if (mean.length != order) {091 throw new DimensionMismatchException(mean.length, order);092 }093 this.mean = mean.clone();094095 final RectangularCholeskyDecomposition decomposition =096 new RectangularCholeskyDecomposition(covariance, small);097 root = decomposition.getRootMatrix();098099 this.generator = generator;100 normalized = new double[decomposition.getRank()];101102 }103104 /**105 * Builds a null mean random correlated vector generator from its106 * covariance matrix.107 *108 * @param covariance Covariance matrix.109 * @param small Diagonal elements threshold under which column are110 * considered to be dependent on previous ones and are discarded.111 * @param generator Underlying generator for uncorrelated normalized112 * components.113 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException114 * if the covariance matrix is not strictly positive definite.115 */116 public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,117 NormalizedRandomGenerator generator) {118 int order = covariance.getRowDimension();119 mean = new double[order];120 for (int i = 0; i < order; ++i) {121 mean[i] = 0;122 }123124 final RectangularCholeskyDecomposition decomposition =125 new RectangularCholeskyDecomposition(covariance, small);126 root = decomposition.getRootMatrix();127128 this.generator = generator;129 normalized = new double[decomposition.getRank()];130131 }132133 /** Get the underlying normalized components generator.134 * @return underlying uncorrelated components generator135 */136 public NormalizedRandomGenerator getGenerator() {137 return generator;138 }139140 /** Get the rank of the covariance matrix.141 * The rank is the number of independent rows in the covariance142 * matrix, it is also the number of columns of the root matrix.143 * @return rank of the square matrix.144 * @see #getRootMatrix()145 */146 public int getRank() {147 return normalized.length;148 }149150 /** Get the root of the covariance matrix.151 * The root is the rectangular matrix <code>B</code> such that152 * the covariance matrix is equal to <code>B.B<sup>T</sup></code>153 * @return root of the square matrix154 * @see #getRank()155 */156 public RealMatrix getRootMatrix() {157 return root;158 }159160 /** Generate a correlated random vector.161 * @return a random vector as an array of double. The returned array162 * is created at each call, the caller can do what it wants with it.163 */164 public double[] nextVector() {165166 // generate uncorrelated vector167 for (int i = 0; i < normalized.length; ++i) {168 normalized[i] = generator.nextNormalizedDouble();169 }170171 // compute correlated vector172 double[] correlated = new double[mean.length];173 for (int i = 0; i < correlated.length; ++i) {174 correlated[i] = mean[i];175 for (int j = 0; j < root.getColumnDimension(); ++j) {176 correlated[i] += root.getEntry(i, j) * normalized[j];177 }178 }179180 return correlated;181182 }183184}