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: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]