package umontreal.iro.lecuyer.randvarmulti;
import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.randvar.NormalGen;
import umontreal.iro.lecuyer.rng.RandomStream;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.linalg.CholeskyDecomposition;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
/**
* Extends {@link RandomMultivariateGen} for a
* multivariate normal (or multinormal) distribution.
* The d-dimensional multivariate normal distribution
* with mean vector
* μ∈Rd and (symmetric positive-definite)
* covariance matrix
* Σ, denoted
* N(μ, Σ), has density
*
*
*
* f (X) = exp(- (X - μ)tΣ-1(X - μ)/2)/((2π)^d )1/2,
*
* for all
* X∈Rd, and
* Xt is the transpose vector of
* X.
* If
* Z∼N( 0,I) where
* I is the
* identity matrix,
* Z is said to have the standard multinormal
* distribution.
*
*
* For the special case d = 2, if the random vector
* X = (X1, X2)t
* has a bivariate normal distribution, then it has mean
*
* μ = (μ1, μ2)t, and covariance matrix
*
*
*
* Σ = [
]
*
* if and only if
* Var[X1] = σ12,
* Var[X2] = σ22, and the
* linear correlation between X1 and X2 is ρ, where
* σ1 > 0,
*
* σ2 > 0, and
* -1 <= ρ <= 1.
*
*/
public class MultinormalGen extends RandomMultivariateGen {
protected double[] mu;
protected DoubleMatrix2D sigma;
protected DoubleMatrix2D sqrtSigma;
protected double[] temp;
protected static final double MYINF = 37.54;
private void initMN (NormalGen gen1, double[] mu, int d) {
if (gen1 == null)
throw new NullPointerException ("gen1 is null");
NormalDist dist = (NormalDist) gen1.getDistribution();
if (dist.getMu() != 0.0)
throw new IllegalArgumentException ("mu != 0");
if (dist.getSigma() != 1.0)
throw new IllegalArgumentException ("sigma != 1");
this.gen1 = gen1;
dist = null;
if (mu == null) { // d is the dimension
dimension = d;
this.mu = new double[d];
} else { // d is unused
dimension = mu.length;
this.mu = (double[])mu.clone();
}
temp = new double[dimension];
}
/**
* Constructs a generator with the standard multinormal distribution
* (with
* μ = 0 and
* Σ = I) in d dimensions.
* Each vector
* Z will be generated via d successive calls to
* gen1, which must be a standard normal generator.
*
* @param gen1 the one-dimensional generator
*
* @param d the dimension of the generated vectors
*
* @exception IllegalArgumentException if the one-dimensional normal
* generator uses a normal distribution with μ not equal to 0, or
* σ not equal to 1.
*
* @exception IllegalArgumentException if d
* is negative.
*
* @exception NullPointerException if gen1 is null.
*
*
*/
public MultinormalGen (NormalGen gen1, int d) {
initMN (gen1, null, d);
sigma = new DenseDoubleMatrix2D (d, d);
sqrtSigma = new DenseDoubleMatrix2D (d, d);
for (int i = 0; i < d; i++) {
sigma.setQuick (i, i, 1.0);
sqrtSigma.setQuick (i, i, 1.0);
}
}
/**
* Constructs a multinormal generator with mean vector
* mu and covariance matrix sigma.
* The mean vector must have the same length as the dimensions
* of the covariance matrix, which must be symmetric and positive-definite.
* If any of the above conditions is violated, an exception is thrown.
* The vector
* Z is generated by calling d times the generator gen1,
* which must be standard normal.
*
* @param gen1 the one-dimensional generator
*
* @param mu the mean vector.
*
* @param sigma the covariance matrix.
*
* @exception NullPointerException if any argument is null.
*
* @exception IllegalArgumentException if the length of the mean
* vector is incompatible with the dimensions of the covariance matrix.
*
*
*/
protected MultinormalGen (NormalGen gen1, double[] mu,
DoubleMatrix2D sigma) {
initMN (gen1, mu, -1);
this.sigma = sigma.copy();
}
/**
* Equivalent to
* {@link #MultinormalGen((NormalGen, double[], DoubleMatrix2D)) MultinormalGen} (gen1, mu, new DenseDoubleMatrix2D (sigma)).
*
*/
protected MultinormalGen (NormalGen gen1, double[] mu, double[][] sigma) {
initMN (gen1, mu, -1);
this.sigma = new DenseDoubleMatrix2D (sigma);
}
/**
* Returns the mean vector used by this generator.
*
* @return the current mean vector.
*
*/
public double[] getMu() {
return mu;
}
/**
* Returns the i-th component of the mean vector
* for this generator.
*
* @param i the index of the required component.
*
* @return the value of μi.
* @exception ArrayIndexOutOfBoundsException if
* i is negative or greater than or equal to {@link #getDimension(()) getDimension}.
*
*
*/
public double getMu (int i) {
return mu[i];
}
/**
* Sets the mean vector to mu.
*
* @param mu the new mean vector.
*
* @exception NullPointerException if mu is null.
*
* @exception IllegalArgumentException if the length of mu
* does not correspond to {@link #getDimension(()) getDimension}.
*
*
*/
public void setMu (double[] mu) {
if (mu.length != this.mu.length)
throw new IllegalArgumentException
("Incompatible length of mean vector");
this.mu = mu;
}
/**
* Sets the i-th component of the mean vector to mui.
*
* @param i the index of the modified component.
*
* @param mui the new value of μi.
*
* @exception ArrayIndexOutOfBoundsException if i
* is negative or greater than or equal to {@link #getDimension(()) getDimension}.
*
*
*/
public void setMu (int i, double mui) {
mu[i] = mui;
}
/**
* Returns the covariance matrix
* Σ
* used by this generator.
*
* @return the used covariance matrix.
*
*/
public DoubleMatrix2D getSigma() {
return sigma.copy();
}
/**
* Generates a point from this multinormal distribution.
*
* @param p the array to be filled with the generated point
*
*
*/
public void nextPoint (double[] p) {
int n = dimension;
for (int i = 0; i < n; i++) {
temp[i] = gen1.nextDouble();
if (temp[i] == Double.NEGATIVE_INFINITY)
temp[i] = -MYINF;
if (temp[i] == Double.POSITIVE_INFINITY)
temp[i] = MYINF;
}
for (int i = 0; i < n; i++) {
p[i] = 0;
for (int c = 0; c < n; c++)
p[i] += sqrtSigma.getQuick (i, c)*temp[c];
p[i] += mu[i];
}
}
}