/*
* Class: MultinormalCholeskyGen
* Description: multivariate normal random variable generator
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* A copy of the GNU General Public License is available at
<a href="http://www.gnu.org/licenses">GPL licence site</a>.
*/
package umontreal.iro.lecuyer.randvarmulti;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.CholeskyDecomposition;
import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.randvar.NormalGen;
import umontreal.iro.lecuyer.rng.RandomStream;
/**
* Extends {@link MultinormalGen} for a <SPAN CLASS="textit">multivariate normal</SPAN> distribution, generated via a Cholesky decomposition of the covariance
* matrix. The covariance matrix
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN> is decomposed (by the constructor)
* as
* <SPAN CLASS="MATH"><I><B>Σ</B></I> = <B>A</B><B>A</B><SUP>t</SUP></SPAN> where
* <SPAN CLASS="MATH"><B>A</B></SPAN> is a lower-triangular matrix
* (this is the Cholesky decomposition), and
* <SPAN CLASS="MATH"><B>X</B></SPAN> is generated via
*
* <P></P>
* <DIV ALIGN="CENTER" CLASS="mathdisplay">
* <B>X</B> = <I><B>μ</B></I> + <B>A</B><B>Z</B>,
* </DIV><P></P>
* where
* <SPAN CLASS="MATH"><B>Z</B></SPAN> is a <SPAN CLASS="MATH"><I>d</I></SPAN>-dimensional vector of independent standard normal random
* variates, and
* <SPAN CLASS="MATH"><B>A</B><SUP>t</SUP></SPAN> is the transpose of
* <SPAN CLASS="MATH"><B>A</B></SPAN>.
* The covariance matrix
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN> must be positive-definite, otherwise the
* Cholesky decomposition will fail. The decomposition method
* uses the <TT>CholeskyDecomposition</TT> class in <TT>colt</TT>.
*
*/
public class MultinormalCholeskyGen extends MultinormalGen {
private void initL() {
if (mu.length != sigma.rows() || mu.length != sigma.columns())
throw new IllegalArgumentException
("Incompatible mean vector and covariance matrix");
CholeskyDecomposition decomp = new CholeskyDecomposition (sigma);
//if (!decomp.isSymmetricPositiveDefinite())
// throw new IllegalArgumentException
// ("The covariance matrix must be symmetric and positive-definite");
sqrtSigma = decomp.getL();
}
/**
* Equivalent to
* {@link #MultinormalCholeskyGen((NormalGen, double[], DoubleMatrix2D)) MultinormalCholeskyGen}<TT>(gen1, mu, new DenseDoubleMatrix2D(sigma))</TT>.
*
* @param gen1 the one-dimensional generator
*
* @param mu the mean vector.
*
* @param sigma the covariance matrix.
*
* @exception NullPointerException if any argument is <TT>null</TT>.
*
* @exception IllegalArgumentException if the length of the mean
* vector is incompatible with the dimensions of the covariance matrix.
*
*
*/
public MultinormalCholeskyGen (NormalGen gen1, double[] mu,
double[][] sigma) {
super(gen1, mu, sigma);
initL();
}
/**
* Constructs a multinormal generator with mean vector <TT>mu</TT>
* and covariance matrix <TT>sigma</TT>. 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
* <SPAN CLASS="MATH"><B>Z</B></SPAN> is generated by calling <SPAN CLASS="MATH"><I>d</I></SPAN> times the generator <TT>gen1</TT>,
* which must be a <SPAN CLASS="textit">standard normal</SPAN> 1-dimensional generator.
*
* @param gen1 the one-dimensional generator
*
* @param mu the mean vector.
*
* @param sigma the covariance matrix.
*
* @exception NullPointerException if any argument is <TT>null</TT>.
*
* @exception IllegalArgumentException if the length of the mean
* vector is incompatible with the dimensions of the covariance matrix.
*
*/
public MultinormalCholeskyGen (NormalGen gen1, double[] mu,
DoubleMatrix2D sigma) {
super(gen1, mu, sigma);
initL();
}
/**
* Returns the lower-triangular matrix
* <SPAN CLASS="MATH"><B>A</B></SPAN> in the
* Cholesky decomposition of
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN>.
*
* @return the Cholesky decomposition of the covariance matrix.
*
*/
public DoubleMatrix2D getCholeskyDecompSigma() {
return sqrtSigma.copy();
}
/**
* Sets the covariance matrix
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN> of this multinormal generator
* to <TT>sigma</TT> (and recomputes
* <SPAN CLASS="MATH"><B>A</B></SPAN>).
*
* @param sigma the new covariance matrix.
*
* @exception IllegalArgumentException if <TT>sigma</TT> has
* incorrect dimensions.
*
*
*/
public void setSigma (DoubleMatrix2D sigma) {
if (sigma.rows() != mu.length || sigma.columns() != mu.length)
throw new IllegalArgumentException
("Invalid dimensions of covariance matrix");
CholeskyDecomposition decomp = new CholeskyDecomposition (sigma);
//if (!decomp.isSymmetricPositiveDefinite())
// throw new IllegalArgumentException
// ("The new covariance matrix must be symmetric and positive-definite");
this.sigma.assign (sigma);
this.sqrtSigma = decomp.getL();
}
/**
* Equivalent to
* {@link #nextPoint((NormalGen, double[], DoubleMatrix2D, double[])) nextPoint}<TT>(gen1, mu, new DenseDoubleMatrix2D(sigma), p)</TT>.
*
*/
public static void nextPoint (NormalGen gen1, double[] mu,
double[][] sigma, double[] p) {
nextPoint (gen1, mu, new DenseDoubleMatrix2D (sigma), p);
}
/**
* Generates a <SPAN CLASS="MATH"><I>d</I></SPAN>-dimensional vector from the multinormal
* distribution with mean vector <TT>mu</TT> and covariance matrix
* <TT>sigma</TT>, using the one-dimensional normal generator <TT>gen1</TT> to
* generate the coordinates of
* <SPAN CLASS="MATH"><B>Z</B></SPAN>, and using the Cholesky decomposition of
*
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN>. The resulting vector is put into <TT>p</TT>.
* Note that this static method will be very slow for large dimensions, since
* it computes the Cholesky decomposition at every call. It is therefore
* recommended to use a <TT>MultinormalCholeskyGen</TT> object instead,
* if the method is to be called more than once.
*
* @param p the array to be filled with the generated point.
*
* @exception IllegalArgumentException if the one-dimensional normal
* generator uses a normal distribution with <SPAN CLASS="MATH"><I>μ</I></SPAN> not equal to 0, or
* <SPAN CLASS="MATH"><I>σ</I></SPAN> not equal to 1.
*
* @exception IllegalArgumentException if the length of the mean
* vector is different from the dimensions of the covariance matrix,
* or if the covariance matrix is not symmetric and positive-definite.
*
* @exception NullPointerException if any argument is <TT>null</TT>.
*
*
*/
public static void nextPoint (NormalGen gen1, double[] mu,
DoubleMatrix2D sigma, double[] p) {
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");
if (mu.length != sigma.rows() ||
mu.length != sigma.columns())
throw new IllegalArgumentException
("Incompatible mean vector and covariance matrix dimensions");
CholeskyDecomposition decomp = new CholeskyDecomposition (sigma);
double[] temp = new double[mu.length];
DoubleMatrix2D sqrtSigma = decomp.getL();
for (int i = 0; i < temp.length; 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 < temp.length; i++) {
p[i] = 0;
for (int c = 0; c <= i; c++)
p[i] += sqrtSigma.getQuick (i, c)*temp[c];
p[i] += mu[i];
}
}
/**
* Generates a point from this multinormal distribution. This is much
* faster than the static method as it computes the singular value decomposition
* matrix only once in the constructor.
*
* @param p the array to be filled with the generated point
*
*/
public void nextPoint (double[] p) {
int n = mu.length;
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;
// Matrix is lower-triangular
for (int c = 0; c <= i; c++)
p[i] += sqrtSigma.getQuick (i, c)*temp[c];
p[i] += mu[i];
}
}
}