Java Code Examples for org.apache.commons.math3.util.FastMath#log1p()

/** Computes shifted logarithm of a derivative structure.
 * @param operand array holding the operand
 * @param operandOffset offset of the operand in its array
 * @param result array where result must be stored (for
 * shifted logarithm the result array <em>cannot</em> be the input array)
 * @param resultOffset offset of the result in its array
 */
public void log1p(final double[] operand, final int operandOffset,
                  final double[] result, final int resultOffset) {

    // create the function value and derivatives
    double[] function = new double[1 + order];
    function[0] = FastMath.log1p(operand[operandOffset]);
    if (order > 0) {
        double inv = 1.0 / (1.0 + operand[operandOffset]);
        double xk  = inv;
        for (int i = 1; i <= order; ++i) {
            function[i] = xk;
            xk *= -i * inv;
        }
    }

    // apply function composition
    compose(operand, operandOffset, function, result, resultOffset);

}
 

/** {@inheritDoc} */
public double density(double x) {
    recomputeZ();
    if (x < 0 || x > 1) {
        return 0;
    } else if (x == 0) {
        if (alpha < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
        }
        return 0;
    } else if (x == 1) {
        if (beta < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
        }
        return 0;
    } else {
        double logX = FastMath.log(x);
        double log1mX = FastMath.log1p(-x);
        return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
    }
}
 

/**
 * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
 * <em>NSWC Library of Mathematics Subroutines</em> double precision
 * implementation, {@code DGSMLN}. In {@link BetaTest#testLogGammaSum()},
 * this private method is accessed through reflection.
 *
 * @param a First argument.
 * @param b Second argument.
 * @return the value of {@code log(Gamma(a + b))}.
 * @throws OutOfRangeException if {@code a} or {@code b} is lower than
 * {@code 1.0} or greater than {@code 2.0}.
 */
private static double logGammaSum(final double a, final double b)
    throws OutOfRangeException {

    if ((a < 1.0) || (a > 2.0)) {
        throw new OutOfRangeException(a, 1.0, 2.0);
    }
    if ((b < 1.0) || (b > 2.0)) {
        throw new OutOfRangeException(b, 1.0, 2.0);
    }

    final double x = (a - 1.0) + (b - 1.0);
    if (x <= 0.5) {
        return Gamma.logGamma1p(1.0 + x);
    } else if (x <= 1.5) {
        return Gamma.logGamma1p(x) + FastMath.log1p(x);
    } else {
        return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
    }
}
 

/** Computes shifted logarithm of a derivative structure.
 * @param operand array holding the operand
 * @param operandOffset offset of the operand in its array
 * @param result array where result must be stored (for
 * shifted logarithm the result array <em>cannot</em> be the input array)
 * @param resultOffset offset of the result in its array
 */
public void log1p(final double[] operand, final int operandOffset,
                  final double[] result, final int resultOffset) {

    // create the function value and derivatives
    double[] function = new double[1 + order];
    function[0] = FastMath.log1p(operand[operandOffset]);
    if (order > 0) {
        double inv = 1.0 / (1.0 + operand[operandOffset]);
        double xk  = inv;
        for (int i = 1; i <= order; ++i) {
            function[i] = xk;
            xk *= -i * inv;
        }
    }

    // apply function composition
    compose(operand, operandOffset, function, result, resultOffset);

}
 

/** {@inheritDoc} */
public double density(double x) {
    recomputeZ();
    if (x < 0 || x > 1) {
        return 0;
    } else if (x == 0) {
        if (alpha < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
        }
        return 0;
    } else if (x == 1) {
        if (beta < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
        }
        return 0;
    } else {
        double logX = FastMath.log(x);
        double log1mX = FastMath.log1p(-x);
        return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
    }
}
 

/** {@inheritDoc} */
public double density(double x) {
    recomputeZ();
    if (x < 0 || x > 1) {
        return 0;
    } else if (x == 0) {
        if (alpha < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
        }
        return 0;
    } else if (x == 1) {
        if (beta < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
        }
        return 0;
    } else {
        double logX = FastMath.log(x);
        double log1mX = FastMath.log1p(-x);
        return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
    }
}
 

/**
 * Returns the value of log Γ(a + b) for 1 ≤ a, b ≤ 2. Based on the
 * <em>NSWC Library of Mathematics Subroutines</em> double precision
 * implementation, {@code DGSMLN}. In {@code BetaTest.testLogGammaSum()},
 * this private method is accessed through reflection.
 *
 * @param a First argument.
 * @param b Second argument.
 * @return the value of {@code log(Gamma(a + b))}.
 * @throws OutOfRangeException if {@code a} or {@code b} is lower than
 * {@code 1.0} or greater than {@code 2.0}.
 */
private static double logGammaSum(final double a, final double b)
    throws OutOfRangeException {

    if ((a < 1.0) || (a > 2.0)) {
        throw new OutOfRangeException(a, 1.0, 2.0);
    }
    if ((b < 1.0) || (b > 2.0)) {
        throw new OutOfRangeException(b, 1.0, 2.0);
    }

    final double x = (a - 1.0) + (b - 1.0);
    if (x <= 0.5) {
        return Gamma.logGamma1p(1.0 + x);
    } else if (x <= 1.5) {
        return Gamma.logGamma1p(x) + FastMath.log1p(x);
    } else {
        return Gamma.logGamma1p(x - 1.0) + FastMath.log(x * (1.0 + x));
    }
}
 

/** Computes shifted logarithm of a derivative structure.
 * @param operand array holding the operand
 * @param operandOffset offset of the operand in its array
 * @param result array where result must be stored (for
 * shifted logarithm the result array <em>cannot</em> be the input array)
 * @param resultOffset offset of the result in its array
 */
public void log1p(final double[] operand, final int operandOffset,
                  final double[] result, final int resultOffset) {

    // create the function value and derivatives
    double[] function = new double[1 + order];
    function[0] = FastMath.log1p(operand[operandOffset]);
    if (order > 0) {
        double inv = 1.0 / (1.0 + operand[operandOffset]);
        double xk  = inv;
        for (int i = 1; i <= order; ++i) {
            function[i] = xk;
            xk *= -i * inv;
        }
    }

    // apply function composition
    compose(operand, operandOffset, function, result, resultOffset);

}
 

/** {@inheritDoc} **/
@Override
public double logDensity(double x) {
    recomputeZ();
    if (x < 0 || x > 1) {
        return Double.NEGATIVE_INFINITY;
    } else if (x == 0) {
        if (alpha < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
        }
        return Double.NEGATIVE_INFINITY;
    } else if (x == 1) {
        if (beta < 1) {
            throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
        }
        return Double.NEGATIVE_INFINITY;
    } else {
        double logX = FastMath.log(x);
        double log1mX = FastMath.log1p(-x);
        return (alpha - 1) * logX + (beta - 1) * log1mX - z;
    }
}
 

/** Computes shifted logarithm of a derivative structure.
 * @param operand array holding the operand
 * @param operandOffset offset of the operand in its array
 * @param result array where result must be stored (for
 * shifted logarithm the result array <em>cannot</em> be the input array)
 * @param resultOffset offset of the result in its array
 */
public void log1p(final double[] operand, final int operandOffset,
                  final double[] result, final int resultOffset) {

    // create the function value and derivatives
    double[] function = new double[1 + order];
    function[0] = FastMath.log1p(operand[operandOffset]);
    if (order > 0) {
        double inv = 1.0 / (1.0 + operand[operandOffset]);
        double xk  = inv;
        for (int i = 1; i <= order; ++i) {
            function[i] = xk;
            xk *= -i * inv;
        }
    }

    // apply function composition
    compose(operand, operandOffset, function, result, resultOffset);

}
 

/**
 * Returns the value of log Γ(1 + x) for -0.5 ≤ x ≤ 1.5.
 * This implementation is based on the double precision implementation in
 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
 *
 * @param x Argument.
 * @return The value of {@code log(Gamma(1 + x))}.
 * @throws NumberIsTooSmallException if {@code x < -0.5}.
 * @throws NumberIsTooLargeException if {@code x > 1.5}.
 * @since 3.1
 */
public static double logGamma1p(final double x)
    throws NumberIsTooSmallException, NumberIsTooLargeException {

    if (x < -0.5) {
        throw new NumberIsTooSmallException(x, -0.5, true);
    }
    if (x > 1.5) {
        throw new NumberIsTooLargeException(x, 1.5, true);
    }

    return -FastMath.log1p(invGamma1pm1(x));
}
 

/** {@inheritDoc} */
@Override
public double logProbability(int x) {
    double ret;
    if (x < 0) {
        ret = Double.NEGATIVE_INFINITY;
    } else {
        final double p = probabilityOfSuccess;
        ret = x * FastMath.log1p(-p) + FastMath.log(p);
    }
    return ret;
}
 

/**
 * Returns the value of log Γ(1 + x) for -0.5 ≤ x ≤ 1.5.
 * This implementation is based on the double precision implementation in
 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}.
 *
 * @param x Argument.
 * @return The value of {@code log(Gamma(1 + x))}.
 * @throws NumberIsTooSmallException if {@code x < -0.5}.
 * @throws NumberIsTooLargeException if {@code x > 1.5}.
 */
public static double logGamma1p(final double x)
    throws NumberIsTooSmallException, NumberIsTooLargeException {

    if (x < -0.5) {
        throw new NumberIsTooSmallException(x, -0.5, true);
    }
    if (x > 1.5) {
        throw new NumberIsTooLargeException(x, 1.5, true);
    }

    return -FastMath.log1p(invGamma1pm1(x));
}
 

/**
 * Returns the value of log[Γ(b) / Γ(a + b)] for a ≥ 0 and b ≥ 10. Based on
 * the <em>NSWC Library of Mathematics Subroutines</em> double precision
 * implementation, {@code DLGDIV}. In
 * {@link BetaTest#testLogGammaMinusLogGammaSum()}, this private method is
 * accessed through reflection.
 *
 * @param a First argument.
 * @param b Second argument.
 * @return the value of {@code log(Gamma(b) / Gamma(a + b))}.
 * @throws NumberIsTooSmallException if {@code a < 0.0} or {@code b < 10.0}.
 */
private static double logGammaMinusLogGammaSum(final double a,
                                               final double b)
    throws NumberIsTooSmallException {

    if (a < 0.0) {
        throw new NumberIsTooSmallException(a, 0.0, true);
    }
    if (b < 10.0) {
        throw new NumberIsTooSmallException(b, 10.0, true);
    }

    /*
     * d = a + b - 0.5
     */
    final double d;
    final double w;
    if (a <= b) {
        d = b + (a - 0.5);
        w = deltaMinusDeltaSum(a, b);
    } else {
        d = a + (b - 0.5);
        w = deltaMinusDeltaSum(b, a);
    }

    final double u = d * FastMath.log1p(a / b);
    final double v = a * (FastMath.log(b) - 1.0);

    return u <= v ? (w - u) - v : (w - v) - u;
}
 

/** {@inheritDoc} */
@Override
public double logProbability(int x) {
    double ret;
    if (x < 0) {
        ret = Double.NEGATIVE_INFINITY;
    } else {
        final double p = probabilityOfSuccess;
        ret = x * FastMath.log1p(-p) + FastMath.log(p);
    }
    return ret;
}
 

/** {@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 SparseGradient log1p() {
    return new SparseGradient(FastMath.log1p(value), 1.0 / (1.0 + value), derivatives);
}
 

/** {@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 value(double x) {
    return FastMath.log1p(x);
}