Java Code Examples for org.apache.commons.math3.special.Gamma#LANCZOS_G

The following examples show how to use org.apache.commons.math3.special.Gamma#LANCZOS_G . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** {@inheritDoc} **/
@Override
public double logDensity(double x) {
    /*
     * see the comment in {@link #density(double)} for computation details
     */
    if (x < 0) {
        return Double.NEGATIVE_INFINITY;
    }
    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 logDensityPrefactor2 - FastMath.log(x) + aux3;
    }
    /*
     * Natural calculation.
     */
    return logDensityPrefactor1 - y + FastMath.log(y) * (shape - 1);
}
 
Example 2
Source File: GammaDistributionTest.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
public static double logGamma(double x) {
    /*
     * This is a copy of
     * double Gamma.logGamma(double)
     * prior to MATH-849
     */
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double sum = Gamma.lanczos(x);
        double tmp = x + Gamma.LANCZOS_G + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 
Example 3
Source File: GammaDistributionTest.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
public static double logGamma(double x) {
    /*
     * This is a copy of
     * double Gamma.logGamma(double)
     * prior to MATH-849
     */
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double sum = Gamma.lanczos(x);
        double tmp = x + Gamma.LANCZOS_G + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 
Example 4
Source File: GammaDistributionTest.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
public static double logGamma(double x) {
    /*
     * This is a copy of
     * double Gamma.logGamma(double)
     * prior to MATH-849
     */
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double sum = Gamma.lanczos(x);
        double tmp = x + Gamma.LANCZOS_G + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 
Example 5
Source File: GammaDistributionTest.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
public static double logGamma(double x) {
    /*
     * This is a copy of
     * double Gamma.logGamma(double)
     * prior to MATH-849
     */
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double sum = Gamma.lanczos(x);
        double tmp = x + Gamma.LANCZOS_G + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 
Example 6
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
/** {@inheritDoc} **/
@Override
public double logDensity(double x) {
    /*
     * see the comment in {@link #density(double)} for computation details
     */
    if (x < 0) {
        return Double.NEGATIVE_INFINITY;
    }
    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 logDensityPrefactor2 - FastMath.log(x) + aux3;
    }
    /*
     * Natural calculation.
     */
    return logDensityPrefactor1 - y + FastMath.log(y) * (shape - 1);
}
 
Example 7
Source File: GammaDistributionTest.java    From astor with GNU General Public License v2.0 6 votes vote down vote up
public static double logGamma(double x) {
    /*
     * This is a copy of
     * double Gamma.logGamma(double)
     * prior to MATH-849
     */
    double ret;

    if (Double.isNaN(x) || (x <= 0.0)) {
        ret = Double.NaN;
    } else {
        double sum = Gamma.lanczos(x);
        double tmp = x + Gamma.LANCZOS_G + .5;
        ret = ((x + .5) * FastMath.log(tmp)) - tmp +
            HALF_LOG_2_PI + FastMath.log(sum / x);
    }

    return ret;
}
 
Example 8
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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);
}
 
Example 9
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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.logDensityPrefactor2 = FastMath.log(shape) + 0.5 * FastMath.log(aux) -
                                FastMath.log(Gamma.lanczos(shape));
    this.densityPrefactor1 = this.densityPrefactor2 / scale *
            FastMath.pow(shiftedShape, -shape) *
            FastMath.exp(shape + Gamma.LANCZOS_G);
    this.logDensityPrefactor1 = this.logDensityPrefactor2 - FastMath.log(scale) -
            FastMath.log(shiftedShape) * shape +
            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);
}
 
Example 10
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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);
}
 
Example 11
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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.logDensityPrefactor2 = FastMath.log(shape) + 0.5 * FastMath.log(aux) -
                                FastMath.log(Gamma.lanczos(shape));
    this.densityPrefactor1 = this.densityPrefactor2 / scale *
            FastMath.pow(shiftedShape, -shape) *
            FastMath.exp(shape + Gamma.LANCZOS_G);
    this.logDensityPrefactor1 = this.logDensityPrefactor2 - FastMath.log(scale) -
            FastMath.log(shiftedShape) * shape +
            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);
}
 
Example 12
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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);
}
 
Example 13
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 5 votes vote down vote up
/**
 * Creates a Gamma distribution.
 *
 * @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);
}
 
Example 14
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 15
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 16
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 17
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 18
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 19
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}
 
Example 20
Source File: GammaDistribution.java    From astor with GNU General Public License v2.0 4 votes vote down vote up
/** {@inheritDoc} */
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);
}