← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
 
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
 
 
18
package org.apache.commons.math.random;
 
19
 
 
20
import org.apache.commons.math.DimensionMismatchException;
 
21
import org.apache.commons.math.linear.MatrixUtils;
 
22
import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
 
23
import org.apache.commons.math.linear.RealMatrix;
 
24
 
 
25
/** 
 
26
 * A {@link RandomVectorGenerator} that generates vectors with with 
 
27
 * correlated components.
 
28
 * <p>Random vectors with correlated components are built by combining
 
29
 * the uncorrelated components of another random vector in such a way that
 
30
 * the resulting correlations are the ones specified by a positive
 
31
 * definite covariance matrix.</p>
 
32
 * <p>The main use for correlated random vector generation is for Monte-Carlo
 
33
 * simulation of physical problems with several variables, for example to
 
34
 * generate error vectors to be added to a nominal vector. A particularly
 
35
 * interesting case is when the generated vector should be drawn from a <a
 
36
 * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
 
37
 * Multivariate Normal Distribution</a>. The approach using a Cholesky
 
38
 * decomposition is quite usual in this case. However, it cas be extended
 
39
 * to other cases as long as the underlying random generator provides
 
40
 * {@link NormalizedRandomGenerator normalized values} like {@link
 
41
 * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
 
42
 * <p>Sometimes, the covariance matrix for a given simulation is not
 
43
 * strictly positive definite. This means that the correlations are
 
44
 * not all independent from each other. In this case, however, the non
 
45
 * strictly positive elements found during the Cholesky decomposition
 
46
 * of the covariance matrix should not be negative either, they
 
47
 * should be null. Another non-conventional extension handling this case
 
48
 * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
 
49
 * where <code>C</code> is the covariance matrix and <code>U</code>
 
50
 * is an uppertriangular matrix, we compute <code>C = B.B<sup>T</sup></code>
 
51
 * where <code>B</code> is a rectangular matrix having
 
52
 * more rows than columns. The number of columns of <code>B</code> is
 
53
 * the rank of the covariance matrix, and it is the dimension of the
 
54
 * uncorrelated random vector that is needed to compute the component
 
55
 * of the correlated vector. This class handles this situation
 
56
 * automatically.</p>
 
57
 *
 
58
 * @version $Revision: 781122 $ $Date: 2009-06-02 14:53:23 -0400 (Tue, 02 Jun 2009) $
 
59
 * @since 1.2
 
60
 */
 
61
 
 
62
public class CorrelatedRandomVectorGenerator
 
63
    implements RandomVectorGenerator {
 
64
 
 
65
    /** Simple constructor.
 
66
     * <p>Build a correlated random vector generator from its mean
 
67
     * vector and covariance matrix.</p>
 
68
     * @param mean expected mean values for all components
 
69
     * @param covariance covariance matrix
 
70
     * @param small diagonal elements threshold under which  column are
 
71
     * considered to be dependent on previous ones and are discarded
 
72
     * @param generator underlying generator for uncorrelated normalized
 
73
     * components
 
74
     * @exception IllegalArgumentException if there is a dimension
 
75
     * mismatch between the mean vector and the covariance matrix
 
76
     * @exception NotPositiveDefiniteMatrixException if the
 
77
     * covariance matrix is not strictly positive definite
 
78
     * @exception DimensionMismatchException if the mean and covariance
 
79
     * arrays dimensions don't match
 
80
     */
 
81
    public CorrelatedRandomVectorGenerator(double[] mean,
 
82
                                           RealMatrix covariance, double small,
 
83
                                           NormalizedRandomGenerator generator)
 
84
    throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
 
85
 
 
86
        int order = covariance.getRowDimension();
 
87
        if (mean.length != order) {
 
88
            throw new DimensionMismatchException(mean.length, order);
 
89
        }
 
90
        this.mean = mean.clone();
 
91
 
 
92
        decompose(covariance, small);
 
93
 
 
94
        this.generator = generator;
 
95
        normalized = new double[rank];
 
96
 
 
97
    }
 
98
 
 
99
    /** Simple constructor.
 
100
     * <p>Build a null mean random correlated vector generator from its
 
101
     * covariance matrix.</p>
 
102
     * @param covariance covariance matrix
 
103
     * @param small diagonal elements threshold under which  column are
 
104
     * considered to be dependent on previous ones and are discarded
 
105
     * @param generator underlying generator for uncorrelated normalized
 
106
     * components
 
107
     * @exception NotPositiveDefiniteMatrixException if the
 
108
     * covariance matrix is not strictly positive definite
 
109
     */
 
110
    public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
 
111
                                           NormalizedRandomGenerator generator)
 
112
    throws NotPositiveDefiniteMatrixException {
 
113
 
 
114
        int order = covariance.getRowDimension();
 
115
        mean = new double[order];
 
116
        for (int i = 0; i < order; ++i) {
 
117
            mean[i] = 0;
 
118
        }
 
119
 
 
120
        decompose(covariance, small);
 
121
 
 
122
        this.generator = generator;
 
123
        normalized = new double[rank];
 
124
 
 
125
    }
 
126
 
 
127
    /** Get the underlying normalized components generator.
 
128
     * @return underlying uncorrelated components generator
 
129
     */
 
130
    public NormalizedRandomGenerator getGenerator() {
 
131
        return generator;
 
132
    }
 
133
 
 
134
    /** Get the root of the covariance matrix.
 
135
     * The root is the rectangular matrix <code>B</code> such that
 
136
     * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
 
137
     * @return root of the square matrix
 
138
     * @see #getRank()
 
139
     */
 
140
    public RealMatrix getRootMatrix() {
 
141
        return root;
 
142
    }
 
143
 
 
144
    /** Get the rank of the covariance matrix.
 
145
     * The rank is the number of independent rows in the covariance
 
146
     * matrix, it is also the number of columns of the rectangular
 
147
     * matrix of the decomposition.
 
148
     * @return rank of the square matrix.
 
149
     * @see #getRootMatrix()
 
150
     */
 
151
    public int getRank() {
 
152
        return rank;
 
153
    }
 
154
 
 
155
    /** Decompose the original square matrix.
 
156
     * <p>The decomposition is based on a Choleski decomposition
 
157
     * where additional transforms are performed:
 
158
     * <ul>
 
159
     *   <li>the rows of the decomposed matrix are permuted</li>
 
160
     *   <li>columns with the too small diagonal element are discarded</li>
 
161
     *   <li>the matrix is permuted</li>
 
162
     * </ul>
 
163
     * This means that rather than computing M = U<sup>T</sup>.U where U
 
164
     * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
 
165
     * where B is a rectangular matrix.
 
166
     * @param covariance covariance matrix
 
167
     * @param small diagonal elements threshold under which  column are
 
168
     * considered to be dependent on previous ones and are discarded
 
169
     * @exception NotPositiveDefiniteMatrixException if the
 
170
     * covariance matrix is not strictly positive definite
 
171
     */
 
172
    private void decompose(RealMatrix covariance, double small)
 
173
    throws NotPositiveDefiniteMatrixException {
 
174
 
 
175
        int order = covariance.getRowDimension();
 
176
        double[][] c = covariance.getData();
 
177
        double[][] b = new double[order][order];
 
178
 
 
179
        int[] swap  = new int[order];
 
180
        int[] index = new int[order];
 
181
        for (int i = 0; i < order; ++i) {
 
182
            index[i] = i;
 
183
        }
 
184
 
 
185
        rank = 0;
 
186
        for (boolean loop = true; loop;) {
 
187
 
 
188
            // find maximal diagonal element
 
189
            swap[rank] = rank;
 
190
            for (int i = rank + 1; i < order; ++i) {
 
191
                int ii  = index[i];
 
192
                int isi = index[swap[i]];
 
193
                if (c[ii][ii] > c[isi][isi]) {
 
194
                    swap[rank] = i;
 
195
                }
 
196
            }
 
197
 
 
198
 
 
199
            // swap elements
 
200
            if (swap[rank] != rank) {
 
201
                int tmp = index[rank];
 
202
                index[rank] = index[swap[rank]];
 
203
                index[swap[rank]] = tmp;
 
204
            }
 
205
 
 
206
            // check diagonal element
 
207
            int ir = index[rank];
 
208
            if (c[ir][ir] < small) {
 
209
 
 
210
                if (rank == 0) {
 
211
                    throw new NotPositiveDefiniteMatrixException();
 
212
                }
 
213
 
 
214
                // check remaining diagonal elements
 
215
                for (int i = rank; i < order; ++i) {
 
216
                    if (c[index[i]][index[i]] < -small) {
 
217
                        // there is at least one sufficiently negative diagonal element,
 
218
                        // the covariance matrix is wrong
 
219
                        throw new NotPositiveDefiniteMatrixException();
 
220
                    }
 
221
                }
 
222
 
 
223
                // all remaining diagonal elements are close to zero,
 
224
                // we consider we have found the rank of the covariance matrix
 
225
                ++rank;
 
226
                loop = false;
 
227
 
 
228
            } else {
 
229
 
 
230
                // transform the matrix
 
231
                double sqrt = Math.sqrt(c[ir][ir]);
 
232
                b[rank][rank] = sqrt;
 
233
                double inverse = 1 / sqrt;
 
234
                for (int i = rank + 1; i < order; ++i) {
 
235
                    int ii = index[i];
 
236
                    double e = inverse * c[ii][ir];
 
237
                    b[i][rank] = e;
 
238
                    c[ii][ii] -= e * e;
 
239
                    for (int j = rank + 1; j < i; ++j) {
 
240
                        int ij = index[j];
 
241
                        double f = c[ii][ij] - e * b[j][rank];
 
242
                        c[ii][ij] = f;
 
243
                        c[ij][ii] = f;
 
244
                    }
 
245
                }
 
246
 
 
247
                // prepare next iteration
 
248
                loop = ++rank < order;
 
249
 
 
250
            }
 
251
 
 
252
        }
 
253
 
 
254
        // build the root matrix
 
255
        root = MatrixUtils.createRealMatrix(order, rank);
 
256
        for (int i = 0; i < order; ++i) {
 
257
            for (int j = 0; j < rank; ++j) {
 
258
                root.setEntry(index[i], j, b[i][j]);
 
259
            }
 
260
        }
 
261
 
 
262
    }
 
263
 
 
264
    /** Generate a correlated random vector.
 
265
     * @return a random vector as an array of double. The returned array
 
266
     * is created at each call, the caller can do what it wants with it.
 
267
     */
 
268
    public double[] nextVector() {
 
269
 
 
270
        // generate uncorrelated vector
 
271
        for (int i = 0; i < rank; ++i) {
 
272
            normalized[i] = generator.nextNormalizedDouble();
 
273
        }
 
274
 
 
275
        // compute correlated vector
 
276
        double[] correlated = new double[mean.length];
 
277
        for (int i = 0; i < correlated.length; ++i) {
 
278
            correlated[i] = mean[i];
 
279
            for (int j = 0; j < rank; ++j) {
 
280
                correlated[i] += root.getEntry(i, j) * normalized[j];
 
281
            }
 
282
        }
 
283
 
 
284
        return correlated;
 
285
 
 
286
    }
 
287
 
 
288
    /** Mean vector. */
 
289
    private double[] mean;
 
290
 
 
291
    /** Permutated Cholesky root of the covariance matrix. */
 
292
    private RealMatrix root;
 
293
 
 
294
    /** Rank of the covariance matrix. */
 
295
    private int rank;
 
296
 
 
297
    /** Underlying generator. */
 
298
    private NormalizedRandomGenerator generator;
 
299
 
 
300
    /** Storage for the normalized vector. */
 
301
    private double[] normalized;
 
302
 
 
303
}