/* *
* Title:
* Description:
* Copyright: Copyright (c)
* Company:
* @author
* @version 1.0
*/
package umontreal.iro.lecuyer.probdistmulti;
import umontreal.iro.lecuyer.util.Num;
import umontreal.iro.lecuyer.util.RootFinder;
import umontreal.iro.lecuyer.functions.MathFunction;
/**
* Implements the abstract class {@link DiscreteDistributionIntMulti} for the
* negative multinomial distribution with parameters
* γ > 0 and
* (
* p1,…, pd), such that all 0 < pi < 1 and
* ∑i=1dpi < 1.
* The probability mass function is
*
*
*
* P[X = (x1,…, xd)] = (Γ(γ+∑i=1dxi)/Γ(γ)∏i=1dxi!)p0γ∏i=1dpixi
*
* where
* p0 = 1 - ∑i=1dpi.
*
*/
public class NegativeMultinomialDist extends DiscreteDistributionIntMulti {
protected double gamma;
protected double p[];
private static class Function implements MathFunction {
protected double Fl[];
protected int ups[];
protected int n;
protected int M;
protected int sumUps;
public Function (int n, int m, int ups[], double Fl[])
{
this.n = n;
this.M = m;
this.Fl = new double[Fl.length];
System.arraycopy (Fl, 0, this.Fl, 0, Fl.length);
this.ups = new int[ups.length];
System.arraycopy (ups, 0, this.ups, 0, ups.length);
sumUps = 0;
for (int i = 0; i < ups.length; i++)
sumUps += ups[i];
}
public double evaluate (double gamma)
{
double sum = 0.0;
for (int l = 0; l < M; l++)
sum += (Fl[l] / (gamma + (double) l));
return (sum - Math.log1p (sumUps / (n * gamma)));
}
}
/**
* Creates a NegativeMultinomialDist object with parameters γ =
* gamma and (p1,...,pd) such that
* ∑i=1dpi < 1,
* as described above. We have pi = p[i-1].
*
*/
public NegativeMultinomialDist (double gamma, double p[]) {
setParams (gamma, p);
}
public double prob (int x[]) {
return prob_ (gamma, p, x);
}
/*
public double cdf (int x[]) {
throw new UnsupportedOperationException ("cdf not implemented");
}
*/
public double[] getMean() {
return getMean_ (gamma, p);
}
public double[][] getCovariance() {
return getCovariance_ (gamma, p);
}
public double[][] getCorrelation() {
return getCorrelation_ (gamma, p);
}
private static void verifParam (double gamma, double p[]) {
double sumPi = 0.0;
if (gamma <= 0.0)
throw new IllegalArgumentException("gamma <= 0");
for (int i = 0; i < p.length;i++)
{
if ((p[i] < 0) || (p[i] >= 1))
throw new IllegalArgumentException("p is not a probability vector");
sumPi += p[i];
}
if (sumPi >= 1.0)
throw new IllegalArgumentException("p is not a probability vector");
}
private static double prob_ (double gamma, double p[], int x[]) {
double p0 = 0.0;
double sumPi= 0.0;
double sumXi= 0.0;
double sumLnXiFact = 0.0;
double sumXiLnPi = 0.0;
if (x.length != p.length)
throw new IllegalArgumentException ("x and p must have the same size");
for (int i = 0; i < p.length;i++)
{
sumPi += p[i];
sumXi += x[i];
sumLnXiFact += Num.lnFactorial (x[i]);
sumXiLnPi += x[i] * Math.log (p[i]);
}
p0 = 1.0 - sumPi;
return Math.exp (Num.lnGamma (gamma + sumXi) - (Num.lnGamma (gamma) + sumLnXiFact) +
gamma * Math.log (p0) + sumXiLnPi);
}
/**
* Computes the probability mass function
* of the negative multinomial distribution with parameters
* γ and (p1,...,pd), evaluated at x.
*
*/
public static double prob (double gamma, double p[], int x[]) {
verifParam (gamma, p);
return prob_ (gamma, p, x);
}
private static double cdf_ (double gamma, double p[], int x[]) {
throw new UnsupportedOperationException ("cdf not implemented");
}
/**
* Computes the cumulative probability function F of the
* negative multinomial distribution with parameters γ
* and (p1,...,pk), evaluated at x.
*
*/
public static double cdf (double gamma, double p[], int x[]) {
verifParam (gamma, p);
return cdf_ (gamma, p, x);
}
private static double[] getMean_ (double gamma, double p[]) {
double p0 = 0.0;
double sumPi= 0.0;
double mean[] = new double[p.length];
for (int i = 0; i < p.length;i++)
sumPi += p[i];
p0 = 1.0 - sumPi;
for (int i = 0; i < p.length; i++)
mean[i] = gamma * p[i] / p0;
return mean;
}
/**
* Computes the mean
* E[X] = γpi/p0 of the negative multinomial distribution
* with parameters γ and (p1,...,pd).
*
*/
public static double[] getMean (double gamma, double p[]) {
verifParam (gamma, p);
return getMean_ (gamma, p);
}
private static double[][] getCovariance_ (double gamma, double p[]) {
double p0 = 0.0;
double sumPi= 0.0;
double cov[][] = new double[p.length][p.length];
for (int i = 0; i < p.length;i++)
sumPi += p[i];
p0 = 1.0 - sumPi;
for (int i = 0; i < p.length; i++)
{
for (int j = 0; j < p.length; j++)
cov[i][j] = gamma * p[i] * p[j] / (p0 * p0);
cov[i][i] = gamma * p[i] * (p[i] + p0) / (p0 * p0);
}
return cov;
}
/**
* Computes the covariance matrix of the negative multinomial distribution
* with parameters γ and (p1,...,pd).
*
*/
public static double[][] getCovariance (double gamma, double p[]) {
verifParam (gamma, p);
return getCovariance_ (gamma, p);
}
private static double[][] getCorrelation_ (double gamma, double[] p) {
double corr[][] = new double[p.length][p.length];
double sumPi= 0.0;
double p0;
for (int i = 0; i < p.length;i++)
sumPi += p[i];
p0 = 1.0 - sumPi;
for (int i = 0; i < p.length; i++) {
for (int j = 0; j < p.length; j++)
corr[i][j] = Math.sqrt(p[i] * p[j] /((p0 + p[i]) * (p0 + p[j])));
corr[i][i] = 1.0;
}
return corr;
}
/**
* Computes the correlation matrix of the negative multinomial distribution
* with parameters γ and (p1,...,pd).
*
*/
public static double[][] getCorrelation (double gamma, double[] p) {
verifParam (gamma, p);
return getCorrelation_ (gamma, p);
}
/**
* It is now called getMLE.
*
*/
@Deprecated
public static double[] getMaximumLikelihoodEstimate (int x[][], int n,
int d) {
return getMLE(x, n, d);
}
/**
* Estimates and returns the parameters [
* hat(γ), hat(p_1),...,
* hat(p_d)]
* of the negative multinomial 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 vector
*
* @return returns the parameters [
* hat(γ), hat(p_1),...,hat(p_d)]
*
*/
public static double[] getMLE (int x[][], int n, int d) {
double parameters[] = new double[d + 1];
int ups[] = new int[n];
double mean[] = new double[d];
int i, j, l;
int M;
int prop;
// Initialization
for (i = 0; i < d; i++)
mean[i] = 0;
// Ups_j = Sum_k x_ji
// mean_i = Sum_n x_ji / n
for (j = 0; j < n; j++) {
ups[j] = 0;
for (i = 0; i < d; i++) {
ups[j] += x[j][i];
mean[i] += x[j][i];
}
}
for (i = 0; i < d; i++)
mean[i] /= n;
/*
double var = 0.0;
if (d > 1) {
// Calcule la covariance 0,1
for (j = 0; j < n; j++)
var += (x[j][0] - mean[0])*(x[j][1] - mean[1]);
var /= n;
} else {
// Calcule la variance 0
for (j = 0; j < n; j++)
var += (x[j][0] - mean[0])*(x[j][0] - mean[0]);
var /= n;
}
*/
// M = Max(Ups_j)
M = ups[0];
for (j = 1; j < n; j++)
if (ups[j] > M)
M = ups[j];
if (M >= Integer.MAX_VALUE)
throw new IllegalArgumentException("gamma/p_i too large");
double Fl[] = new double[M];
for (l = 0; l < M; l++) {
prop = 0;
for (j = 0; j < n; j++)
if (ups[j] > l)
prop++;
Fl[l] = (double) prop / (double) n;
}
/*
// Estime la valeur initiale de gamma pour brentDekker: pourrait
// accélérer brentDekker (gam0/1000, gam0*1000, f, 1e-5).
// Reste à bien tester.
if (d > 1) {
double gam0 = mean[0] * mean[1] / var;
System.out.println ("gam0 = " + gam0);
} else {
double t = var/mean[0] - 1.0;
double gam0 = mean[0] / t;
System.out.println ("gam0 = " + gam0);
}
*/
Function f = new Function (n, (int)M, ups, Fl);
parameters[0] = RootFinder.brentDekker (1e-9, 1e9, f, 1e-5);
double lambda[] = new double[d];
double sumLambda = 0.0;
for (i = 0; i < d; i++) {
lambda[i] = mean[i] / parameters[0];
sumLambda += lambda[i];
}
for (i = 0; i < d; i++) {
parameters[i + 1] = lambda[i] / (1.0 + sumLambda);
if (parameters[i + 1] > 1.0)
throw new IllegalArgumentException("p_i > 1");
}
return parameters;
}
/**
* Returns the parameter γ of this object.
*
*/
public double getGamma() {
return gamma;
}
/**
* Returns the parameters (p1,...,pd) of this object.
*
*/
public double[] getP() {
return p;
}
/**
* Sets the parameters γ and (p1,...,pd) of this object.
*
*/
public void setParams (double gamma, double p[]) {
double sumPi = 0.0;
if (gamma <= 0.0)
throw new IllegalArgumentException("gamma <= 0");
this.gamma = gamma;
this.dimension = p.length;
this.p = new double[dimension];
for (int i = 0; i < dimension; i++)
{
if ((p[i] < 0) || (p[i] >= 1))
throw new IllegalArgumentException("p is not a probability vector");
sumPi += p[i];
this.p[i] = p[i];
}
if (sumPi >= 1.0)
throw new IllegalArgumentException("p is not a probability vector");
}
}
Ads help maintain this website.