1 files changed, 178 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java
new file mode 100644
index 0000000..6668356
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/random/CorrelatedRandomVectorGenerator.java
@@ -0,0 +1,178 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math3.random;
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.linear.RectangularCholeskyDecomposition;
+
+/**
+ * A {@link RandomVectorGenerator} that generates vectors with with correlated components.
+ *
+ * <p>Random vectors with correlated components are built by combining the uncorrelated components
+ * of another random vector in such a way that the resulting correlations are the ones specified by
+ * a positive definite covariance matrix.
+ *
+ * <p>The main use for correlated random vector generation is for Monte-Carlo simulation of physical
+ * problems with several variables, for example to generate error vectors to be added to a nominal
+ * vector. A particularly interesting case is when the generated vector should be drawn from a <a
+ * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">Multivariate Normal
+ * Distribution</a>. The approach using a Cholesky decomposition is quite usual in this case.
+ * However, it can be extended to other cases as long as the underlying random generator provides
+ * {@link NormalizedRandomGenerator normalized values} like {@link GaussianRandomGenerator} or
+ * {@link UniformRandomGenerator}.
+ *
+ * <p>Sometimes, the covariance matrix for a given simulation is not strictly positive definite.
+ * This means that the correlations are not all independent from each other. In this case, however,
+ * the non strictly positive elements found during the Cholesky decomposition of the covariance
+ * matrix should not be negative either, they should be null. Another non-conventional extension
+ * handling this case is used here. Rather than computing <code>C = U<sup>T</sup>.U</code> where
+ * <code>C</code> is the covariance matrix and <code>U</code> is an upper-triangular matrix, we
+ * compute <code>C = B.B<sup>T</sup></code> where <code>B</code> is a rectangular matrix having more
+ * rows than columns. The number of columns of <code>B</code> is the rank of the covariance matrix,
+ * and it is the dimension of the uncorrelated random vector that is needed to compute the component
+ * of the correlated vector. This class handles this situation automatically.
+ *
+ * @since 1.2
+ */
+public class CorrelatedRandomVectorGenerator implements RandomVectorGenerator {
+    /** Mean vector. */
+    private final double[] mean;
+
+    /** Underlying generator. */
+    private final NormalizedRandomGenerator generator;
+
+    /** Storage for the normalized vector. */
+    private final double[] normalized;
+
+    /** Root of the covariance matrix. */
+    private final RealMatrix root;
+
+    /**
+     * Builds a correlated random vector generator from its mean vector and covariance matrix.
+     *
+     * @param mean Expected mean values for all components.
+     * @param covariance Covariance matrix.
+     * @param small Diagonal elements threshold under which column are considered to be dependent on
+     *     previous ones and are discarded
+     * @param generator underlying generator for uncorrelated normalized components.
+     * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException if the covariance
+     *     matrix is not strictly positive definite.
+     * @throws DimensionMismatchException if the mean and covariance arrays dimensions do not match.
+     */
+    public CorrelatedRandomVectorGenerator(
+            double[] mean,
+            RealMatrix covariance,
+            double small,
+            NormalizedRandomGenerator generator) {
+        int order = covariance.getRowDimension();
+        if (mean.length != order) {
+            throw new DimensionMismatchException(mean.length, order);
+        }
+        this.mean = mean.clone();
+
+        final RectangularCholeskyDecomposition decomposition =
+                new RectangularCholeskyDecomposition(covariance, small);
+        root = decomposition.getRootMatrix();
+
+        this.generator = generator;
+        normalized = new double[decomposition.getRank()];
+    }
+
+    /**
+     * Builds a null mean random correlated vector generator from its covariance matrix.
+     *
+     * @param covariance Covariance matrix.
+     * @param small Diagonal elements threshold under which column are considered to be dependent on
+     *     previous ones and are discarded.
+     * @param generator Underlying generator for uncorrelated normalized components.
+     * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException if the covariance
+     *     matrix is not strictly positive definite.
+     */
+    public CorrelatedRandomVectorGenerator(
+            RealMatrix covariance, double small, NormalizedRandomGenerator generator) {
+        int order = covariance.getRowDimension();
+        mean = new double[order];
+        for (int i = 0; i < order; ++i) {
+            mean[i] = 0;
+        }
+
+        final RectangularCholeskyDecomposition decomposition =
+                new RectangularCholeskyDecomposition(covariance, small);
+        root = decomposition.getRootMatrix();
+
+        this.generator = generator;
+        normalized = new double[decomposition.getRank()];
+    }
+
+    /**
+     * Get the underlying normalized components generator.
+     *
+     * @return underlying uncorrelated components generator
+     */
+    public NormalizedRandomGenerator getGenerator() {
+        return generator;
+    }
+
+    /**
+     * Get the rank of the covariance matrix. The rank is the number of independent rows in the
+     * covariance matrix, it is also the number of columns of the root matrix.
+     *
+     * @return rank of the square matrix.
+     * @see #getRootMatrix()
+     */
+    public int getRank() {
+        return normalized.length;
+    }
+
+    /**
+     * Get the root of the covariance matrix. The root is the rectangular matrix <code>B</code> such
+     * that the covariance matrix is equal to <code>B.B<sup>T</sup></code>
+     *
+     * @return root of the square matrix
+     * @see #getRank()
+     */
+    public RealMatrix getRootMatrix() {
+        return root;
+    }
+
+    /**
+     * Generate a correlated random vector.
+     *
+     * @return a random vector as an array of double. The returned array is created at each call,
+     *     the caller can do what it wants with it.
+     */
+    public double[] nextVector() {
+
+        // generate uncorrelated vector
+        for (int i = 0; i < normalized.length; ++i) {
+            normalized[i] = generator.nextNormalizedDouble();
+        }
+
+        // compute correlated vector
+        double[] correlated = new double[mean.length];
+        for (int i = 0; i < correlated.length; ++i) {
+            correlated[i] = mean[i];
+            for (int j = 0; j < root.getColumnDimension(); ++j) {
+                correlated[i] += root.getEntry(i, j) * normalized[j];
+            }
+        }
+
+        return correlated;
+    }
+}