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.math3.random; 019 020import org.apache.commons.math3.exception.DimensionMismatchException; 021import org.apache.commons.math3.linear.RealMatrix; 022import org.apache.commons.math3.linear.RectangularCholeskyDecomposition; 023 024/** 025 * A {@link RandomVectorGenerator} that generates vectors with with 026 * correlated components. 027 * <p>Random vectors with correlated components are built by combining 028 * the uncorrelated components of another random vector in such a way that 029 * the resulting correlations are the ones specified by a positive 030 * definite covariance matrix.</p> 031 * <p>The main use for correlated random vector generation is for Monte-Carlo 032 * simulation of physical problems with several variables, for example to 033 * generate error vectors to be added to a nominal vector. A particularly 034 * interesting case is when the generated vector should be drawn from a <a 035 * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution"> 036 * Multivariate Normal Distribution</a>. The approach using a Cholesky 037 * decomposition is quite usual in this case. However, it can be extended 038 * to other cases as long as the underlying random generator provides 039 * {@link NormalizedRandomGenerator normalized values} like {@link 040 * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p> 041 * <p>Sometimes, the covariance matrix for a given simulation is not 042 * strictly positive definite. This means that the correlations are 043 * not all independent from each other. In this case, however, the non 044 * strictly positive elements found during the Cholesky decomposition 045 * of the covariance matrix should not be negative either, they 046 * should be null. Another non-conventional extension handling this case 047 * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code> 048 * where <code>C</code> is the covariance matrix and <code>U</code> 049 * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code> 050 * where <code>B</code> is a rectangular matrix having 051 * more rows than columns. The number of columns of <code>B</code> is 052 * the rank of the covariance matrix, and it is the dimension of the 053 * uncorrelated random vector that is needed to compute the component 054 * of the correlated vector. This class handles this situation 055 * automatically.</p> 056 * 057 * @since 1.2 058 */ 059 060public class CorrelatedRandomVectorGenerator 061 implements RandomVectorGenerator { 062 /** Mean vector. */ 063 private final double[] mean; 064 /** Underlying generator. */ 065 private final NormalizedRandomGenerator generator; 066 /** Storage for the normalized vector. */ 067 private final double[] normalized; 068 /** Root of the covariance matrix. */ 069 private final RealMatrix root; 070 071 /** 072 * Builds a correlated random vector generator from its mean 073 * vector and covariance matrix. 074 * 075 * @param mean Expected mean values for all components. 076 * @param covariance Covariance matrix. 077 * @param small Diagonal elements threshold under which column are 078 * considered to be dependent on previous ones and are discarded 079 * @param generator underlying generator for uncorrelated normalized 080 * components. 081 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException 082 * if the covariance matrix is not strictly positive definite. 083 * @throws DimensionMismatchException if the mean and covariance 084 * arrays dimensions do not match. 085 */ 086 public CorrelatedRandomVectorGenerator(double[] mean, 087 RealMatrix covariance, double small, 088 NormalizedRandomGenerator generator) { 089 int order = covariance.getRowDimension(); 090 if (mean.length != order) { 091 throw new DimensionMismatchException(mean.length, order); 092 } 093 this.mean = mean.clone(); 094 095 final RectangularCholeskyDecomposition decomposition = 096 new RectangularCholeskyDecomposition(covariance, small); 097 root = decomposition.getRootMatrix(); 098 099 this.generator = generator; 100 normalized = new double[decomposition.getRank()]; 101 102 } 103 104 /** 105 * Builds a null mean random correlated vector generator from its 106 * covariance matrix. 107 * 108 * @param covariance Covariance matrix. 109 * @param small Diagonal elements threshold under which column are 110 * considered to be dependent on previous ones and are discarded. 111 * @param generator Underlying generator for uncorrelated normalized 112 * components. 113 * @throws org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException 114 * if the covariance matrix is not strictly positive definite. 115 */ 116 public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small, 117 NormalizedRandomGenerator generator) { 118 int order = covariance.getRowDimension(); 119 mean = new double[order]; 120 for (int i = 0; i < order; ++i) { 121 mean[i] = 0; 122 } 123 124 final RectangularCholeskyDecomposition decomposition = 125 new RectangularCholeskyDecomposition(covariance, small); 126 root = decomposition.getRootMatrix(); 127 128 this.generator = generator; 129 normalized = new double[decomposition.getRank()]; 130 131 } 132 133 /** Get the underlying normalized components generator. 134 * @return underlying uncorrelated components generator 135 */ 136 public NormalizedRandomGenerator getGenerator() { 137 return generator; 138 } 139 140 /** Get the rank of the covariance matrix. 141 * The rank is the number of independent rows in the covariance 142 * matrix, it is also the number of columns of the root matrix. 143 * @return rank of the square matrix. 144 * @see #getRootMatrix() 145 */ 146 public int getRank() { 147 return normalized.length; 148 } 149 150 /** Get the root of the covariance matrix. 151 * The root is the rectangular matrix <code>B</code> such that 152 * the covariance matrix is equal to <code>B.B<sup>T</sup></code> 153 * @return root of the square matrix 154 * @see #getRank() 155 */ 156 public RealMatrix getRootMatrix() { 157 return root; 158 } 159 160 /** Generate a correlated random vector. 161 * @return a random vector as an array of double. The returned array 162 * is created at each call, the caller can do what it wants with it. 163 */ 164 public double[] nextVector() { 165 166 // generate uncorrelated vector 167 for (int i = 0; i < normalized.length; ++i) { 168 normalized[i] = generator.nextNormalizedDouble(); 169 } 170 171 // compute correlated vector 172 double[] correlated = new double[mean.length]; 173 for (int i = 0; i < correlated.length; ++i) { 174 correlated[i] = mean[i]; 175 for (int j = 0; j < root.getColumnDimension(); ++j) { 176 correlated[i] += root.getEntry(i, j) * normalized[j]; 177 } 178 } 179 180 return correlated; 181 182 } 183 184}