umontreal.iro.lecuyer.probdistmulti.MultiNormalDist Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: MultiNormalDist
* Description: multinormal distribution
* 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
GPL licence site.
*/
package umontreal.iro.lecuyer.probdistmulti;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;
/**
* Implements the abstract class {@link ContinuousDistributionMulti} for the
* multinormal distribution with mean vector μ and covariance
* matrix
* Σ.
* The probability density is
*
*
*
* f (x = x1,…, xd) = exp(- (x - μ)TΣ-1(x - μ)/2)/((2π)^d det())1/2
*
* where
* 
= (x1,…, xd).
*
*/
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());
}
/**
* Computes the density of the multinormal distribution
* with parameters μ = mu and
* Σ = sigma,
* evaluated at x.
*
*/
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)));
}
/**
* Returns the dimension d of the distribution.
*
*/
public int getDimension() {
return dim;
}
/**
* Returns the mean
* E[X] = μ of the multinormal distribution
* with parameters μ and
* Σ.
*
*/
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;
}
/**
* Computes the covariance matrix of the multinormal distribution
* with parameters μ and
* Σ.
*
*/
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;
}
/**
* Computes the correlation matrix of the multinormal distribution
* with parameters μ and
* Σ).
*
*/
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);
}
/**
* Estimates the parameters μ of the multinormal distribution using
* the maximum likelihood method. It uses the n observations of d
* components in table x[i][j],
* i = 0, 1,…, n - 1 and
*
* j = 0, 1,…, d - 1.
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @param d the dimension of each observation
*
* @return returns the parameters [μ1,...,μd]
*
*/
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;
}
/**
* Estimates the parameters
* Σ of the multinormal distribution using
* the maximum likelihood method. It uses the n observations of d
* components in table x[i][j],
* i = 0, 1,…, n - 1 and
*
* j = 0, 1,…, d - 1.
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @param d the dimension of each observation
*
* @return returns the covariance matrix
* Σ
*
*/
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;
}
/**
* Returns the parameter μ of this object.
*
*/
public double[] getMu() {
return mu;
}
/**
* Returns the i-th component of the parameter μ of this object.
*
*/
public double getMu (int i) {
return mu[i];
}
/**
* Returns the parameter
* Σ of this object.
*
*/
public double[][] getSigma() {
return sigma.toArray();
}
/**
* Sets the parameters μ and
* Σ of this object.
*
*/
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;
}
}