implemented correlated random vectors generation
git-svn-id: https://svn.apache.org/repos/asf/jakarta/commons/proper/math/trunk@512933 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
12dd062aec
commit
9ddf255a72
|
|
@ -22,6 +22,7 @@ import java.util.Comparator;
|
|||
|
||||
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 @@ public abstract class DirectSearchOptimizer {
|
|||
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);
|
||||
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +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.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 {
|
||||
|
||||
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])));
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
private double[] mean;
|
||||
private RealMatrixImpl covariance;
|
||||
private CorrelatedRandomVectorGenerator generator;
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue