MultinormalPCAGen.java

/*
 * Class:        MultinormalPCAGen
 * Description:  multivariate normal random variable generator
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * 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 umontreal.ssj.randvarmulti;
 
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.linalg.SingularValueDecomposition;
import umontreal.ssj.probdist.NormalDist;
import umontreal.ssj.randvar.NormalGen;
import umontreal.ssj.rng.RandomStream;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;

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));
   }

   public MultinormalPCAGen(NormalGen gen1, double[] mu, double[][] sigma) {
      super(gen1, mu, sigma);
      initL();
   }

   public MultinormalPCAGen(NormalGen gen1, double[] mu, DoubleMatrix2D sigma) {
      super(gen1, mu, sigma);
      initL();
   }

   public static DoubleMatrix2D decompPCA(double[][] sigma) {
      return decompPCA(new DenseDoubleMatrix2D(sigma));
   }

   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);
   }

   public DoubleMatrix2D getPCADecompSigma() {
      return sqrtSigma.copy();
   }

   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;
   }

   public double[] getLambda() {
      return lambda;
   }

   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));
   }

   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];
      }
   }

   public static void nextPoint(NormalGen gen1, double[] mu, double[][] sigma, double[] p) {
      nextPoint(gen1, mu, new DenseDoubleMatrix2D(sigma), p);
   }

   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];
      }
   }
}