Java Gamma.LANCZOS_G属性代码示例

/** {@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);
}
 

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;
}
 

/**
 * 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 CustomGammaDistribution(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);
}
 

/**
 * 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);
}
 

/**
 * 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);
}
 

/**
 * Copied from Apache Commons FastMath. Removed support for NaN and x < 0.5.
 * <p>
 * Returns the value of log&nbsp;&Gamma;(x) for x&nbsp;&gt;&nbsp;0.
 * </p>
 * <p>
 * For x &le; 8, the implementation is based on the double precision
 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
 * {@code DGAMLN}. For x &gt; 8, the implementation is based on
 * </p>
 * <ul>
 * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
 * Function</a>, equation (28).</li>
 * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
 * Lanczos Approximation</a>, equations (1) through (5).</li>
 * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
 * the computation of the convergent Lanczos complex Gamma
 * approximation</a></li>
 * </ul>
 *
 * @param x
 *            Argument.
 * @return the value of {@code log(Gamma(x))}, {@code Double.NaN} if
 *         {@code x <= 0.0}.
 */
private static double logGamma(double x)
{
	//		if (Double.isNaN(x) || (x <= 0.0))
	//		{
	//			return Double.NaN;
	//		}
	//		else if (x < 0.5)
	//		{
	//			return Gamma.logGamma1p(x) - Math.log(x);
	//		}
	//		else 
	if (x <= 2.5)
	{
		return Gamma.logGamma1p((x - 0.5) - 0.5);
	}
	else if (x <= 8.0)
	{
		final int n = (int) FastMath.floor(x - 1.5);
		double prod = 1.0;
		for (int i = 1; i <= n; i++)
		{
			prod *= x - i;
		}
		return Gamma.logGamma1p(x - (n + 1)) + Math.log(prod);
	}
	else
	{
		double sum = Gamma.lanczos(x);
		double tmp = x + Gamma.LANCZOS_G + .5;
		return ((x + .5) * Math.log(tmp)) - tmp + HALF_LOG_2_PI + Math.log(sum / x);
	}
}
 

public void init() {
    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);
}
 

/** {@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);
}
 

/** {@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);
}
 

/** {@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);
}