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 18 package org.apache.commons.math4.legacy.random; 19 20 import java.util.function.Supplier; 21 22 import org.apache.commons.rng.UniformRandomProvider; 23 import org.apache.commons.rng.sampling.distribution.ContinuousSampler; 24 import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler; 25 import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; 26 import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 27 import org.apache.commons.math4.core.jdkmath.JdkMath; 28 import org.apache.commons.math4.legacy.linear.RealMatrix; 29 import org.apache.commons.math4.legacy.linear.RectangularCholeskyDecomposition; 30 31 /** 32 * Generates vectors with with correlated components. 33 * 34 * <p>Random vectors with correlated components are built by combining 35 * the uncorrelated components of another random vector in such a way 36 * that the resulting correlations are the ones specified by a positive 37 * definite covariance matrix.</p> 38 * 39 * <p>The main use of correlated random vector generation is for Monte-Carlo 40 * simulation of physical problems with several variables (for example to 41 * generate error vectors to be added to a nominal vector). A particularly 42 * common case is when the generated vector should be drawn from a 43 * <a href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution"> 44 * Multivariate Normal Distribution</a>, usually using Cholesky decomposition. 45 * Other distributions are possible as long as the underlying sampler provides 46 * normalized values (unit standard deviation).</p> 47 * 48 * <p>Sometimes, the covariance matrix for a given simulation is not 49 * strictly positive definite. This means that the correlations are 50 * not all independent from each other. In this case, however, the non 51 * strictly positive elements found during the Cholesky decomposition 52 * of the covariance matrix should not be negative either, they 53 * should be null. Another non-conventional extension handling this case 54 * is used here. Rather than computing <code>C = U<sup>T</sup> U</code> 55 * where <code>C</code> is the covariance matrix and <code>U</code> 56 * is an upper-triangular matrix, we compute <code>C = B B<sup>T</sup></code> 57 * where <code>B</code> is a rectangular matrix having more rows than 58 * columns. The number of columns of <code>B</code> is the rank of the 59 * covariance matrix, and it is the dimension of the uncorrelated 60 * random vector that is needed to compute the component of the 61 * correlated vector. This class handles this situation automatically.</p> 62 */ 63 public class CorrelatedVectorFactory { 64 /** Square root of three. */ 65 private static final double SQRT3 = JdkMath.sqrt(3); 66 /** Mean vector. */ 67 private final double[] mean; 68 /** Root of the covariance matrix. */ 69 private final RealMatrix root; 70 /** Size of uncorrelated vector. */ 71 private final int lengthUncorrelated; 72 /** Size of correlated vector. */ 73 private final int lengthCorrelated; 74 75 /** 76 * Correlated vector factory. 77 * 78 * @param mean Expected mean values of the components. 79 * @param covariance Covariance matrix. 80 * @param small Diagonal elements threshold under which columns are 81 * considered to be dependent on previous ones and are discarded. 82 * @throws org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException 83 * if the covariance matrix is not strictly positive definite. 84 * @throws DimensionMismatchException if the mean and covariance 85 * arrays dimensions do not match. 86 */ 87 public CorrelatedVectorFactory(double[] mean, 88 RealMatrix covariance, 89 double small) { 90 lengthCorrelated = covariance.getRowDimension(); 91 if (mean.length != lengthCorrelated) { 92 throw new DimensionMismatchException(mean.length, lengthCorrelated); 93 } 94 this.mean = mean.clone(); 95 96 final RectangularCholeskyDecomposition decomposition 97 = new RectangularCholeskyDecomposition(covariance, small); 98 root = decomposition.getRootMatrix(); 99 100 lengthUncorrelated = decomposition.getRank(); 101 } 102 103 /** 104 * Null mean correlated vector factory. 105 * 106 * @param covariance Covariance matrix. 107 * @param small Diagonal elements threshold under which columns are 108 * considered to be dependent on previous ones and are discarded. 109 * @throws org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException 110 * if the covariance matrix is not strictly positive definite. 111 */ 112 public CorrelatedVectorFactory(RealMatrix covariance, 113 double small) { 114 this(new double[covariance.getRowDimension()], 115 covariance, 116 small); 117 } 118 119 /** 120 * @param rng RNG. 121 * @return a generator of vectors with correlated components sampled 122 * from a uniform distribution. 123 */ 124 public Supplier<double[]> uniform(UniformRandomProvider rng) { 125 return with(new ContinuousUniformSampler(rng, -SQRT3, SQRT3)); 126 } 127 128 /** 129 * @param rng RNG. 130 * @return a generator of vectors with correlated components sampled 131 * from a normal distribution. 132 */ 133 public Supplier<double[]> gaussian(UniformRandomProvider rng) { 134 return with(new ZigguratNormalizedGaussianSampler(rng)); 135 } 136 137 /** 138 * @param sampler Generator of samples from a normalized distribution. 139 * @return a generator of vectors with correlated components. 140 */ 141 private Supplier<double[]> with(final ContinuousSampler sampler) { 142 return new Supplier<double[]>() { 143 @Override 144 public double[] get() { 145 // Uncorrelated vector. 146 final double[] uncorrelated = new double[lengthUncorrelated]; 147 for (int i = 0; i < lengthUncorrelated; i++) { 148 uncorrelated[i] = sampler.sample(); 149 } 150 151 // Correlated vector. 152 final double[] correlated = mean.clone(); 153 for (int i = 0; i < correlated.length; i++) { 154 for (int j = 0; j < lengthUncorrelated; j++) { 155 correlated[i] += root.getEntry(i, j) * uncorrelated[j]; 156 } 157 } 158 159 return correlated; 160 } 161 }; 162 } 163 }