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.math3.random;
019
020 import org.apache.commons.math3.exception.DimensionMismatchException;
021 import org.apache.commons.math3.linear.RealMatrix;
022 import org.apache.commons.math3.linear.RectangularCholeskyDecomposition;
023
024 /**
025 * A {@link RandomVectorGenerator} that generates vectors with with
026 * correlated components.
027 * <p>Random vectors with correlated components are built by combining
028 * the uncorrelated components of another random vector in such a way that
029 * the resulting correlations are the ones specified by a positive
030 * definite covariance matrix.</p>
031 * <p>The main use for correlated random vector generation is for Monte-Carlo
032 * simulation of physical problems with several variables, for example to
033 * generate error vectors to be added to a nominal vector. A particularly
034 * interesting case is when the generated vector should be drawn from a <a
035 * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
036 * Multivariate Normal Distribution</a>. The approach using a Cholesky
037 * decomposition is quite usual in this case. However, it can be extended
038 * to other cases as long as the underlying random generator provides
039 * {@link NormalizedRandomGenerator normalized values} like {@link
040 * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
041 * <p>Sometimes, the covariance matrix for a given simulation is not
042 * strictly positive definite. This means that the correlations are
043 * not all independent from each other. In this case, however, the non
044 * strictly positive elements found during the Cholesky decomposition
045 * of the covariance matrix should not be negative either, they
046 * should be null. Another non-conventional extension handling this case
047 * 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 having
051 * more rows than columns. The number of columns of <code>B</code> is
052 * the rank of the covariance matrix, and it is the dimension of the
053 * uncorrelated random vector that is needed to compute the component
054 * of the correlated vector. This class handles this situation
055 * automatically.</p>
056 *
057 * @version $Id: CorrelatedRandomVectorGenerator.java 1416643 2012-12-03 19:37:14Z tn $
058 * @since 1.2
059 */
060
061 public class CorrelatedRandomVectorGenerator
062 implements RandomVectorGenerator {
063 /** Mean vector. */
064 private final double[] mean;
065 /** Underlying generator. */
066 private final NormalizedRandomGenerator generator;
067 /** Storage for the normalized vector. */
068 private final double[] normalized;
069 /** Root of the covariance matrix. */
070 private final RealMatrix root;
071
072 /**
073 * Builds a correlated random vector generator from its mean
074 * vector and covariance matrix.
075 *
076 * @param mean Expected mean values for all components.
077 * @param covariance Covariance matrix.
078 * @param small Diagonal elements threshold under which column are
079 * considered to be dependent on previous ones and are discarded
080 * @param generator underlying generator for uncorrelated normalized
081 * components.
082 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException
083 * if the covariance matrix is not strictly positive definite.
084 * @throws DimensionMismatchException if the mean and covariance
085 * arrays dimensions do not match.
086 */
087 public CorrelatedRandomVectorGenerator(double[] mean,
088 RealMatrix covariance, double small,
089 NormalizedRandomGenerator generator) {
090 int order = covariance.getRowDimension();
091 if (mean.length != order) {
092 throw new DimensionMismatchException(mean.length, order);
093 }
094 this.mean = mean.clone();
095
096 final RectangularCholeskyDecomposition decomposition =
097 new RectangularCholeskyDecomposition(covariance, small);
098 root = decomposition.getRootMatrix();
099
100 this.generator = generator;
101 normalized = new double[decomposition.getRank()];
102
103 }
104
105 /**
106 * Builds a null mean random correlated vector generator from its
107 * covariance matrix.
108 *
109 * @param covariance Covariance matrix.
110 * @param small Diagonal elements threshold under which column are
111 * considered to be dependent on previous ones and are discarded.
112 * @param generator Underlying generator for uncorrelated normalized
113 * components.
114 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException
115 * if the covariance matrix is not strictly positive definite.
116 */
117 public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
118 NormalizedRandomGenerator generator) {
119 int order = covariance.getRowDimension();
120 mean = new double[order];
121 for (int i = 0; i < order; ++i) {
122 mean[i] = 0;
123 }
124
125 final RectangularCholeskyDecomposition decomposition =
126 new RectangularCholeskyDecomposition(covariance, small);
127 root = decomposition.getRootMatrix();
128
129 this.generator = generator;
130 normalized = new double[decomposition.getRank()];
131
132 }
133
134 /** Get the underlying normalized components generator.
135 * @return underlying uncorrelated components generator
136 */
137 public NormalizedRandomGenerator getGenerator() {
138 return generator;
139 }
140
141 /** Get the rank of the covariance matrix.
142 * The rank is the number of independent rows in the covariance
143 * matrix, it is also the number of columns of the root matrix.
144 * @return rank of the square matrix.
145 * @see #getRootMatrix()
146 */
147 public int getRank() {
148 return normalized.length;
149 }
150
151 /** Get the root of the covariance matrix.
152 * The root is the rectangular matrix <code>B</code> such that
153 * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
154 * @return root of the square matrix
155 * @see #getRank()
156 */
157 public RealMatrix getRootMatrix() {
158 return root;
159 }
160
161 /** Generate a correlated random vector.
162 * @return a random vector as an array of double. The returned array
163 * is created at each call, the caller can do what it wants with it.
164 */
165 public double[] nextVector() {
166
167 // generate uncorrelated vector
168 for (int i = 0; i < normalized.length; ++i) {
169 normalized[i] = generator.nextNormalizedDouble();
170 }
171
172 // compute correlated vector
173 double[] correlated = new double[mean.length];
174 for (int i = 0; i < correlated.length; ++i) {
175 correlated[i] = mean[i];
176 for (int j = 0; j < root.getColumnDimension(); ++j) {
177 correlated[i] += root.getEntry(i, j) * normalized[j];
178 }
179 }
180
181 return correlated;
182
183 }
184
185 }