You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by lu...@apache.org on 2007/02/28 20:51:49 UTC
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/
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