MultinormalCholeskyGen.java

/*
 * Class:        MultinormalCholeskyGen
 * 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.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.CholeskyDecomposition;
import umontreal.ssj.probdist.NormalDist;
import umontreal.ssj.randvar.NormalGen;
import umontreal.ssj.rng.RandomStream;

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

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

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

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

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

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

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

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