View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math4.legacy.random;
18  
19  import java.util.Arrays;
20  import java.util.function.Supplier;
21  
22  import org.junit.Test;
23  import org.junit.Assert;
24  
25  import org.apache.commons.rng.simple.RandomSource;
26  import org.apache.commons.math4.legacy.TestUtils;
27  import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
28  import org.apache.commons.math4.legacy.linear.MatrixUtils;
29  import org.apache.commons.math4.legacy.linear.RealMatrix;
30  import org.apache.commons.math4.legacy.stat.correlation.StorelessCovariance;
31  import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialCovariance;
32  import org.apache.commons.math4.legacy.stat.descriptive.moment.VectorialMean;
33  import org.apache.commons.math4.core.jdkmath.JdkMath;
34  
35  public class CorrelatedVectorFactoryTest {
36      private double[] mean;
37      private RealMatrix covariance;
38      private Supplier<double[]> generator;
39  
40      public CorrelatedVectorFactoryTest() {
41          mean = new double[] { 0.0, 1.0, -3.0, 2.3 };
42  
43          final RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
44          int counter = 0;
45          for (int i = 0; i < b.getRowDimension(); ++i) {
46              for (int j = 0; j < b.getColumnDimension(); ++j) {
47                  b.setEntry(i, j, 1.0 + 0.1 * ++counter);
48              }
49          }
50          final RealMatrix bbt = b.multiply(b.transpose());
51          covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
52          for (int i = 0; i < covariance.getRowDimension(); ++i) {
53              covariance.setEntry(i, i, bbt.getEntry(i, i));
54              for (int j = 0; j < covariance.getColumnDimension(); ++j) {
55                  double s = bbt.getEntry(i, j);
56                  covariance.setEntry(i, j, s);
57                  covariance.setEntry(j, i, s);
58              }
59          }
60  
61          generator = new CorrelatedVectorFactory(mean,
62                                                  covariance,
63                                                  1e-12 * covariance.getNorm())
64              .gaussian(RandomSource.KISS.create());
65      }
66  
67      @Test
68      public void testMath226() {
69          final double[] mean = { 1, 1, 10, 1 };
70          final double[][] cov = {
71              { 1, 3, 2, 6 },
72              { 3, 13, 16, 2 },
73              { 2, 16, 38, -1 },
74              { 6, 2, -1, 197 }
75          };
76          final RealMatrix covRM = MatrixUtils.createRealMatrix(cov);
77          final Supplier<double[]> sg = new CorrelatedVectorFactory(mean, covRM, 1e-5)
78              .gaussian(RandomSource.WELL_1024_A.create());
79  
80          final double[] min = new double[mean.length];
81          Arrays.fill(min, Double.POSITIVE_INFINITY);
82          final double[] max = new double[mean.length];
83          Arrays.fill(max, Double.NEGATIVE_INFINITY);
84          for (int i = 0; i < 10; i++) {
85              double[] generated = sg.get();
86              for (int j = 0; j < generated.length; ++j) {
87                  min[j] = JdkMath.min(min[j], generated[j]);
88                  max[j] = JdkMath.max(max[j], generated[j]);
89              }
90          }
91          for (int j = 0; j < min.length; ++j) {
92              Assert.assertTrue(max[j] - min[j] > 2.0);
93          }
94      }
95  
96      @Test
97      public void testMeanAndCovariance() {
98          final VectorialMean meanStat = new VectorialMean(mean.length);
99          final VectorialCovariance covStat = new VectorialCovariance(mean.length, true);
100         for (int i = 0; i < 5000; ++i) {
101             final double[] v = generator.get();
102             meanStat.increment(v);
103             covStat.increment(v);
104         }
105 
106         final double[] estimatedMean = meanStat.getResult();
107         final RealMatrix estimatedCovariance = covStat.getResult();
108         for (int i = 0; i < estimatedMean.length; ++i) {
109             Assert.assertEquals(mean[i], estimatedMean[i], 0.07);
110             for (int j = 0; j <= i; ++j) {
111                 Assert.assertEquals(covariance.getEntry(i, j),
112                                     estimatedCovariance.getEntry(i, j),
113                                     1e-1 * (1 + JdkMath.abs(mean[i])) * (1 + JdkMath.abs(mean[j])));
114             }
115         }
116     }
117 
118     @Test
119     public void testSampleWithZeroCovariance() {
120         final double[][] covMatrix1 = new double[][]{
121             {0.013445532, 0.010394690, 0.009881156, 0.010499559},
122             {0.010394690, 0.023006616, 0.008196856, 0.010732709},
123             {0.009881156, 0.008196856, 0.019023866, 0.009210099},
124             {0.010499559, 0.010732709, 0.009210099, 0.019107243}
125         };
126 
127         final double[][] covMatrix2 = new double[][]{
128             {0.0, 0.0, 0.0, 0.0, 0.0},
129             {0.0, 0.013445532, 0.010394690, 0.009881156, 0.010499559},
130             {0.0, 0.010394690, 0.023006616, 0.008196856, 0.010732709},
131             {0.0, 0.009881156, 0.008196856, 0.019023866, 0.009210099},
132             {0.0, 0.010499559, 0.010732709, 0.009210099, 0.019107243}
133         };
134 
135         final double[][] covMatrix3 = new double[][]{
136             {0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
137             {0.010394690, 0.023006616, 0.0, 0.008196856, 0.010732709},
138             {0.0, 0.0, 0.0, 0.0, 0.0},
139             {0.009881156, 0.008196856, 0.0, 0.019023866, 0.009210099},
140             {0.010499559, 0.010732709, 0.0, 0.009210099, 0.019107243}
141         };
142 
143         testSampler(covMatrix1, 10000, 1e-3);
144         testSampler(covMatrix2, 10000, 1e-3);
145         testSampler(covMatrix3, 10000, 1e-3);
146     }
147 
148     private Supplier<double[]> createSampler(double[][] cov) {
149         final RealMatrix matrix = new Array2DRowRealMatrix(cov);
150         final double small = 1e-12 * matrix.getNorm();
151         return new CorrelatedVectorFactory(matrix, small)
152             .gaussian(RandomSource.WELL_1024_A.create());
153     }
154 
155     private void testSampler(final double[][] covMatrix,
156                              int samples,
157                              double epsilon) {
158         final Supplier<double[]> sampler = createSampler(covMatrix);
159 
160         final StorelessCovariance cov = new StorelessCovariance(covMatrix.length);
161         for (int i = 0; i < samples; ++i) {
162             cov.increment(sampler.get());
163         }
164 
165         final double[][] sampleCov = cov.getData();
166         for (int r = 0; r < covMatrix.length; ++r) {
167             TestUtils.assertEquals(covMatrix[r], sampleCov[r], epsilon);
168         }
169     }
170 }