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 multivariate normal distribution, generated via the method of principal components analysis
* (PCA) of the covariance matrix. The covariance matrix
* Σ is
* decomposed (by the constructor) as
* Σ = VΛVt where
* V is an orthogonal matrix and
* Λ is the diagonal matrix made up
* of the eigenvalues of
* Σ. Vt is the transpose
* matrix of V. The eigenvalues are ordered from the
* largest (λ1) to the smallest (λd). The random multinormal
* vector
* X is generated via
*
*
*
* X = μ + AZ,
*
* where
* A = V()1/2, and
* Z is a d-dimensional vector of
* independent standard normal random variates. The decomposition method
* uses the SingularValueDecomposition class in colt.
*
*/
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}(gen1, mu, new DenseDoubleMatrix2D(sigma)).
*
*/
public MultinormalPCAGen (NormalGen gen1, double[] mu, double[][] sigma) {
super(gen1, mu, sigma);
initL();
}
/**
* 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 semi-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 a standard normal 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 null.
*
* @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 sigma =
*
* Σ = VΛVt. Returns
* A = V()1/2.
*
*/
public static DoubleMatrix2D decompPCA (double[][] sigma) {
return decompPCA (new DenseDoubleMatrix2D (sigma));
}
/**
* Computes the decomposition sigma =
* Σ = VΛVt. Returns
* A = V()1/2.
*
*/
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
* A = V()1/2 of this object.
*
* @return the PCA square root of the covariance matrix
*
*/
public DoubleMatrix2D getPCADecompSigma() {
return sqrtSigma.copy();
}
/**
* Computes and returns the eigenvalues of sigma 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
* Σ in decreasing order.
*
*/
public double[] getLambda() {
return lambda;
}
/**
* Sets the covariance matrix
* Σ of this multinormal generator
* to sigma (and recomputes
* A).
*
* @param sigma the new covariance matrix.
*
* @exception IllegalArgumentException if sigma 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 d-dimensional vector from the multinormal
* distribution with mean vector mu and covariance matrix
* sigma, using the one-dimensional normal generator gen1 to
* generate the coordinates of
* Z, and using the PCA decomposition of
*
* Σ. The resulting vector is put into p.
* 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 MultinormalPCAGen 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 μ not equal to 0, or
* σ 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 null.
*
*
*/
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}(gen1, mu, new DenseDoubleMatrix2D(sigma), p).
*
*/
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];
}
}
}
Ads help maintain this website.