de.lmu.ifi.dbs.elki.math.statistics.distribution

Class GammaDistribution

java.lang.Object

de.lmu.ifi.dbs.elki.math.statistics.distribution.AbstractDistribution

de.lmu.ifi.dbs.elki.math.statistics.distribution.GammaDistribution

All Implemented Interfaces:

Distribution

Direct Known Subclasses:

ChiSquaredDistribution

public class GammaDistribution
extends AbstractDistribution

Gamma Distribution, with random generation and density functions.

Author:

Erich Schubert

Nested Class Summary

Nested Classes

Modifier and Type	Class and Description
static class	GammaDistribution.Parameterizer Parameterization class TODO: allow alternate parameterization, with alpha+beta?

Field Summary

Fields

Modifier and Type	Field and Description
static double	EULERS_CONST Euler–Mascheroni constant
private double	k Alpha == k
(package private) static double[]	LANCZOS LANCZOS-Coefficients for Gamma approximation.
(package private) static int	MAX_ITERATIONS Maximum number of iterations for regularizedGammaP.
(package private) static double	NUM_PRECISION Numerical precision to use (data type dependent!)
private double	theta Theta == 1 / Beta

Fields inherited from class de.lmu.ifi.dbs.elki.math.statistics.distribution.AbstractDistribution

random

Constructor Summary

Constructors

Constructor and Description
GammaDistribution(double k, double theta) Constructor for Gamma distribution.
GammaDistribution(double k, double theta, Random random) Constructor for Gamma distribution.
GammaDistribution(double k, double theta, RandomFactory random) Constructor for Gamma distribution.

Method Summary

Methods

Modifier and Type	Method and Description
double	cdf(double val) Return the cumulative density function at the given value.
static double	cdf(double val, double k, double theta) The CDF, static version.
protected static double	chisquaredProbitApproximation(double p, double nu, double g) Approximate probit for chi squared distribution Based on first half of algorithm AS 91 Reference: Algorithm AS 91: The percentage points of the $\chi^2$ distribution D.J.
static double	digamma(double x) Compute the Psi / Digamma function Reference: J.
static double	gamma(double x) Compute the regular Gamma function.
protected static double	gammaQuantileNewtonRefinement(double logpt, double k, double theta, int maxit, double x) Refinement of ChiSquared probit using Newton iterations.
double	getK()
double	getTheta()
static double	logcdf(double val, double k, double theta) The log CDF, static version.
static double	logGamma(double x) Compute logGamma.
static double	logpdf(double x, double k, double theta) Gamma distribution PDF (with 0.0 for x < 0)
static double	logregularizedGammaP(double a, double x) Returns the regularized gamma function log P(a, x).
double	nextRandom() Generate a new random value
static double	nextRandom(double k, double theta, Random random) Generate a random value with the generators parameters.
double	pdf(double val) Return the density of an existing value
static double	pdf(double x, double k, double theta) Gamma distribution PDF (with 0.0 for x < 0)
double	quantile(double val) Quantile aka probit (for normal) aka inverse CDF (invcdf, cdf^-1) function.
static double	quantile(double p, double k, double theta) Compute probit (inverse cdf) for Gamma distributions.
static double	regularizedGammaP(double a, double x) Returns the regularized gamma function P(a, x).
static double	regularizedGammaQ(double a, double x) Returns the regularized gamma function Q(a, x) = 1 - P(a, x).
String	toString() Simple toString explaining the distribution parameters.
static double	trigamma(double x) Compute the Trigamma function.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, wait, wait, wait

Field Detail

EULERS_CONST

public static final double EULERS_CONST

Euler–Mascheroni constant

See Also:

Constant Field Values

LANCZOS

static final double[] LANCZOS

LANCZOS-Coefficients for Gamma approximation. These are said to have higher precision than those in "Numerical Recipes". They probably come from Paul Godfrey: http://my.fit.edu/~gabdo/gamma.txt

NUM_PRECISION

static final double NUM_PRECISION

Numerical precision to use (data type dependent!) If you change this, make sure to test exhaustively!

See Also:

Constant Field Values

MAX_ITERATIONS

static final int MAX_ITERATIONS

Maximum number of iterations for regularizedGammaP. To prevent degeneration for extreme values. FIXME: is this too high, too low? Can we improve behavior for extreme cases?

See Also:

Constant Field Values

k

private final double k

Alpha == k

theta

private final double theta

Theta == 1 / Beta

Constructor Detail

GammaDistribution

public GammaDistribution(double k,
                 double theta,
                 Random random)

Constructor for Gamma distribution.

Parameters:

k - k, alpha aka. "shape" parameter

theta - Theta = 1.0/Beta aka. "scaling" parameter

random - Random generator

GammaDistribution

public GammaDistribution(double k,
                 double theta,
                 RandomFactory random)

Constructor for Gamma distribution.

Parameters:

k - k, alpha aka. "shape" parameter

theta - Theta = 1.0/Beta aka. "scaling" parameter

random - Random generator

GammaDistribution

public GammaDistribution(double k,
                 double theta)

Constructor for Gamma distribution.

Parameters:

k - k, alpha aka. "shape" parameter

theta - Theta = 1.0/Beta aka. "scaling" parameter

Method Detail

pdf

public double pdf(double val)

Description copied from interface: Distribution

Return the density of an existing value

Parameters:

val - existing value

Returns:

distribution density

cdf

public double cdf(double val)

Description copied from interface: Distribution

Return the cumulative density function at the given value.

Parameters:

val - existing value

Returns:

cumulative density

quantile

public double quantile(double val)

Description copied from interface: Distribution

Quantile aka probit (for normal) aka inverse CDF (invcdf, cdf^-1) function.

Parameters:

val - Quantile to find

Returns:

Quantile position

nextRandom

public double nextRandom()

Description copied from interface: Distribution

Generate a new random value

Specified by:

nextRandom in interface Distribution

Overrides:

nextRandom in class AbstractDistribution

Returns:

new random value

toString

public String toString()

Simple toString explaining the distribution parameters. Used in producing a model description.

Specified by:

toString in interface Distribution

Overrides:

toString in class Object

Returns:

description

getK

public double getK()

Returns:

the value of k

getTheta

public double getTheta()

Returns:

the standard deviation

cdf

public static double cdf(double val,
         double k,
         double theta)

The CDF, static version.

Parameters:

val - Value

k - Shape k

theta - Theta = 1.0/Beta aka. "scaling" parameter

Returns:

cdf value

logcdf

public static double logcdf(double val,
            double k,
            double theta)

The log CDF, static version.

Parameters:

val - Value

k - Shape k

theta - Theta = 1.0/Beta aka. "scaling" parameter

Returns:

cdf value

pdf

public static double pdf(double x,
         double k,
         double theta)

Gamma distribution PDF (with 0.0 for x < 0)

Parameters:

x - query value

k - Alpha

theta - Theta = 1 / Beta

Returns:

probability density

logpdf

public static double logpdf(double x,
            double k,
            double theta)

Gamma distribution PDF (with 0.0 for x < 0)

Parameters:

x - query value

k - Alpha

theta - Theta = 1 / Beta

Returns:

probability density

logGamma

public static double logGamma(double x)

Compute logGamma. Based loosely on "Numerical Recpies" and the work of Paul Godfrey at http://my.fit.edu/~gabdo/gamma.txt TODO: find out which approximation really is the best...

Parameters:

x - Parameter x

Returns:

log(Γ(x))

gamma

public static double gamma(double x)

Compute the regular Gamma function. Note: for numerical reasons, it is preferrable to use logGamma(double) when possible! In particular, this method just computes Math.exp(logGamma(x)) anyway. Try to postpone the Math.exp call to preserve numeric range!

Parameters:

x - Position

Returns:

Gamma at this position

regularizedGammaP

public static double regularizedGammaP(double a,
                       double x)

Returns the regularized gamma function P(a, x). Includes the quadrature way of computing. TODO: find "the" most accurate version of this. We seem to agree with others for the first 10+ digits, but diverge a bit later than that.

Parameters:

a - Parameter a

x - Parameter x

Returns:

Gamma value

logregularizedGammaP

public static double logregularizedGammaP(double a,
                          double x)

Returns the regularized gamma function log P(a, x). Includes the quadrature way of computing. TODO: find "the" most accurate version of this. We seem to agree with others for the first 10+ digits, but diverge a bit later than that.

Parameters:

a - Parameter a

x - Parameter x

Returns:

Gamma value

regularizedGammaQ

public static double regularizedGammaQ(double a,
                       double x)

Returns the regularized gamma function Q(a, x) = 1 - P(a, x). Includes the continued fraction way of computing, based loosely on the book "Numerical Recipes"; but probably not with the exactly same precision, since we reimplemented this in our coding style, not literally. TODO: find "the" most accurate version of this. We seem to agree with others for the first 10+ digits, but diverge a bit later than that.

Parameters:

a - parameter a

x - parameter x

Returns:

Result

nextRandom

public static double nextRandom(double k,
                double theta,
                Random random)

Generate a random value with the generators parameters. Along the lines of - J.H. Ahrens, U. Dieter (1974): Computer methods for sampling from gamma, beta, Poisson and binomial distributions, Computing 12, 223-246. - J.H. Ahrens, U. Dieter (1982): Generating gamma variates by a modified rejection technique, Communications of the ACM 25, 47-54.

Parameters:

k - K parameter

theta - Theta parameter

random - Random generator

chisquaredProbitApproximation

@Reference(title="Algorithm AS 91: The percentage points of the $\\chi^2$ distribution",
           authors="D.J. Best, D. E. Roberts",
           booktitle="Journal of the Royal Statistical Society. Series C (Applied Statistics)")
protected static double chisquaredProbitApproximation(double p,
                                                                  double nu,
                                                                  double g)

Approximate probit for chi squared distribution Based on first half of algorithm AS 91 Reference:

Algorithm AS 91: The percentage points of the $\chi^2$ distribution
D.J. Best, D. E. Roberts
Journal of the Royal Statistical Society. Series C (Applied Statistics)

Parameters:

p - Probit value

nu - Shape parameter for Chi, nu = 2 * k

g - log(nu)

Returns:

Probit for chi squared

quantile

@Reference(title="Algorithm AS 91: The percentage points of the $\\chi^2$ distribution",
           authors="D.J. Best, D. E. Roberts",
           booktitle="Journal of the Royal Statistical Society. Series C (Applied Statistics)")
public static double quantile(double p,
                                             double k,
                                             double theta)

Compute probit (inverse cdf) for Gamma distributions. Based on algorithm AS 91: Reference:

Algorithm AS 91: The percentage points of the $\chi^2$ distribution
D.J. Best, D. E. Roberts
Journal of the Royal Statistical Society. Series C (Applied Statistics)

Parameters:

p - Probability

k - k, alpha aka. "shape" parameter

theta - Theta = 1.0/Beta aka. "scaling" parameter

Returns:

Probit for Gamma distribution

gammaQuantileNewtonRefinement

protected static double gammaQuantileNewtonRefinement(double logpt,
                                   double k,
                                   double theta,
                                   int maxit,
                                   double x)

Refinement of ChiSquared probit using Newton iterations. A trick used by GNU R to improve precision.

Parameters:

logpt - Target value of log p

k - Alpha

theta - Theta = 1 / Beta

maxit - Maximum number of iterations to do

x - Initial estimate

Returns:

Refined value

digamma

@Reference(authors="J. M. Bernando",
           title="Algorithm AS 103: Psi (Digamma) Function",
           booktitle="Statistical Algorithms")
public static double digamma(double x)

Compute the Psi / Digamma function Reference:

J. M. Bernando
Algorithm AS 103: Psi (Digamma) Function
Statistical Algorithms

TODO: is there a more accurate version maybe in R?

Parameters:

x - Position

Returns:

digamma value

trigamma

public static double trigamma(double x)

Compute the Trigamma function. Based on digamma. TODO: is there a more accurate version maybe in R?

Parameters:

x - Position

Returns:

trigamma value