Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-1601: Simplified and more robust API. 

Factory methods ensure correct use (removed dependency on "NormalizedRandomGenerator").
This commit is contained in:
Gilles Sadowski 2021-06-02 02:53:26 +02:00
parent 8afd815000
commit 456de1bf98
6 changed files with 346 additions and 421 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/RectangularCholeskyDecomposition.java
View File

@ -29,7 +29,7 @@ import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
* is that rows/columns may be permuted (hence the rectangular shape instead
* of the traditional triangular shape) and there is a threshold to ignore
* small diagonal elements. This is used for example to generate {@link
* org.apache.commons.math4.legacy.random.CorrelatedRandomVectorGenerator correlated
* org.apache.commons.math4.legacy.random.CorrelatedVectorFactory correlated
* random n-dimensions vectors} in a p-dimension subspace (p < n).
* In other words, it allows generating random vectors from a covariance
* matrix that is only positive semidefinite, and not positive definite.</p>

184
commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGenerator.java
View File

@ -1,184 +0,0 @@
/*
* 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.math4.legacy.random;
import java.util.function.Supplier;
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RectangularCholeskyDecomposition;
/**
* 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>
* <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 UniformRandomGenerator}.</p>
* <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.</p>
*
* @since 1.2
*/
public class CorrelatedRandomVectorGenerator implements Supplier<double[]> {
/** 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.math4.legacy.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.math4.legacy.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.
*/
@Override
public double[] get() {
// 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;
}
}

163
commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactory.java Normal file
View File

@ -0,0 +1,163 @@
/*
* 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.math4.legacy.random;
import java.util.function.Supplier;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RectangularCholeskyDecomposition;
/**
* 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>
*
* <p>The main use of 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
* common 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>, usually using Cholesky decomposition.
* Other distributions are possible as long as the underlying sampler provides
* normalized values (unit standard deviation).</p>
*
* <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.</p>
*/
public class CorrelatedVectorFactory {
/** Square root of three. */
private static final double SQRT3 = AccurateMath.sqrt(3);
/** Mean vector. */
private final double[] mean;
/** Root of the covariance matrix. */
private final RealMatrix root;
/** Size of uncorrelated vector. */
private final int lengthUncorrelated;
/** Size of correlated vector. */
private final int lengthCorrelated;
/**
* Correlated vector factory.
*
* @param mean Expected mean values of the components.
* @param covariance Covariance matrix.
* @param small Diagonal elements threshold under which columns are
* considered to be dependent on previous ones and are discarded.
* @throws org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException
* if the covariance matrix is not strictly positive definite.
* @throws DimensionMismatchException if the mean and covariance
* arrays dimensions do not match.
*/
public CorrelatedVectorFactory(double[] mean,
RealMatrix covariance,
double small) {
lengthCorrelated = covariance.getRowDimension();
if (mean.length != lengthCorrelated) {
throw new DimensionMismatchException(mean.length, lengthCorrelated);
}
this.mean = mean.clone();
final RectangularCholeskyDecomposition decomposition
= new RectangularCholeskyDecomposition(covariance, small);
root = decomposition.getRootMatrix();
lengthUncorrelated = decomposition.getRank();
}
/**
* Null mean correlated vector factory.
*
* @param covariance Covariance matrix.
* @param small Diagonal elements threshold under which columns are
* considered to be dependent on previous ones and are discarded.
* @throws org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException
* if the covariance matrix is not strictly positive definite.
*/
public CorrelatedVectorFactory(RealMatrix covariance,
double small) {
this(new double[covariance.getRowDimension()],
covariance,
small);
}
/**
* @param rng RNG.
* @return a generator of vectors with correlated components sampled
* from a uniform distribution.
*/
public Supplier<double[]> uniform(UniformRandomProvider rng) {
return with(new ContinuousUniformSampler(rng, -SQRT3, SQRT3));
}
/**
* @param rng RNG.
* @return a generator of vectors with correlated components sampled
* from a normal distribution.
*/
public Supplier<double[]> gaussian(UniformRandomProvider rng) {
return with(new ZigguratNormalizedGaussianSampler(rng));
}
/**
* @param sampler Generator of samples from a normalized distribution.
* @return a generator of vectors with correlated components.
*/
private Supplier<double[]> with(final ContinuousSampler sampler) {
return new Supplier<double[]>() {
@Override
public double[] get() {
// Uncorrelated vector.
final double[] uncorrelated = new double[lengthUncorrelated];
for (int i = 0; i < lengthUncorrelated; i++) {
uncorrelated[i] = sampler.sample();
}
// Correlated vector.
final double[] correlated = mean.clone();
for (int i = 0; i < correlated.length; i++) {
for (int j = 0; j < lengthUncorrelated; j++) {
correlated[i] += root.getEntry(i, j) * uncorrelated[j];
}
}
return correlated;
}
};
}
}

208
commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGeneratorTest.java
View File

@ -1,208 +0,0 @@
//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.math4.legacy.random;
import java.util.Arrays;
import org.apache.commons.math4.legacy.TestUtils;
import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
import org.apache.commons.math4.legacy.linear.MatrixUtils;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
import org.apache.commons.math4.legacy.stat.correlation.StorelessCovariance;
import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialCovariance;
import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialMean;
import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
import org.junit.Test;
import org.junit.Assert;
public class CorrelatedRandomVectorGeneratorTest {
private double[] mean;
private RealMatrix covariance;
private CorrelatedRandomVectorGenerator generator;
public CorrelatedRandomVectorGeneratorTest() {
mean = new double[] { 0.0, 1.0, -3.0, 2.3 };
RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
int counter = 0;
for (int i = 0; i < b.getRowDimension(); ++i) {
for (int j = 0; j < b.getColumnDimension(); ++j) {
b.setEntry(i, j, 1.0 + 0.1 * ++counter);
}
}
RealMatrix bbt = b.multiply(b.transpose());
covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
for (int i = 0; i < covariance.getRowDimension(); ++i) {
covariance.setEntry(i, i, bbt.getEntry(i, i));
for (int j = 0; j < covariance.getColumnDimension(); ++j) {
double s = bbt.getEntry(i, j);
covariance.setEntry(i, j, s);
covariance.setEntry(j, i, s);
}
}
generator = new CorrelatedRandomVectorGenerator(mean,
covariance,
1e-12 * covariance.getNorm(),
gaussianRandom(RandomSource.create(RandomSource.WELL_1024_A)));
}
@Test
public void testRank() {
Assert.assertEquals(2, generator.getRank());
}
@Test
public void testMath226() {
double[] mean = { 1, 1, 10, 1 };
double[][] cov = {
{ 1, 3, 2, 6 },
{ 3, 13, 16, 2 },
{ 2, 16, 38, -1 },
{ 6, 2, -1, 197 }
};
RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
CorrelatedRandomVectorGenerator sg =
new CorrelatedRandomVectorGenerator(mean, covRM, 0.00001,
gaussianRandom(RandomSource.create(RandomSource.WELL_1024_A)));
double[] min = new double[mean.length];
Arrays.fill(min, Double.POSITIVE_INFINITY);
double[] max = new double[mean.length];
Arrays.fill(max, Double.NEGATIVE_INFINITY);
for (int i = 0; i < 10; i++) {
double[] generated = sg.get();
for (int j = 0; j < generated.length; ++j) {
min[j] = AccurateMath.min(min[j], generated[j]);
max[j] = AccurateMath.max(max[j], generated[j]);
}
}
for (int j = 0; j < min.length; ++j) {
Assert.assertTrue(max[j] - min[j] > 2.0);
}
}
@Test
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) {
Assert.assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
}
}
}
@Test
public void testMeanAndCovariance() {
VectorialMean meanStat = new VectorialMean(mean.length);
VectorialCovariance covStat = new VectorialCovariance(mean.length, true);
for (int i = 0; i < 5000; ++i) {
double[] v = generator.get();
meanStat.increment(v);
covStat.increment(v);
}
double[] estimatedMean = meanStat.getResult();
RealMatrix estimatedCovariance = covStat.getResult();
for (int i = 0; i < estimatedMean.length; ++i) {
Assert.assertEquals(mean[i], estimatedMean[i], 0.07);
for (int j = 0; j <= i; ++j) {
Assert.assertEquals(covariance.getEntry(i, j),
estimatedCovariance.getEntry(i, j),
0.1 * (1.0 + AccurateMath.abs(mean[i])) * (1.0 + AccurateMath.abs(mean[j])));
}
}
}
@Test
public void testSampleWithZeroCovariance() {
final double[][] covMatrix1 = new double[][]{
{0.013445532, 0.010394690, 0.009881156, 0.010499559},
{0.010394690, 0.023006616, 0.008196856, 0.010732709},
{0.009881156, 0.008196856, 0.019023866, 0.009210099},
{0.010499559, 0.010732709, 0.009210099, 0.019107243}
};
final double[][] covMatrix2 = new double[][]{
{0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 0.013445532, 0.010394690, 0.009881156, 0.010499559},
{0.0, 0.010394690, 0.023006616, 0.008196856, 0.010732709},
{0.0, 0.009881156, 0.008196856, 0.019023866, 0.009210099},
{0.0, 0.010499559, 0.010732709, 0.009210099, 0.019107243}
};
final double[][] covMatrix3 = new double[][]{
{0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
{0.010394690, 0.023006616, 0.0, 0.008196856, 0.010732709},
{0.0, 0.0, 0.0, 0.0, 0.0},
{0.009881156, 0.008196856, 0.0, 0.019023866, 0.009210099},
{0.010499559, 0.010732709, 0.0, 0.009210099, 0.019107243}
};
testSampler(covMatrix1, 10000, 0.001);
testSampler(covMatrix2, 10000, 0.001);
testSampler(covMatrix3, 10000, 0.001);
}
private CorrelatedRandomVectorGenerator createSampler(double[][] cov) {
RealMatrix matrix = new Array2DRowRealMatrix(cov);
double small = 1e-12 * matrix.getNorm();
return new CorrelatedRandomVectorGenerator(new double[cov.length],
matrix,
small,
gaussianRandom(RandomSource.create(RandomSource.WELL_1024_A)));
}
private void testSampler(final double[][] covMatrix, int samples, double epsilon) {
CorrelatedRandomVectorGenerator sampler = createSampler(covMatrix);
StorelessCovariance cov = new StorelessCovariance(covMatrix.length);
for (int i = 0; i < samples; ++i) {
cov.increment(sampler.get());
}
double[][] sampleCov = cov.getData();
for (int r = 0; r < covMatrix.length; ++r) {
TestUtils.assertEquals(covMatrix[r], sampleCov[r], epsilon);
}
}
/**
* @param rng RNG.
* @return a N(0,1) sampler.
*/
private NormalizedRandomGenerator gaussianRandom(final UniformRandomProvider rng) {
final ZigguratNormalizedGaussianSampler n = new ZigguratNormalizedGaussianSampler(rng);
return new NormalizedRandomGenerator() {
/** {@inheritDoc} */
@Override
public double nextNormalizedDouble() {
return n.sample();
}
};
}
}

172
commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactoryTest.java Normal file
View File

@ -0,0 +1,172 @@
/*
* 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.math4.legacy.random;
import java.util.Arrays;
import java.util.function.Supplier;
import org.junit.Test;
import org.junit.Assert;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
import org.apache.commons.math4.legacy.TestUtils;
import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
import org.apache.commons.math4.legacy.linear.MatrixUtils;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.stat.correlation.StorelessCovariance;
import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialCovariance;
import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialMean;
import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
public class CorrelatedVectorFactoryTest {
private double[] mean;
private RealMatrix covariance;
private Supplier<double[]> generator;
public CorrelatedVectorFactoryTest() {
mean = new double[] { 0.0, 1.0, -3.0, 2.3 };
final RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
int counter = 0;
for (int i = 0; i < b.getRowDimension(); ++i) {
for (int j = 0; j < b.getColumnDimension(); ++j) {
b.setEntry(i, j, 1.0 + 0.1 * ++counter);
}
}
final RealMatrix bbt = b.multiply(b.transpose());
covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
for (int i = 0; i < covariance.getRowDimension(); ++i) {
covariance.setEntry(i, i, bbt.getEntry(i, i));
for (int j = 0; j < covariance.getColumnDimension(); ++j) {
double s = bbt.getEntry(i, j);
covariance.setEntry(i, j, s);
covariance.setEntry(j, i, s);
}
}
generator = new CorrelatedVectorFactory(mean,
covariance,
1e-12 * covariance.getNorm())
.gaussian(RandomSource.create(RandomSource.KISS));
}
@Test
public void testMath226() {
final double[] mean = { 1, 1, 10, 1 };
final double[][] cov = {
{ 1, 3, 2, 6 },
{ 3, 13, 16, 2 },
{ 2, 16, 38, -1 },
{ 6, 2, -1, 197 }
};
final RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
final Supplier<double[]> sg = new CorrelatedVectorFactory(mean, covRM, 1e-5)
.gaussian(RandomSource.create(RandomSource.WELL_1024_A));
final double[] min = new double[mean.length];
Arrays.fill(min, Double.POSITIVE_INFINITY);
final double[] max = new double[mean.length];
Arrays.fill(max, Double.NEGATIVE_INFINITY);
for (int i = 0; i < 10; i++) {
double[] generated = sg.get();
for (int j = 0; j < generated.length; ++j) {
min[j] = AccurateMath.min(min[j], generated[j]);
max[j] = AccurateMath.max(max[j], generated[j]);
}
}
for (int j = 0; j < min.length; ++j) {
Assert.assertTrue(max[j] - min[j] > 2.0);
}
}
@Test
public void testMeanAndCovariance() {
final VectorialMean meanStat = new VectorialMean(mean.length);
final VectorialCovariance covStat = new VectorialCovariance(mean.length, true);
for (int i = 0; i < 5000; ++i) {
final double[] v = generator.get();
meanStat.increment(v);
covStat.increment(v);
}
final double[] estimatedMean = meanStat.getResult();
final RealMatrix estimatedCovariance = covStat.getResult();
for (int i = 0; i < estimatedMean.length; ++i) {
Assert.assertEquals(mean[i], estimatedMean[i], 0.07);
for (int j = 0; j <= i; ++j) {
Assert.assertEquals(covariance.getEntry(i, j),
estimatedCovariance.getEntry(i, j),
1e-1 * (1 + AccurateMath.abs(mean[i])) * (1 + AccurateMath.abs(mean[j])));
}
}
}
@Test
public void testSampleWithZeroCovariance() {
final double[][] covMatrix1 = new double[][]{
{0.013445532, 0.010394690, 0.009881156, 0.010499559},
{0.010394690, 0.023006616, 0.008196856, 0.010732709},
{0.009881156, 0.008196856, 0.019023866, 0.009210099},
{0.010499559, 0.010732709, 0.009210099, 0.019107243}
};
final double[][] covMatrix2 = new double[][]{
{0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 0.013445532, 0.010394690, 0.009881156, 0.010499559},
{0.0, 0.010394690, 0.023006616, 0.008196856, 0.010732709},
{0.0, 0.009881156, 0.008196856, 0.019023866, 0.009210099},
{0.0, 0.010499559, 0.010732709, 0.009210099, 0.019107243}
};
final double[][] covMatrix3 = new double[][]{
{0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
{0.010394690, 0.023006616, 0.0, 0.008196856, 0.010732709},
{0.0, 0.0, 0.0, 0.0, 0.0},
{0.009881156, 0.008196856, 0.0, 0.019023866, 0.009210099},
{0.010499559, 0.010732709, 0.0, 0.009210099, 0.019107243}
};
testSampler(covMatrix1, 10000, 1e-3);
testSampler(covMatrix2, 10000, 1e-3);
testSampler(covMatrix3, 10000, 1e-3);
}
private Supplier<double[]> createSampler(double[][] cov) {
final RealMatrix matrix = new Array2DRowRealMatrix(cov);
final double small = 1e-12 * matrix.getNorm();
return new CorrelatedVectorFactory(matrix, small)
.gaussian(RandomSource.create(RandomSource.WELL_1024_A));
}
private void testSampler(final double[][] covMatrix,
int samples,
double epsilon) {
final Supplier<double[]> sampler = createSampler(covMatrix);
final StorelessCovariance cov = new StorelessCovariance(covMatrix.length);
for (int i = 0; i < samples; ++i) {
cov.increment(sampler.get());
}
final double[][] sampleCov = cov.getData();
for (int r = 0; r < covMatrix.length; ++r) {
TestUtils.assertEquals(covMatrix[r], sampleCov[r], epsilon);
}
}
}

38
commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/regression/GLSMultipleLinearRegressionTest.java
View File

@ -16,22 +16,23 @@
*/
package org.apache.commons.math4.legacy.stat.regression;
import java.util.function.Supplier;
import org.junit.Assert;
import org.junit.Before;
import org.junit.Test;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.statistics.distribution.ContinuousDistribution;
import org.apache.commons.statistics.distribution.NormalDistribution;
import org.apache.commons.math4.legacy.TestUtils;
import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
import org.apache.commons.math4.legacy.exception.NullArgumentException;
import org.apache.commons.math4.legacy.linear.MatrixUtils;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RealVector;
import org.apache.commons.math4.legacy.random.NormalizedRandomGenerator;
import org.apache.commons.math4.legacy.random.CorrelatedRandomVectorGenerator;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler;
import org.apache.commons.statistics.distribution.ContinuousDistribution;
import org.apache.commons.statistics.distribution.NormalDistribution;
import org.apache.commons.math4.legacy.random.CorrelatedVectorFactory;
import org.apache.commons.math4.legacy.stat.correlation.Covariance;
import org.apache.commons.math4.legacy.stat.descriptive.DescriptiveStatistics;
@ -246,11 +247,8 @@ public class GLSMultipleLinearRegressionTest extends MultipleLinearRegressionAbs
RealMatrix cov = (new Covariance(errorSeeds)).getCovarianceMatrix();
// Create a CorrelatedRandomVectorGenerator to use to generate correlated errors
double[] errorMeans = new double[nObs]; // Counting on init to 0 here
CorrelatedRandomVectorGenerator gen
= new CorrelatedRandomVectorGenerator(errorMeans, cov,
1e-12 * cov.getNorm(),
gaussianRandom(rg));
final Supplier<double[]> gen
= new CorrelatedVectorFactory(cov, 1e-12 * cov.getNorm()).gaussian(rg);
// Now start generating models. Use Longley X matrix on LHS
// and Longley OLS beta vector as "true" beta. Generate
@ -299,20 +297,4 @@ public class GLSMultipleLinearRegressionTest extends MultipleLinearRegressionAbs
assert(olsBetaStats.getMean() > 1.5 * glsBetaStats.getMean());
assert(olsBetaStats.getStandardDeviation() > glsBetaStats.getStandardDeviation());
}
/**
* @param rng RNG.
* @return a N(0,1) sampler.
*/
private NormalizedRandomGenerator gaussianRandom(final UniformRandomProvider rng) {
final ZigguratNormalizedGaussianSampler n = new ZigguratNormalizedGaussianSampler(rng);
return new NormalizedRandomGenerator() {
/** {@inheritDoc} */
@Override
public double nextNormalizedDouble() {
return n.sample();
}
};
}
}