← Back to branch summary

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

« back to all changes in this revision

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