diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/RectangularCholeskyDecomposition.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/RectangularCholeskyDecomposition.java index 6614cf15f..68882a51b 100644 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/RectangularCholeskyDecomposition.java +++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/RectangularCholeskyDecomposition.java @@ -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.

diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGenerator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGenerator.java deleted file mode 100644 index 789c05aad..000000000 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGenerator.java +++ /dev/null @@ -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. - * - *

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.

- *

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 - * Multivariate Normal Distribution. 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}.

- *

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 C = UT.U - * where C is the covariance matrix and U - * is an upper-triangular matrix, we compute C = B.BT - * where B is a rectangular matrix having - * more rows than columns. The number of columns of B 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.

- * - * @since 1.2 - */ - -public class CorrelatedRandomVectorGenerator implements Supplier { - /** 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 B such that - * the covariance matrix is equal to B.BT - * @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; - - } -} diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactory.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactory.java new file mode 100644 index 000000000..2fc3a600c --- /dev/null +++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactory.java @@ -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. + * + *

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.

+ * + *

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 + * + * Multivariate Normal Distribution, usually using Cholesky decomposition. + * Other distributions are possible as long as the underlying sampler provides + * normalized values (unit standard deviation).

+ * + *

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 C = UT U + * where C is the covariance matrix and U + * is an upper-triangular matrix, we compute C = B BT + * where B is a rectangular matrix having more rows than + * columns. The number of columns of B 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.

+ */ +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 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 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 with(final ContinuousSampler sampler) { + return new Supplier() { + @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; + } + }; + } +} diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGeneratorTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGeneratorTest.java deleted file mode 100644 index a6622c6f0..000000000 --- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedRandomVectorGeneratorTest.java +++ /dev/null @@ -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(); - } - }; - } -} diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactoryTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactoryTest.java new file mode 100644 index 000000000..89aaf5dfa --- /dev/null +++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/CorrelatedVectorFactoryTest.java @@ -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 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 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 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 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); + } + } +} diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/regression/GLSMultipleLinearRegressionTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/regression/GLSMultipleLinearRegressionTest.java index 0e7bf3ffa..12f080720 100644 --- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/regression/GLSMultipleLinearRegressionTest.java +++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/stat/regression/GLSMultipleLinearRegressionTest.java @@ -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 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(); - } - }; - } }