MultiNormalDist.java

/*
 * Class:        MultiNormalDist
 * Description:  multinormal distribution
 * 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.probdistmulti;
 
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;

public class MultiNormalDist extends ContinuousDistributionMulti {
   protected int dim;
   protected double[] mu;
   protected DoubleMatrix2D sigma;
   protected DoubleMatrix2D invSigma;
 
   protected static Algebra algebra = new Algebra();
 
   public MultiNormalDist(double[] mu, double[][] sigma) {
      setParams(mu, sigma);
   }
 
   public double density(double[] x) {
      double sum = 0.0;
 
      if (invSigma == null)
         invSigma = algebra.inverse(sigma);
 
      double[] temp = new double[mu.length];
      for (int i = 0; i < mu.length; i++) {
         sum = 0.0;
         for (int j = 0; j < mu.length; j++)
            sum += ((x[j] - mu[j]) * invSigma.getQuick(j, i));
         temp[i] = sum;
      }
 
      sum = 0.0;
      for (int i = 0; i < mu.length; i++)
         sum += temp[i] * (x[i] - mu[i]);
 
      return (Math.exp(-0.5 * sum) / Math.sqrt(Math.pow(2 * Math.PI, mu.length) * algebra.det(sigma)));
   }
 
   public double[] getMean() {
      return mu;
   }
 
   public double[][] getCovariance() {
      return sigma.toArray();
   }
 
   public double[][] getCorrelation() {
      return getCorrelation_(mu, sigma.toArray());
   }

   public static double density(double[] mu, double[][] sigma, double[] x) {
      double sum = 0.0;
      DoubleMatrix2D sig;
      DoubleMatrix2D inv;
 
      if (sigma.length != sigma[0].length)
         throw new IllegalArgumentException("sigma must be a square matrix");
      if (mu.length != sigma.length)
         throw new IllegalArgumentException("mu and sigma must have the same dimension");
 
      sig = new DenseDoubleMatrix2D(sigma);
      inv = algebra.inverse(sig);
 
      double[] temp = new double[mu.length];
      for (int i = 0; i < mu.length; i++) {
         sum = 0.0;
         for (int j = 0; j < mu.length; j++)
            sum += ((x[j] - mu[j]) * inv.getQuick(j, i));
         temp[i] = sum;
      }
 
      sum = 0.0;
      for (int i = 0; i < mu.length; i++)
         sum += temp[i] * (x[i] - mu[i]);
 
      return (Math.exp(-0.5 * sum) / Math.sqrt(Math.pow(2 * Math.PI, mu.length) * algebra.det(sig)));
   }

   public int getDimension() {
      return dim;
   }

   public static double[] getMean(double[] mu, double[][] sigma) {
      if (sigma.length != sigma[0].length)
         throw new IllegalArgumentException("sigma must be a square matrix");
      if (mu.length != sigma.length)
         throw new IllegalArgumentException("mu and sigma must have the same dimension");
 
      return mu;
   }

   public static double[][] getCovariance(double[] mu, double[][] sigma) {
      if (sigma.length != sigma[0].length)
         throw new IllegalArgumentException("sigma must be a square matrix");
      if (mu.length != sigma.length)
         throw new IllegalArgumentException("mu and sigma must have the same dimension");
 
      return sigma;
   }
 
   private static double[][] getCorrelation_(double[] mu, double[][] sigma) {
      double corr[][] = new double[mu.length][mu.length];
 
      for (int i = 0; i < mu.length; i++) {
         for (int j = 0; j < mu.length; j++)
            corr[i][j] = -sigma[i][j] / Math.sqrt(sigma[i][i] * sigma[j][j]);
         corr[i][i] = 1.0;
      }
      return corr;
   }

   public static double[][] getCorrelation(double[] mu, double[][] sigma) {
      if (sigma.length != sigma[0].length)
         throw new IllegalArgumentException("sigma must be a square matrix");
      if (mu.length != sigma.length)
         throw new IllegalArgumentException("mu and sigma must have the same dimension");
 
      return getCorrelation_(mu, sigma);
   }

   public static double[] getMLEMu(double[][] x, int n, int d) {
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
 
      double[] parameters = new double[d];
      for (int i = 0; i < parameters.length; i++)
         parameters[i] = 0.0;
 
      for (int i = 0; i < n; i++)
         for (int j = 0; j < d; j++)
            parameters[j] += x[i][j];
 
      for (int i = 0; i < parameters.length; i++)
         parameters[i] = parameters[i] / (double) n;
 
      return parameters;
   }

   public static double[][] getMLESigma(double[][] x, int n, int d) {
      double sum = 0.0;
 
      if (n <= 0)
         throw new IllegalArgumentException("n <= 0");
      if (d <= 0)
         throw new IllegalArgumentException("d <= 0");
 
      double[] mean = getMLEMu(x, n, d);
      double[][] parameters = new double[d][d];
      for (int i = 0; i < parameters.length; i++)
         for (int j = 0; j < parameters.length; j++)
            parameters[i][j] = 0.0;
 
      for (int i = 0; i < parameters.length; i++) {
         for (int j = 0; j < parameters.length; j++) {
            sum = 0.0;
            for (int t = 0; t < n; t++)
               sum += (x[t][i] - mean[i]) * (x[t][j] - mean[j]);
            parameters[i][j] = sum / (double) n;
         }
      }
 
      return parameters;
   }

   public double[] getMu() {
      return mu;
   }

   public double getMu(int i) {
      return mu[i];
   }

   public double[][] getSigma() {
      return sigma.toArray();
   }

   public void setParams(double[] mu, double[][] sigma) {
      if (sigma.length != sigma[0].length)
         throw new IllegalArgumentException("sigma must be a square matrix");
      if (mu.length != sigma.length)
         throw new IllegalArgumentException("mu and sigma must have the same dimension");
 
      this.mu = new double[mu.length];
      this.dimension = mu.length;
      System.arraycopy(mu, 0, this.mu, 0, mu.length);
      this.sigma = new DenseDoubleMatrix2D(sigma);
 
      invSigma = null;
   }
 
}