umontreal.iro.lecuyer.probdist.GammaDist 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: GammaDist
* Description: gamma 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.probdist;
import umontreal.iro.lecuyer.util.Num;
import umontreal.iro.lecuyer.functions.MathFunction;
import umontreal.iro.lecuyer.util.RootFinder;
/**
* Extends the class {@link ContinuousDistribution} for
* the gamma distribution with
* shape parameter
* α > 0 and scale parameter
*
* λ > 0.
* The density is
*
*
*
* f (x) = λαxα-1e-λx/Γ(α), for x > 0,
*
* where Γ is the gamma function, defined by
*
*
*
* Γ(α) = ∫0∞xα-1e-xdx.
*
* In particular,
* Γ(n) = (n - 1)! when n is a positive integer.
*
*/
public class GammaDist extends ContinuousDistribution {
private double alpha;
private double lambda;
private double logFactor; // Log (lambda^alpha / Gamma (alpha))
private static final double ALIM = 1.0E5;
private static class Function implements MathFunction {
// For MLE
private int n;
private double empiricalMean;
private double sumLn;
public Function (int n, double empiricalMean, double sumLn) {
this.n = n;
this.empiricalMean = empiricalMean;
this.sumLn = sumLn;
}
public double evaluate (double x) {
if (x <= 0.0) return 1.0e200;
return (n * Math.log (empiricalMean / x) + n * Num.digamma (x) - sumLn);
}
}
private static class myFunc implements MathFunction {
// For inverseF
protected int d;
protected double alp, u;
public myFunc (double alp, int d, double u) {
this.alp = alp;
this.d = d;
this.u = u;
}
public double evaluate (double x) {
return u - GammaDist.cdf(alp, d, x);
}
}
/**
* Constructs a GammaDist object with parameters
* α = alpha and λ = 1.
*
*/
public GammaDist (double alpha) {
setParams (alpha, 1.0, decPrec);
}
/**
* Constructs a GammaDist object with parameters
* α = alpha and λ = lambda.
*
*/
public GammaDist (double alpha, double lambda) {
setParams (alpha, lambda, decPrec);
}
/**
* Constructs a GammaDist object with parameters α =
* alpha and λ = lambda, and approximations of
* roughly d decimal digits of precision when computing functions.
*
*/
public GammaDist (double alpha, double lambda, int d) {
setParams (alpha, lambda, d);
}
private static double mybelog (double x)
/*
* This is the function 1 + (1 - x*x + 2*x*log(x)) / ((1 - x)*(1 - x))
*/
{
if (x < 1.0e-30)
return 0.0;
if (x > 1.0e30)
return 2.0*(Math.log(x) - 1.0) / x;
if (x == 1.0)
return 1.0;
double t = 1.0 - x;
if (x < 0.9 || x > 1.1) {
double w = (t + x*Math.log(x)) / (t*t);
return 2.0 * w;
}
// For x near 1, use a series expansion to avoid loss of precision.
double term;
final double EPS = 1.0e-12;
double tpow = 1.0;
double sum = 0.5;
int j = 3;
do {
tpow *= t;
term = tpow / (j * (j - 1));
sum += term;
j++;
} while (Math.abs (term / sum) > EPS);
return 2.0*sum;
}
public double density (double x) {
if (x <= 0)
return 0.0;
double z = logFactor + (alpha - 1.0) * Math.log(x) - lambda * x;
if (z > -XBIGM)
return Math.exp (z);
else
return 0.0;
}
public double cdf (double x) {
return cdf (alpha, lambda, decPrec, x);
}
public double barF (double x) {
return barF (alpha, lambda, decPrec, x);
}
public double inverseF (double u) {
return inverseF (alpha, decPrec, u)/lambda;
}
public double getMean() {
return GammaDist.getMean (alpha, lambda);
}
public double getVariance() {
return GammaDist.getVariance (alpha, lambda);
}
public double getStandardDeviation() {
return GammaDist.getStandardDeviation (alpha, lambda);
}
/**
* Computes the density function at x.
*
*/
public static double density (double alpha, double lambda, double x) {
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0)
throw new IllegalArgumentException ("lambda <= 0");
if (x <= 0)
return 0.0;
double z = alpha * Math.log (lambda*x) - lambda*x - Num.lnGamma (alpha);
if (z > -XBIGM)
return Math.exp (z) / x;
else
return 0.0;
}
/**
* Returns an approximation of the gamma distribution
* function with parameters α = alpha and
* λ = lambda.
* The function tries to return d decimals digits of precision.
* For α not too large (e.g.,
* α <= 1000),
* d gives a good idea of the precision attained.
*
*/
public static double cdf (double alpha, double lambda, int d, double x) {
return cdf (alpha, d, lambda*x);
}
/**
* Equivalent to cdf (alpha, 1.0, d, x).
*
*/
public static double cdf (double alpha, int d, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (d <= 0)
throw new IllegalArgumentException ("d <= 0");
if (x <= 0.0)
return 0.0;
if (1.0 == alpha)
return ExponentialDist.cdf (1.0, x);
if (alpha > 10.0) {
if (x > alpha * 10.0)
return 1.0;
} else {
if (x > XBIG)
return 1.0;
}
if (alpha >= ALIM) {
double d2 = x + 1.0/3.0 - alpha - 0.02/alpha;
double S = alpha - 1.0/2.0;
double w = mybelog(S/x);
double y = d2 * Math.sqrt(w/x);
return NormalDist.cdf01 (y);
}
if (x <= 1.0 || x < alpha) {
double factor, z, rn, term;
factor = Math.exp (alpha*Math.log (x) - x - Num.lnGamma (alpha));
final double EPS = EPSARRAY[d];
z = 1.0;
term = 1.0;
rn = alpha;
do {
rn += 1.0;
term *= x/rn;
z += term;
} while (term >= EPS * z);
return z*factor/alpha;
} else
return 1.0 - barF (alpha, d, x);
}
/**
* Computes the complementary distribution function.
*
*/
public static double barF (double alpha, double lambda, int d, double x) {
return barF (alpha, d, lambda*x);
}
/**
* Same as {@link #barF(double,double,int,double) barF} (alpha, 1.0, d, x).
*
*/
public static double barF (double alpha, int d, double x) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (d <= 0)
throw new IllegalArgumentException ("d <= 0");
if (x <= 0.0)
return 1.0;
if (1.0 == alpha)
return ExponentialDist.barF (1.0, x);
if (alpha >= 70.0) {
if (x >= alpha * XBIG)
return 0.0;
} else {
if (x >= XBIGM)
return 0.0;
}
if (alpha >= ALIM) {
double d2 = x + 1.0/3.0 - alpha - 0.02/alpha;
double S = alpha - 1.0/2.0;
double w = mybelog(S/x);
double y = d2 * Math.sqrt(w/x);
return NormalDist.barF01 (y);
}
if (x <= 1.0 || x < alpha)
return 1.0 - cdf (alpha, d, x);
double[] V = new double[6];
final double EPS = EPSARRAY[d];
final double RENORM = 1.0E100;
double R, dif;
int i;
double factor = Math.exp (alpha*Math.log (x) - x - Num.lnGamma (alpha));
double A = 1.0 - alpha;
double B = A + x + 1.0;
double term = 0.0;
V[0] = 1.0;
V[1] = x;
V[2] = x + 1.0;
V[3] = x * B;
double res = V[2]/V[3];
do {
A += 1.0;
B += 2.0;
term += 1.0;
V[4] = B * V[2] - A * term * V[0];
V[5] = B * V[3] - A * term * V[1];
if (V[5] != 0.0) {
R = V[4]/V[5];
dif = Math.abs (res - R);
if (dif <= EPS*R)
return factor*res;
res = R;
}
for (i = 0; i < 4; i++)
V[i] = V[i + 2];
if (Math.abs (V[4]) >= RENORM) {
for (i = 0; i < 4; i++)
V[i] /= RENORM;
}
} while (true);
}
/**
* Computes the inverse distribution function.
*
*/
public static double inverseF (double alpha, double lambda, int d,
double u) {
return inverseF (alpha, d, u)/lambda;
}
/**
* Same as
* {@link #inverseF(double,double,int,double) inverseF} (alpha, 1, d, u).
*
*/
public static double inverseF (double alpha, int d, double u) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (u > 1.0 || u < 0.0)
throw new IllegalArgumentException ("u not in [0,1]");
if (u <= 0.0)
return 0;
if (u >= 1.0)
return Double.POSITIVE_INFINITY;
if (d <= 0)
throw new IllegalArgumentException ("d <= 0");
if (d > 15)
d = 15;
final double EPS = Math.pow (10.0, -d);
double sigma = GammaDist.getStandardDeviation (alpha, 1.0);
double x = NormalDist.inverseF (alpha, sigma, u);
double v = GammaDist.cdf (alpha, d, x);
double xmax;
if (alpha < 1.0)
xmax = 100.0;
else
xmax = alpha + 40.0 * sigma;
myFunc f = new myFunc (alpha, d, u);
if (u <= 1.0e-8 || alpha <= 1.5) {
if (v < u)
return RootFinder.bisection (x, xmax, f, EPS);
else
return RootFinder.bisection (0, x, f, EPS);
} else {
if (v < u)
return RootFinder.brentDekker (x, xmax, f, EPS);
else
return RootFinder.brentDekker (0, x, f, EPS);
}
}
/**
* Estimates the parameters
* (α, λ) of the gamma distribution
* using the maximum likelihood method, from the n observations
* x[i],
* i = 0, 1,…, n - 1. The estimates are returned in a two-element
* array, in regular order: [α, λ].
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
* @return returns the parameters [
* hat(α),
* hat(λ)]
*
*/
public static double[] getMLE (double[] x, int n) {
double parameters[];
double sum = 0.0;
double sumLn = 0.0;
double empiricalMean;
double alphaMME;
double a;
final double LN_EPS = Num.LN_DBL_MIN - Num.LN2;
parameters = new double[2];
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
for (int i = 0; i < n; i++)
{
sum += x[i];
if (x[i] <= 0.0)
sumLn += LN_EPS;
else
sumLn += Math.log (x[i]);
}
empiricalMean = sum / (double) n;
sum = 0.0;
for (int i = 0; i < n; i++) {
sum += (x[i] - empiricalMean) * (x[i] - empiricalMean);
}
alphaMME = (empiricalMean * empiricalMean * (double) n) / sum;
if ((a = alphaMME - 10.0) <= 0) {
a = 1.0e-5;
}
Function f = new Function (n, empiricalMean, sumLn);
parameters[0] = RootFinder.brentDekker (a, alphaMME + 10.0, f, 1e-7);
parameters[1] = parameters[0] / empiricalMean;
return parameters;
}
/**
* Creates a new instance of a gamma distribution with parameters α and λ
* estimated using the maximum likelihood method based on the n observations
* x[i],
* i = 0, 1,…, n - 1.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
*
*/
public static GammaDist getInstanceFromMLE (double[] x, int n) {
double parameters[] = getMLE (x, n);
return new GammaDist (parameters[0], parameters[1]);
}
/**
* Computes and returns the mean
* E[X] = α/λ
* of the gamma distribution with parameters α and λ.
*
* @return the mean of the gamma distribution
* E[X] = α/λ
*
*/
public static double getMean (double alpha, double lambda) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
return (alpha / lambda);
}
/**
* Computes and returns the variance
* Var[X] = α/λ2
* of the gamma distribution with parameters α and λ.
*
* @return the variance of the gamma distribution
* Var[X] = α/λ2
*
*/
public static double getVariance (double alpha, double lambda) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
return (alpha / (lambda * lambda));
}
/**
* Computes and returns the standard deviation of the gamma
* distribution with parameters α and λ.
*
* @return the standard deviation of the gamma distribution
*
*/
public static double getStandardDeviation (double alpha, double lambda) {
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0.0)
throw new IllegalArgumentException ("lambda <= 0");
return (Math.sqrt(alpha) / lambda);
}
/**
* Return the parameter α for this object.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Return the parameter λ for this object.
*
*/
public double getLambda() {
return lambda;
}
public void setParams (double alpha, double lambda, int d) {
if (alpha <= 0)
throw new IllegalArgumentException ("alpha <= 0");
if (lambda <= 0)
throw new IllegalArgumentException ("lambda <= 0");
this.alpha = alpha;
this.lambda = lambda;
this.decPrec = d;
logFactor = alpha * Math.log(lambda) - Num.lnGamma (alpha);
supportA = 0.0;
}
/**
* Return a table containing the parameters of the current distribution.
* This table is put in regular order: [α, λ].
*
*
*/
public double[] getParams () {
double[] retour = {alpha, lambda};
return retour;
}
public String toString () {
return getClass().getSimpleName() + " : alpha = " + alpha + ", lambda = " + lambda;
}
}