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
18package org.apache.commons.math.random;
19
20import org.apache.commons.math.DimensionMismatchException;
21import org.apache.commons.math.linear.MatrixUtils;
22import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
23import org.apache.commons.math.linear.RealMatrix;
24import 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
63public 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}
305