Example usage for org.apache.commons.math3.special Gamma LANCZOS_G

Introduction

In this page you can find the example usage for org.apache.commons.math3.special Gamma LANCZOS_G.

Prototype

double LANCZOS_G

To view the source code for org.apache.commons.math3.special Gamma LANCZOS_G.

Click Source Link

Document

The value of the g constant in the Lanczos approximation, see #lanczos(double) .

Usage

From source file:statalign.utils.GammaDistribution.java

/**
 * Creates a Gamma distribution.//from   ww  w.  j  ava 2s  . com
 *
 * @param rng Random number generator.
 * @param shape the shape parameter
 * @param scale the scale parameter
 * @param inverseCumAccuracy the maximum absolute error in inverse
 * cumulative probability estimates (defaults to
 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 * @since 3.1
 */
public GammaDistribution(RandomGenerator rng, double shape, double scale, double inverseCumAccuracy)
        throws NotStrictlyPositiveException {
    super(rng);

    if (shape <= 0) {
        throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE, shape);
    }
    if (scale <= 0) {
        throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale);
    }

    this.shape = shape;
    this.scale = scale;
    this.solverAbsoluteAccuracy = inverseCumAccuracy;
    this.shiftedShape = shape + Gamma.LANCZOS_G + 0.5;
    final double aux = FastMath.E / (2.0 * FastMath.PI * shiftedShape);
    this.densityPrefactor2 = shape * FastMath.sqrt(aux) / Gamma.lanczos(shape);
    this.densityPrefactor1 = this.densityPrefactor2 / scale * FastMath.pow(shiftedShape, -shape)
            * FastMath.exp(shape + Gamma.LANCZOS_G);
    this.minY = shape + Gamma.LANCZOS_G - FastMath.log(Double.MAX_VALUE);
    this.maxLogY = FastMath.log(Double.MAX_VALUE) / (shape - 1.0);
}

---

From source file:statalign.utils.GammaDistribution.java

/** {@inheritDoc} */
@Override// w  w  w . java 2 s.co m
public double density(double x) {
    /* The present method must return the value of
     *
     *     1       x a     - x
     * ---------- (-)  exp(---)
     * x Gamma(a)  b        b
     *
     * where a is the shape parameter, and b the scale parameter.
     * Substituting the Lanczos approximation of Gamma(a) leads to the
     * following expression of the density
     *
     * a              e            1         y      a
     * - sqrt(------------------) ---- (-----------)  exp(a - y + g),
     * x      2 pi (a + g + 0.5)  L(a)  a + g + 0.5
     *
     * where y = x / b. The above formula is the "natural" computation, which
     * is implemented when no overflow is likely to occur. If overflow occurs
     * with the natural computation, the following identity is used. It is
     * based on the BOOST library
     * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/igamma.html
     * Formula (15) needs adaptations, which are detailed below.
     *
     *       y      a
     * (-----------)  exp(a - y + g)
     *  a + g + 0.5
     *                              y - a - g - 0.5    y (g + 0.5)
     *               = exp(a log1pm(---------------) - ----------- + g),
     *                                a + g + 0.5      a + g + 0.5
     *
     *  where log1pm(z) = log(1 + z) - z. Therefore, the value to be
     *  returned is
     *
     * a              e            1
     * - sqrt(------------------) ----
     * x      2 pi (a + g + 0.5)  L(a)
     *                              y - a - g - 0.5    y (g + 0.5)
     *               * exp(a log1pm(---------------) - ----------- + g).
     *                                a + g + 0.5      a + g + 0.5
     */
    if (x < 0) {
        return 0;
    }
    final double y = x / scale;
    if ((y <= minY) || (FastMath.log(y) >= maxLogY)) {
        /*
         * Overflow.
         */
        final double aux1 = (y - shiftedShape) / shiftedShape;
        final double aux2 = shape * (FastMath.log1p(aux1) - aux1);
        final double aux3 = -y * (Gamma.LANCZOS_G + 0.5) / shiftedShape + Gamma.LANCZOS_G + aux2;
        return densityPrefactor2 / x * FastMath.exp(aux3);
    }
    /*
     * Natural calculation.
     */
    return densityPrefactor1 * FastMath.exp(-y) * FastMath.pow(y, shape - 1);
}

---