/*
* Class: MultinormalPCAGen
* 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.linalg.SingularValueDecomposition;
import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.randvar.NormalGen;
import umontreal.iro.lecuyer.rng.RandomStream;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
/**
* Extends {@link MultinormalGen} for a <SPAN CLASS="textit">multivariate normal</SPAN> distribution, generated via the method of principal components analysis
* (PCA) 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>V</B><I><B>Λ</B></I><B>V</B><SUP>t</SUP></SPAN> where
* <SPAN CLASS="MATH"><B>V</B></SPAN> is an orthogonal matrix and
* <SPAN CLASS="MATH"><I><B>Λ</B></I></SPAN> is the diagonal matrix made up
* of the eigenvalues of
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN>. <SPAN CLASS="MATH"><B>V</B><SUP>t</SUP></SPAN> is the transpose
* matrix of <SPAN CLASS="MATH"><B>V</B></SPAN>. The eigenvalues are ordered from the
* largest (<SPAN CLASS="MATH"><I>λ</I><SUB>1</SUB></SPAN>) to the smallest (<SPAN CLASS="MATH"><I>λ</I><SUB>d</SUB></SPAN>). The random multinormal
* vector
* <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>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>, and
* <SPAN CLASS="MATH"><B>Z</B></SPAN> is a <SPAN CLASS="MATH"><I>d</I></SPAN>-dimensional vector of
* independent standard normal random variates. The decomposition method
* uses the <TT>SingularValueDecomposition</TT> class in <TT>colt</TT>.
*
*/
public class MultinormalPCAGen extends MultinormalGen {
private double[] lambda;
private static SingularValueDecomposition getSvd (DoubleMatrix2D sigma) {
return (new SingularValueDecomposition (sigma));
}
private DoubleMatrix2D decompPCA (SingularValueDecomposition svd) {
DoubleMatrix2D D = svd.getS ();
// Calculer la racine carree des valeurs propres
for (int i = 0; i < D.rows(); ++i) {
lambda[i] = D.getQuick (i, i);
D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
}
DoubleMatrix2D P = svd.getV();
return P.zMult (D, null);
}
private void initL() {
if (mu.length != sigma.rows() || mu.length != sigma.columns())
throw new IllegalArgumentException
("Incompatible mean vector and covariance matrix");
lambda = new double[mu.length];
sqrtSigma = decompPCA (getSvd(sigma));
}
/**
* Equivalent to
* {@link #MultinormalPCAGen((NormalGen, double[], DoubleMatrix2D)) MultinormalPCAGen}<TT>(gen1, mu, new DenseDoubleMatrix2D(sigma))</TT>.
*
*/
public MultinormalPCAGen (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 semi-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 MultinormalPCAGen (NormalGen gen1, double[] mu,
DoubleMatrix2D sigma) {
super(gen1, mu, sigma);
initL();
}
/**
* Computes the decomposition <TT>sigma</TT> =
*
* <SPAN CLASS="MATH"><I><B>Σ</B></I> = <B>V</B><I><B>Λ</B></I><B>V</B><SUP>t</SUP></SPAN>. Returns
* <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>.
*
*/
public static DoubleMatrix2D decompPCA (double[][] sigma) {
return decompPCA (new DenseDoubleMatrix2D (sigma));
}
/**
* Computes the decomposition <TT>sigma</TT> =
* <SPAN CLASS="MATH"><I><B>Σ</B></I> = <B>V</B><I><B>Λ</B></I><B>V</B><SUP>t</SUP></SPAN>. Returns
* <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN>.
*
*/
public static DoubleMatrix2D decompPCA (DoubleMatrix2D sigma) {
// L'objet SingularValueDecomposition permet de recuperer la matrice
// des valeurs propres en ordre decroissant et celle des vecteurs propres de
// sigma (pour une matrice symetrique et definie-positive seulement).
SingularValueDecomposition sv = getSvd (sigma);
// D contient les valeurs propres sur la diagonale
DoubleMatrix2D D = sv.getS ();
// Calculer la racine carree des valeurs propres
for (int i = 0; i < D.rows(); ++i)
D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
DoubleMatrix2D P = sv.getV();
// Multiplier par la matrice orthogonale (ici P)
return P.zMult (D, null);
}
/**
* Returns the matrix
* <SPAN CLASS="MATH"><B>A</B> = <B>V</B>()<SUP>1/2</SUP></SPAN> of this object.
*
* @return the PCA square root of the covariance matrix
*
*/
public DoubleMatrix2D getPCADecompSigma() {
return sqrtSigma.copy();
}
/**
* Computes and returns the eigenvalues of <TT>sigma</TT> in
* decreasing order.
*
*/
public static double[] getLambda (DoubleMatrix2D sigma) {
SingularValueDecomposition sv = getSvd (sigma);
DoubleMatrix2D D = sv.getS ();
double[] lambd = new double[D.rows()];
for (int i = 0; i < D.rows(); ++i)
lambd[i] = D.getQuick (i, i);
return lambd;
}
/**
* Returns the eigenvalues of
* <SPAN CLASS="MATH"><I><B>Σ</B></I></SPAN> in decreasing order.
*
*/
public double[] getLambda() {
return lambda;
}
/**
* 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");
this.sigma.assign (sigma);
this.sqrtSigma = decompPCA(getSvd (sigma));
}
/**
* 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 PCA 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, because
* it recomputes the singular value decomposition at every call. It is therefore
* recommended to use a <TT>MultinormalPCAGen</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");
double[] temp = new double[mu.length];
DoubleMatrix2D sqrtSigma = decompPCA(sigma);
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 < temp.length; c++)
p[i] += sqrtSigma.getQuick (i, c)*temp[c];
p[i] += mu[i];
}
}
/**
* 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 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;
for (int c = 0; c < n; c++)
p[i] += sqrtSigma.getQuick (i, c)*temp[c];
p[i] += mu[i];
}
}
}