/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
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
*
*
*
* P[X = (x1,..., xd)] = exp(- (x - μ)TΣ-1(x - μ)/2)/((2π)^d |Σ|)1/2
*
*
*/
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);
}
/**
* It is now called getMLEMu.
*
*/
@Deprecated
public static double[] getMaximumLikelihoodEstimateMu (double[][] x,
int n, int d) {
return getMLEMu(x, n, d);
}
/**
* 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;
}
/**
* It is now called getMLESigma.
*
*/
@Deprecated
public static double[][] getMaximumLikelihoodEstimateSigma (double[][] x,
int n, int d) {
return getMLESigma(x, n, d);
}
/**
* 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 = getMaximumLikelihoodEstimateMu (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;
}
}
Ads help maintain this website.