Package elki.math.statistics.distribution

Class GammaDistribution

  • All Implemented Interfaces:
    Distribution
    Direct Known Subclasses:
    ChiSquaredDistribution
    public class GammaDistribution
extends java.lang.Object
implements Distribution
    Gamma Distribution, with random generation and density functions.
    Since:
    0.4.0
    Author:
    Erich Schubert

    • Field Summary

      Fields 
      Modifier and Type Field Description
      static double EULERS_CONST
      Euler–Mascheroni constant
      (package private) static double FPMIN
      Precision threshold in regularizedGammaQ(double, double)
      private double k
      Alpha == k
      (package private) static double[] LANCZOS
      LANCZOS-Coefficients for Gamma approximation.
      (package private) static double LOGGAMMA_G
      Constant used in the logGamma function.
      (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

    • Constructor Summary

      Constructors 
      Constructor Description
      GammaDistribution​(double k, double theta)
      Constructor for Gamma distribution.

    • Method Summary

      All Methods Static Methods Instance Methods Concrete Methods 
      Modifier and Type Method 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
      static double digamma​(double x)
      Compute the Psi / Digamma function
      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.
      double logpdf​(double val)
      Return the log density of an existing value
      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).
      static double nextRandom​(double k, double theta, java.util.Random random)
      Generate a random value with the generators parameters.
      double nextRandom​(java.util.Random random)
      Generate a new random value
      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).
      java.lang.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

      • LOGGAMMA_G

        static final double LOGGAMMA_G
        Constant used in the logGamma function.
        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)
        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
        Specified by:
        pdf in interface Distribution
        Parameters:
        val - existing value
        Returns:
        distribution density

      • logpdf

        public double logpdf​(double val)
        Description copied from interface: Distribution
        Return the log density of an existing value
        Specified by:
        logpdf in interface Distribution
        Parameters:
        val - existing value
        Returns:
        log distribution density

      • cdf

        public double cdf​(double val)
        Description copied from interface: Distribution
        Return the cumulative density function at the given value.
        Specified by:
        cdf in interface Distribution
        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.
        Specified by:
        quantile in interface Distribution
        Parameters:
        val - Quantile to find
        Returns:
        Quantile position

      • nextRandom

        public double nextRandom​(java.util.Random random)
        Description copied from interface: Distribution
        Generate a new random value
        Specified by:
        nextRandom in interface Distribution
        Parameters:
        random - Random number generator
        Returns:
        new random value

      • toString

        public java.lang.String toString()
        Simple toString explaining the distribution parameters. Used in producing a model description.
        Specified by:
        toString in interface Distribution
        Overrides:
        toString in class java.lang.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 Recipes" 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 preferable to use logGamma(double) when possible! In particular, this method just computes FastMath.exp(logGamma(x)) anyway.

        Try to postpone the FastMath.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

        @Reference(authors="J. H. Ahrens, U. Dieter",title="Computer methods for sampling from gamma, beta, Poisson and binomial distributions",booktitle="Computing 12",url="https://doi.org/10.1007/BF02293108",bibkey="DBLP:journals/computing/AhrensD74") @Reference(authors="J. H. Ahrens, U. Dieter",title="Generating gamma variates by a modified rejection technique",booktitle="Communications of the ACM 25",url="https://doi.org/10.1145/358315.358390",bibkey="DBLP:journals/cacm/AhrensD82")
public static double nextRandom​(double k,
                                double theta,
                                java.util.Random random)
        Generate a random value with the generators parameters.

        Along the lines of

        J. H. Ahrens, U. Dieter
        Computer methods for sampling from gamma, beta, Poisson and binomial distributions
        Computing 12

        J. H. Ahrens, U. Dieter
        Generating gamma variates by a modified rejection technique
        Communications of the ACM 25

        Parameters:
        k - K parameter
        theta - Theta parameter
        random - Random generator

      • chisquaredProbitApproximation

        @Reference(authors="D. J. Best, D. E. Roberts",
           title="Algorithm AS 91: The percentage points of the \u03c7\u00b2 distribution",
           booktitle="Journal of the Royal Statistical Society. Series C (Applied Statistics)",
           url="https://doi.org/10.2307/2347113",
           bibkey="doi:10.2307/2347113")
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:

        D. J. Best, D. E. Roberts
        Algorithm AS 91: The percentage points of the χ² distribution
        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(authors="D. J. Best, D. E. Roberts",
           title="Algorithm AS 91: The percentage points of the \u03c7\u00b2 distribution",
           booktitle="Journal of the Royal Statistical Society. Series C (Applied Statistics)",
           url="https://doi.org/10.2307/2347113",
           bibkey="doi:10.2307/2347113")
public static double quantile​(double p,
                              double k,
                              double theta)
        Compute probit (inverse cdf) for Gamma distributions.

        Based on algorithm AS 91:

        Reference:

        D. J. Best, D. E. Roberts
        Algorithm AS 91: The percentage points of the χ² distribution
        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",
           url="https://doi.org/10.2307/2347257",
           bibkey="doi:10.2307/2347257")
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