MultivariateGaussianRandomizer.java
/*
* Copyright (C) 2016 Alberto Irurueta Carro (alberto@irurueta.com)
*
* Licensed 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 com.irurueta.statistics;
import com.irurueta.algebra.AlgebraException;
import com.irurueta.algebra.CholeskyDecomposer;
import com.irurueta.algebra.Matrix;
import com.irurueta.algebra.WrongSizeException;
import java.util.Random;
/**
* Generates pseudo-random values following a multivariate Gaussian distribution
* having the specified mean and covariance. By default, mean is equal to zero
* and the covariance is equal to the identity (unitary standard deviation).
*/
public class MultivariateGaussianRandomizer {
/**
* Contains mean values to be used for random value generation.
*/
private double[] mean;
/**
* Covariance matrix.
*/
private Matrix covariance;
/**
* Lower triangular Cholesky decomposition of covariance matrix.
*/
private Matrix l;
/**
* Instance in charge of generating pseudo-random values. Secure instances
* can be used if the generated values need to be ensured "more" random at
* the expense of higher computational cost.
*/
private final Random internalRandom;
/**
* Constructor.
*/
public MultivariateGaussianRandomizer() {
this(new Random());
}
/**
* Constructor.
* Because neither mean nor covariance are provided, default values will be
* used instead.
*
* @param internalRandom internal random instance in charge of generating
* pseudo-random values.
* @throws NullPointerException thrown if provided internal random is null.
*/
public MultivariateGaussianRandomizer(final Random internalRandom) {
if (internalRandom == null) {
throw new NullPointerException();
}
this.internalRandom = internalRandom;
mean = new double[1];
try {
covariance = l = Matrix.identity(1, 1);
} catch (final WrongSizeException ignore) {
// never thrown
}
}
/**
* Constructor.
*
* @param mean mean to be set.
* @param covariance covariance to be set.
* @throws WrongSizeException if mean length is not compatible with
* covariance size. Mean length must be equal to size of square covariance
* matrix.
* @throws InvalidCovarianceMatrixException if provided covariance matrix is
* not symmetric positive definite.
*/
public MultivariateGaussianRandomizer(final double[] mean, final Matrix covariance)
throws WrongSizeException, InvalidCovarianceMatrixException {
this(new Random(), mean, covariance);
}
/**
* Constructor.
*
* @param internalRandom internal random instance in charge of generating
* pseudo-random values.
* @param mean mean to be set.
* @param covariance covariance to be set.
* @throws WrongSizeException if mean length is not compatible with
* covariance size. Mean length must be equal to size of square covariance
* matrix.
* @throws InvalidCovarianceMatrixException if provided covariance matrix is
* not symmetric positive definite.
* @throws NullPointerException thrown if provided internal random is null.
*/
public MultivariateGaussianRandomizer(
final Random internalRandom, final double[] mean, final Matrix covariance) throws WrongSizeException,
InvalidCovarianceMatrixException {
if (internalRandom == null) {
throw new NullPointerException();
}
this.internalRandom = internalRandom;
setMeanAndCovariance(mean, covariance);
}
/**
* Returns mean value to be used for Gaussian random value generation.
*
* @return mean value.
*/
public double[] getMean() {
return mean;
}
/**
* Returns covariance to be used for Gaussian random value generation.
*
* @return covariance.
*/
public Matrix getCovariance() {
return covariance;
}
/**
* Sets mean and covariance to generate multivariate Gaussian random values.
*
* @param mean mean to be set.
* @param covariance covariance to be set.
* @throws WrongSizeException if mean length is not compatible with
* covariance size. Mean length must be equal to size of square covariance
* matrix.
* @throws InvalidCovarianceMatrixException if provided covariance matrix is
* not symmetric positive definite.
*/
public final void setMeanAndCovariance(final double[] mean, final Matrix covariance)
throws WrongSizeException, InvalidCovarianceMatrixException {
final var length = mean.length;
if (covariance.getRows() != length || covariance.getColumns() != length) {
throw new WrongSizeException("mean must have same covariance size");
}
try {
final var decomposer = new CholeskyDecomposer(covariance);
decomposer.decompose();
if (!decomposer.isSPD()) {
throw new InvalidCovarianceMatrixException(
"covariance matrix must be symmetric positive definite (non singular)");
}
this.mean = mean;
this.covariance = covariance;
l = decomposer.getL();
} catch (final AlgebraException e) {
throw new InvalidCovarianceMatrixException("covariance matrix must be square", e);
}
}
/**
* Generate next set of multivariate Gaussian random values having current
* mean and covariance of this instance.
*
* @param values array where generated random values will be stored.
* @throws IllegalArgumentException if provided array length does not have the same length
* as provided mean array.
*/
public void next(final double[] values) {
final var n = values.length;
if (n != mean.length) {
throw new IllegalArgumentException("values must have mean length");
}
// generate initial Gaussian values with zero mean and unitary standard
// deviation
final var tmp = new double[n];
for (var i = 0; i < n; i++) {
tmp[i] = internalRandom.nextGaussian();
}
// multiply square root of covariance (its Lower triangular Cholesky
// decomposition) by the generated values and add mean
for (var i = 0; i < n; i++) {
values[i] = 0.0;
for (var j = 0; j < n; j++) {
if (i >= j) {
// only evaluate lower triangular portion
values[i] += l.getElementAt(i, j) * tmp[j];
}
}
// add mean
values[i] += mean[i];
}
}
/**
* Generate next set of multivariate Gaussian random values having current
* mean and covariance of this instance.
*
* @return a new array containing generated random values.
*/
public double[] next() {
final var values = new double[mean.length];
next(values);
return values;
}
}