001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math4.legacy.random; 019 020import java.util.function.Supplier; 021 022import org.apache.commons.rng.UniformRandomProvider; 023import org.apache.commons.rng.sampling.distribution.ContinuousSampler; 024import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler; 025import org.apache.commons.rng.sampling.distribution.ZigguratNormalizedGaussianSampler; 026import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 027import org.apache.commons.math4.core.jdkmath.JdkMath; 028import org.apache.commons.math4.legacy.linear.RealMatrix; 029import org.apache.commons.math4.legacy.linear.RectangularCholeskyDecomposition; 030 031/** 032 * Generates vectors with with correlated components. 033 * 034 * <p>Random vectors with correlated components are built by combining 035 * the uncorrelated components of another random vector in such a way 036 * that the resulting correlations are the ones specified by a positive 037 * definite covariance matrix.</p> 038 * 039 * <p>The main use of correlated random vector generation is for Monte-Carlo 040 * simulation of physical problems with several variables (for example to 041 * generate error vectors to be added to a nominal vector). A particularly 042 * common case is when the generated vector should be drawn from a 043 * <a href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution"> 044 * Multivariate Normal Distribution</a>, usually using Cholesky decomposition. 045 * Other distributions are possible as long as the underlying sampler provides 046 * normalized values (unit standard deviation).</p> 047 * 048 * <p>Sometimes, the covariance matrix for a given simulation is not 049 * strictly positive definite. This means that the correlations are 050 * not all independent from each other. In this case, however, the non 051 * strictly positive elements found during the Cholesky decomposition 052 * of the covariance matrix should not be negative either, they 053 * should be null. Another non-conventional extension handling this case 054 * is used here. Rather than computing <code>C = U<sup>T</sup> U</code> 055 * where <code>C</code> is the covariance matrix and <code>U</code> 056 * is an upper-triangular matrix, we compute <code>C = B B<sup>T</sup></code> 057 * where <code>B</code> is a rectangular matrix having more rows than 058 * columns. The number of columns of <code>B</code> is the rank of the 059 * covariance matrix, and it is the dimension of the uncorrelated 060 * random vector that is needed to compute the component of the 061 * correlated vector. This class handles this situation automatically.</p> 062 */ 063public class CorrelatedVectorFactory { 064 /** Square root of three. */ 065 private static final double SQRT3 = JdkMath.sqrt(3); 066 /** Mean vector. */ 067 private final double[] mean; 068 /** Root of the covariance matrix. */ 069 private final RealMatrix root; 070 /** Size of uncorrelated vector. */ 071 private final int lengthUncorrelated; 072 /** Size of correlated vector. */ 073 private final int lengthCorrelated; 074 075 /** 076 * Correlated vector factory. 077 * 078 * @param mean Expected mean values of the components. 079 * @param covariance Covariance matrix. 080 * @param small Diagonal elements threshold under which columns are 081 * considered to be dependent on previous ones and are discarded. 082 * @throws org.apache.commons.math4.legacy.linear.NonPositiveDefiniteMatrixException 083 * if the covariance matrix is not strictly positive definite. 084 * @throws DimensionMismatchException if the mean and covariance 085 * arrays dimensions do not match. 086 */ 087 public CorrelatedVectorFactory(double[] mean, 088 RealMatrix covariance, 089 double small) { 090 lengthCorrelated = covariance.getRowDimension(); 091 if (mean.length != lengthCorrelated) { 092 throw new DimensionMismatchException(mean.length, lengthCorrelated); 093 } 094 this.mean = mean.clone(); 095 096 final RectangularCholeskyDecomposition decomposition 097 = new RectangularCholeskyDecomposition(covariance, small); 098 root = decomposition.getRootMatrix(); 099 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}