commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r512933 - in /jakarta/commons/proper/math/trunk/src: java/org/apache/commons/math/optimization/ java/org/apache/commons/math/random/ test/org/apache/commons/math/random/
Date Wed, 28 Feb 2007 19:51:49 GMT
Author: luc
Date: Wed Feb 28 11:51:48 2007
New Revision: 512933

URL: http://svn.apache.org/viewvc?view=rev&rev=512933
Log:
implemented correlated random vectors generation

Added:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
  (with props)
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/random/CorrelatedRandomVectorGeneratorTest.java
      - copied, changed from r511516, jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/random/CorrelatedRandomVectorGeneratorTest.java
Modified:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/DirectSearchOptimizer.java

Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/DirectSearchOptimizer.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/DirectSearchOptimizer.java?view=diff&rev=512933&r1=512932&r2=512933
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/DirectSearchOptimizer.java
(original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/optimization/DirectSearchOptimizer.java
Wed Feb 28 11:51:48 2007
@@ -22,6 +22,7 @@
 
 import org.apache.commons.math.ConvergenceException;
 import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.linear.RealMatrix;
 import org.apache.commons.math.random.CorrelatedRandomVectorGenerator;
 import org.apache.commons.math.random.JDKRandomGenerator;
 import org.apache.commons.math.random.NotPositiveDefiniteMatrixException;
@@ -251,12 +252,15 @@
                 meanStat.increment(vertices[i]);
                 covStat.increment(vertices[i]);
             }
+            double[] mean = meanStat.getResult();
+            RealMatrix covariance = covStat.getResult();
+            
 
             RandomGenerator rg = new JDKRandomGenerator();
             rg.setSeed(seed);
             RandomVectorGenerator rvg =
-                new CorrelatedRandomVectorGenerator(meanStat.getResult(),
-                                                    covStat.getResult(),
+                new CorrelatedRandomVectorGenerator(mean,
+                                                    covariance, 1.0e-12 * covariance.getNorm(),
                                                     new UniformRandomGenerator(rg));
             setMultiStart(starts, rvg);
 

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java?view=auto&rev=512933
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
(added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
Wed Feb 28 11:51:48 2007
@@ -0,0 +1,289 @@
+//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.math.random;
+
+import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.linear.RealMatrixImpl;
+
+/** This class allows to generate random vectors with correlated components.
+
+ * <p>Random vectors with correlated components are built by combining
+ * the uncorrelated components of another random vector in such a way
+ * the resulting correlations are the ones specified by a positive
+ * definite covariance matrix.</p>
+
+ * <p>Sometimes, the covariance matrix for a given simulation is not
+ * strictly positive definite. This means that the correlations are
+ * not all independant 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. This implies that 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 uppertriangular 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 does handle this situation
+ * automatically.</p>
+
+ * @version $Revision:$ $Date$
+
+ */
+
+public class CorrelatedRandomVectorGenerator
+implements RandomVectorGenerator {
+
+    /** Simple constructor.
+     * <p>Build a correlated random vector generator from its mean
+     * vector and covariance matrix.</p>
+     * @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
+     * @exception IllegalArgumentException if there is a dimension
+     * mismatch between the mean vector and the covariance matrix
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     * @exception DimensionMismatchException if the mean and covariance
+     * arrays dimensions don't match
+     */
+    public CorrelatedRandomVectorGenerator(double[] mean,
+                                           RealMatrix covariance, double small,
+                                           NormalizedRandomGenerator generator)
+    throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
+
+        int order = covariance.getRowDimension();
+        if (mean.length != order) {
+            throw new DimensionMismatchException(mean.length, order);
+        }
+        this.mean = (double[]) mean.clone();
+
+        decompose(covariance, small);
+
+        this.generator = generator;
+        normalized = new double[rank];
+
+    }
+
+    /** Simple constructor.
+     * <p>Build a null mean random correlated vector generator from its
+     * covariance matrix.</p>
+     * @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
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     */
+    public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
+                                           NormalizedRandomGenerator generator)
+    throws NotPositiveDefiniteMatrixException {
+
+        int order = covariance.getRowDimension();
+        mean = new double[order];
+        for (int i = 0; i < order; ++i) {
+            mean[i] = 0;
+        }
+
+        decompose(covariance, small);
+
+        this.generator = generator;
+        normalized = new double[rank];
+
+    }
+
+    /** Get the underlying normalized components generator.
+     * @return underlying uncorrelated components generator
+     */
+    public NormalizedRandomGenerator getGenerator() {
+        return generator;
+    }
+
+    /** 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;
+    }
+
+    /** Get the rank of the covariance matrix.
+     * The rank is the number of independant rows in the covariance
+     * matrix, it is also the number of columns of the rectangular
+     * matrix of the decomposition.
+     * @return rank of the square matrix.
+     * @see #getRootMatrix()
+     */
+    public int getRank() {
+        return rank;
+    }
+
+    /** Decompose the original square matrix.
+     * <p>The decomposition is based on a Choleski decomposition
+     * where additional transforms are performed:
+     * <ul>
+     *   <li>the rows of the decomposed matrix are permuted</li>
+     *   <li>columns with the too small diagonal element are discarded</li>
+     *   <li>the matrix is permuted</li>
+     * </ul>
+     * This means that rather than computing M = U<sup>T</sup>.U where U
+     * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
+     * where B is a rectangular matrix.
+     * @param covariance covariance matrix
+     * @param small diagonal elements threshold under which  column are
+     * considered to be dependent on previous ones and are discarded
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     */
+    private void decompose(RealMatrix covariance, double small)
+    throws NotPositiveDefiniteMatrixException {
+
+        int order = covariance.getRowDimension();
+        double[][] c = covariance.getData();
+        double[][] b = new double[order][order];
+
+        int[] swap  = new int[order];
+        int[] index = new int[order];
+        for (int i = 0; i < order; ++i) {
+            index[i] = i;
+        }
+
+        rank = 0;
+        for (boolean loop = true; loop;) {
+
+            // find maximal diagonal element
+            swap[rank] = rank;
+            for (int i = rank + 1; i < order; ++i) {
+                int ii  = index[i];
+                int isi = index[swap[i]];
+                if (c[ii][ii] > c[isi][isi]) {
+                    swap[rank] = i;
+                }
+            }
+
+
+            // swap elements
+            if (swap[rank] != rank) {
+                int tmp = index[rank];
+                index[rank] = index[swap[rank]];
+                index[swap[rank]] = tmp;
+            }
+
+            // check diagonal element
+            int ir = index[rank];
+            if (c[ir][ir] < small) {
+
+                if (rank == 0) {
+                    throw new NotPositiveDefiniteMatrixException();
+                }
+
+                // check remaining diagonal elements
+                for (int i = rank; i < order; ++i) {
+                    if (c[index[i]][index[i]] < -small) {
+                        // there is at least one sufficiently negative diagonal element,
+                        // the covariance matrix is wrong
+                        throw new NotPositiveDefiniteMatrixException();
+                    }
+                }
+
+                // all remaining diagonal elements are close to zero,
+                // we consider we have found the rank of the covariance matrix
+                ++rank;
+                loop = false;
+
+            } else {
+
+                // transform the matrix
+                double sqrt = Math.sqrt(c[ir][ir]);
+                b[rank][rank] = sqrt;
+                double inverse = 1 / sqrt;
+                for (int i = rank + 1; i < order; ++i) {
+                    int ii = index[i];
+                    double e = inverse * c[ii][ir];
+                    b[i][rank] = e;
+                    c[ii][ii] -= e * e;
+                    for (int j = rank + 1; j < i; ++j) {
+                        int ij = index[j];
+                        double f = c[ii][ij] - e * b[j][rank];
+                        c[ii][ij] = f;
+                        c[ij][ii] = f;
+                    }
+                }
+
+                // prepare next iteration
+                loop = ++rank < order;
+
+            }
+
+        }
+
+        // build the root matrix
+        root = new RealMatrixImpl(order, rank);
+        for (int i = 0; i < order; ++i) {
+            System.arraycopy(b[i], 0, root.getDataRef()[swap[i]], 0, rank);
+        }
+
+    }
+
+    /** 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 < rank; ++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 < rank; ++j) {
+                correlated[i] += root.getEntry(i, j) * normalized[j];
+            }
+        }
+
+        return correlated;
+
+    }
+
+    /** Mean vector. */
+    private double[] mean;
+
+    /** Permutated Cholesky root of the covariance matrix. */
+    private RealMatrixImpl root;
+
+    /** Rank of the covariance matrix. */
+    private int rank;
+
+    /** Underlying generator. */
+    private NormalizedRandomGenerator generator;
+
+    /** Storage for the normalized vector. */
+    private double[] normalized;
+
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Copied: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/random/CorrelatedRandomVectorGeneratorTest.java
(from r511516, jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/random/CorrelatedRandomVectorGeneratorTest.java)
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/random/CorrelatedRandomVectorGeneratorTest.java?view=diff&rev=512933&p1=jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/random/CorrelatedRandomVectorGeneratorTest.java&r1=511516&p2=jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/random/CorrelatedRandomVectorGeneratorTest.java&r2=512933
==============================================================================
--- jakarta/commons/proper/math/trunk/src/mantissa/tests-src/org/spaceroots/mantissa/random/CorrelatedRandomVectorGeneratorTest.java
(original)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/random/CorrelatedRandomVectorGeneratorTest.java
Wed Feb 28 11:51:48 2007
@@ -1,114 +1,128 @@
-// 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.spaceroots.mantissa.random;
-
-import org.spaceroots.mantissa.linalg.Matrix;
-import org.spaceroots.mantissa.linalg.GeneralMatrix;
-import org.spaceroots.mantissa.linalg.SymetricalMatrix;
+//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.math.random;
+
+import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.linear.RealMatrixImpl;
+import org.apache.commons.math.stat.descriptive.moment.VectorialCovariance;
+import org.apache.commons.math.stat.descriptive.moment.VectorialMean;
 
 import junit.framework.*;
 
 public class CorrelatedRandomVectorGeneratorTest
-  extends TestCase {
+extends TestCase {
 
-  public CorrelatedRandomVectorGeneratorTest(String name) {
-    super(name);
-    mean       = null;
-    covariance = null;
-    generator  = null;
-  }
-
-  public void testRank() {
-    assertEquals(3, generator.getRank());
-  }
-
-  public void testRootMatrix() {
-    Matrix b = generator.getRootMatrix();
-    Matrix bbt = b.mul(b.getTranspose());
-    for (int i = 0; i < covariance.getRows(); ++i) {
-      for (int j = 0; j < covariance.getColumns(); ++j) {
-        assertEquals(covariance.getElement(i, j),
-                     bbt.getElement(i, j),
-                     1.0e-12);
-      }
-    }
-  }
-
-  public void testMeanAndCovariance() {
-
-    VectorialSampleStatistics sample = new VectorialSampleStatistics();
-    for (int i = 0; i < 5000; ++i) {
-      sample.add(generator.nextVector());
-    }
-
-    double[] estimatedMean = sample.getMean();
-    SymetricalMatrix estimatedCovariance = sample.getCovarianceMatrix(null);
-    for (int i = 0; i < estimatedMean.length; ++i) {
-      assertEquals(mean[i], estimatedMean[i], 0.07);
-      for (int j = 0; j <= i; ++j) {
-        assertEquals(covariance.getElement(i, j),
-                     estimatedCovariance.getElement(i, j),
-                     0.1 * (1.0 + Math.abs(mean[i])) * (1.0 + Math.abs(mean[j])));
-      }
-    }
-
-  }
-
-  public void setUp() {
-    try {
-      mean = new double[] { 0.0, 1.0, -3.0, 2.3};
-
-      GeneralMatrix b = new GeneralMatrix(4, 3);
-      int counter = 0;
-      for (int i = 0; i < b.getRows(); ++i) {
-        for (int j = 0; j < b.getColumns(); ++j) {
-          b.setElement(i, j, 1.0 + 0.1 * ++counter);
+    public CorrelatedRandomVectorGeneratorTest(String name) {
+        super(name);
+        mean       = null;
+        covariance = null;
+        generator  = null;
+    }
+
+    public void testRank() {
+        assertEquals(3, generator.getRank());
+    }
+
+    public void testRootMatrix() {
+        RealMatrix b = generator.getRootMatrix();
+        RealMatrix bbt = b.multiply(b.transpose());
+        for (int i = 0; i < covariance.getRowDimension(); ++i) {
+            for (int j = 0; j < covariance.getColumnDimension(); ++j) {
+                assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
+            }
+        }
+    }
+
+    public void testMeanAndCovariance() throws DimensionMismatchException {
+
+        VectorialMean meanStat = new VectorialMean(mean.length);
+        VectorialCovariance covStat = new VectorialCovariance(mean.length);
+        for (int i = 0; i < 5000; ++i) {
+            double[] v = generator.nextVector();
+            meanStat.increment(v);
+            covStat.increment(v);
+        }
+
+        double[] estimatedMean = meanStat.getResult();
+        RealMatrix estimatedCovariance = covStat.getResult();
+        for (int i = 0; i < estimatedMean.length; ++i) {
+            assertEquals(mean[i], estimatedMean[i], 0.07);
+            for (int j = 0; j <= i; ++j) {
+                assertEquals(covariance.getEntry(i, j),
+                        estimatedCovariance.getEntry(i, j),
+                        0.1 * (1.0 + Math.abs(mean[i])) * (1.0 + Math.abs(mean[j])));
+            }
         }
-      }
-      Matrix bbt = b.mul(b.getTranspose());
-      covariance = new SymetricalMatrix(mean.length);
-      for (int i = 0; i < covariance.getRows(); ++i) {
-        covariance.setElement(i, i, bbt.getElement(i, i));
-        for (int j = 0; j < covariance.getColumns(); ++j) {
-          covariance.setElementAndSymetricalElement(i, j,
-                                                    bbt.getElement(i, j));
+
+    }
+
+    public void setUp() {
+        try {
+            mean = new double[] { 0.0, 1.0, -3.0, 2.3};
+
+            RealMatrixImpl b = new RealMatrixImpl(4, 3);
+            double[][] bData = b.getDataRef();
+            int counter = 0;
+            for (int i = 0; i < bData.length; ++i) {
+                double[] bi = bData[i];
+                for (int j = 0; j < b.getColumnDimension(); ++j) {
+                    bi[j] = 1.0 + 0.1 * ++counter;
+                }
+            }
+            RealMatrix bbt = b.multiply(b.transpose());
+            covariance = new RealMatrixImpl(mean.length, mean.length);
+            double[][] covData = covariance.getDataRef();
+            for (int i = 0; i < covariance.getRowDimension(); ++i) {
+                covData[i][i] = bbt.getEntry(i, i);
+                for (int j = 0; j < covariance.getColumnDimension(); ++j) {
+                    double s = bbt.getEntry(i, j);
+                    covData[i][j] = s;
+                    covData[j][i] = s;
+                }
+            }
+
+            RandomGenerator rg = new JDKRandomGenerator();
+            rg.setSeed(17399225432l);
+            GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
+            generator = new CorrelatedRandomVectorGenerator(mean,
+                                                            covariance,
+                                                            1.0e-12 * covariance.getNorm(),
+                                                            rawGenerator);
+        } catch (DimensionMismatchException e) {
+            fail(e.getMessage());
+        } catch (NotPositiveDefiniteMatrixException e) {
+            fail("not positive definite matrix");
         }
-      }
+    }
+
+    public void tearDown() {
+        mean       = null;
+        covariance = null;
+        generator  = null;
+    }
+
+    public static Test suite() {
+        return new TestSuite(CorrelatedRandomVectorGeneratorTest.class);
+    }
 
-      GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(17399225432l);
-      generator = new CorrelatedRandomVectorGenerator(mean, covariance, rawGenerator);
-    } catch (NotPositiveDefiniteMatrixException e) {
-      fail("not positive definite matrix");
-    }
-  }
-
-  public void tearDown() {
-    mean       = null;
-    covariance = null;
-    generator  = null;
-  }
-
-  public static Test suite() {
-    return new TestSuite(CorrelatedRandomVectorGeneratorTest.class);
-  }
-
-  private double[] mean;
-  private SymetricalMatrix covariance;
-  private CorrelatedRandomVectorGenerator generator;
+    private double[] mean;
+    private RealMatrixImpl covariance;
+    private CorrelatedRandomVectorGenerator generator;
 
 }



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message