/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdistmulti;
import umontreal.iro.lecuyer.util.Num;
import optimization.*;
/**
* Implements the abstract class {@link ContinuousDistributionMulti} for the
* Dirichlet distribution with parameters
* (α1,...,αd),
* αi > 0.
* The probability density is
*
*
*
* P[X = (x1,..., xd)] = Γ(α0)∏i=1dxiαi-1/(∏i=1dΓ(αi))
*
* where xi >= 0,
* ∑i=1dxi = 1,
* α0 = ∑i=1dαi,
* and Γ is the Gamma function.
*
*/
public class DirichletDist extends ContinuousDistributionMulti {
private static final double LOGMIN = -709.1; // Log(MIN_DOUBLE/2)
protected double[] alpha;
private static class Optim implements Uncmin_methods
{
double[] logP;
int n;
int k;
public Optim (double[] logP, int n) {
this.n = n;
this.k = logP.length;
this.logP = new double[k];
System.arraycopy (logP, 0, this.logP, 0, k);
}
public double f_to_minimize (double[] alpha) {
double sumAlpha = 0.0;
double sumLnGammaAlpha = 0.0;
double sumAlphaLnP = 0.0;
for (int i = 1; i < alpha.length; i++) {
if (alpha[i] <= 0.0)
return 1.0e200;
sumAlpha += alpha[i];
sumLnGammaAlpha += Num.lnGamma (alpha[i]);
sumAlphaLnP += ((alpha[i] - 1.0) * logP[i - 1]);
}
return (- n * (Num.lnGamma (sumAlpha) - sumLnGammaAlpha + sumAlphaLnP));
}
public void gradient (double[] alpha, double[] g)
{
}
public void hessian (double[] alpha, double[][] h)
{
}
}
public DirichletDist (double[] alpha) {
setParams (alpha);
}
public double density (double[] x) {
return density_ (alpha, x);
}
public double[] getMean() {
return getMean_ (alpha);
}
public double[][] getCovariance() {
return getCovariance_ (alpha);
}
public double[][] getCorrelation () {
return getCorrelation_ (alpha);
}
private static void verifParam (double[] alpha) {
for (int i = 0; i < alpha.length;i++)
{
if (alpha[i] <= 0)
throw new IllegalArgumentException("alpha[" + i + "] <= 0");
}
}
private static double density_ (double[] alpha, double[] x) {
double alpha0 = 0.0;
double sumLnGamma = 0.0;
double sumAlphaLnXi = 0.0;
if (alpha.length != x.length)
throw new IllegalArgumentException ("alpha and x must have the same dimension");
for (int i = 0; i < alpha.length; i++) {
alpha0 += alpha[i];
sumLnGamma += Num.lnGamma (alpha[i]);
if (x[i] <= 0.0 || x[i] >= 1.0)
return 0.0;
sumAlphaLnXi += (alpha[i] - 1.0) * Math.log (x[i]);
}
return Math.exp (Num.lnGamma (alpha0) - sumLnGamma + sumAlphaLnXi);
}
/**
* Computes the density of the Dirichlet distribution
* with parameters (α1,...,αd).
*
*/
public static double density (double[] alpha, double[] x) {
verifParam (alpha);
return density_ (alpha, x);
}
private static double[][] getCovariance_ (double[] alpha) {
double[][] cov = new double[alpha.length][alpha.length];
double alpha0 = 0.0;
for (int i =0; i < alpha.length; i++)
alpha0 += alpha[i];
for (int i = 0; i < alpha.length; i++) {
for (int j = 0; j < alpha.length; j++)
cov[i][j] = - (alpha[i] * alpha[j]) / (alpha0 * alpha0 * (alpha0 + 1.0));
cov[i][i] = (alpha[i] / alpha0) * (1.0 - alpha[i] / alpha0) / (alpha0 + 1.0);
}
return cov;
}
/**
* Computes the covariance matrix of the Dirichlet distribution
* with parameters (α1,...,αd).
*
*/
public static double[][] getCovariance (double[] alpha) {
verifParam (alpha);
return getCovariance_ (alpha);
}
private static double[][] getCorrelation_ (double[] alpha) {
double corr[][] = new double[alpha.length][alpha.length];
double alpha0 = 0.0;
for (int i =0; i < alpha.length; i++)
alpha0 += alpha[i];
for (int i = 0; i < alpha.length; i++) {
for (int j = 0; j < alpha.length; j++)
corr[i][j] = - Math.sqrt ((alpha[i] * alpha[j]) /
((alpha0 - alpha[i]) * (alpha0 - alpha[j])));
corr[i][i] = 1.0;
}
return corr;
}
/**
* Computes the correlation matrix of the Dirichlet distribution
* with parameters (α1,...,αd).
*
*/
public static double[][] getCorrelation (double[] alpha) {
verifParam (alpha);
return getCorrelation_ (alpha);
}
/**
* It is now called getMLE.
*
*/
@Deprecated
public static double[] getMaximumLikelihoodEstimate (double[][] x,
int n, int d) {
return getMLE(x, n, d);
}
/**
* Estimates the parameters [
* hat(α_1),…, hat(α_d)]
* of the Dirichlet 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 to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @param d the dimension of each vector
*
* @return returns the parameter [
* hat(α_1),…, hat(α_d)]
*
*/
public static double[] getMLE (double[][] x, int n, int d) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
if (d <= 0)
throw new IllegalArgumentException ("d <= 0");
double[] logP = new double[d];
double mean[] = new double[d];
double var[] = new double[d];
int i;
int j;
for (i = 0; i < d; i++) {
logP[i] = 0.0;
mean[i] = 0.0;
}
for (i = 0; i < n; i++) {
for (j = 0; j < d; j++) {
if (x[i][j] > 0.)
logP[j] += Math.log (x[i][j]);
else
logP[j] += LOGMIN;
mean[j] += x[i][j];
}
}
for (i = 0; i < d; i++) {
logP[i] /= (double) n;
mean[i] /= (double) n;
}
double sum = 0.0;
for (j = 0; j < d; j++) {
sum = 0.0;
for (i = 0; i < n; i++)
sum += (x[i][j] - mean[j]) * (x[i][j] - mean[j]);
var[j] = sum / (double) n;
}
double alpha0 = (mean[0] * (1.0 - mean[0])) / var[0] - 1.0;
Optim system = new Optim (logP, n);
double[] parameters = new double[d];
double[] xpls = new double[d + 1];
double[] alpha = new double[d + 1];
double[] fpls = new double[d + 1];
double[] gpls = new double[d + 1];
int[] itrcmd = new int[2];
double[][] a = new double[d + 1][d + 1];
double[] udiag = new double[d + 1];
for (i = 1; i <= d; i++)
alpha[i] = mean[i - 1] * alpha0;
Uncmin_f77.optif0_f77 (d, alpha, system, xpls, fpls, gpls, itrcmd, a, udiag);
for (i = 0; i < d; i++)
parameters[i] = xpls[i+1];
return parameters;
}
private static double[] getMean_ (double[] alpha) {
double alpha0 = 0.0;
double[] mean = new double[alpha.length];
for (int i = 0; i < alpha.length;i++)
alpha0 += alpha[i];
for (int i = 0; i < alpha.length; i++)
mean[i] = alpha[i] / alpha0;
return mean;
}
/**
* Computes the mean
* E[X] = αi/α0 of the Dirichlet distribution
* with parameters (α1,...,αd), where
* α0 = ∑i=1dαi.
*
*/
public static double[] getMean (double[] alpha) {
verifParam (alpha);
return getMean_ (alpha);
}
/**
* Returns the parameters (α1,...,αd) of this object.
*
*/
public double[] getAlpha() {
return alpha;
}
/**
* Returns the ith component of the alpha vector.
*
*/
public double getAlpha (int i) {
return alpha[i];
}
/**
* Sets the parameters (α1,...,αd) of this object.
*
*/
public void setParams (double[] alpha) {
this.dimension = alpha.length;
this.alpha = new double[dimension];
for (int i = 0; i < dimension; i++) {
if (alpha[i] <= 0)
throw new IllegalArgumentException("alpha[" + i + "] <= 0");
this.alpha[i] = alpha[i];
}
}
}
Ads help maintain this website.