Java Code Examples for org.apache.commons.math.util.FastMath#log()

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

/**
 * A part of the deviance portion of the saddle point approximation.
 * <p>
 * References:
 * <ol>
 * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
 * Probabilities.". <a target="_blank"
 * href="http://www.herine.net/stat/papers/dbinom.pdf">
 * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
 * </ol>
 * </p>
 *
 * @param x the x value.
 * @param mu the average.
 * @return a part of the deviance.
 */
static double getDeviancePart(double x, double mu) {
    double ret;
    if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
        double d = x - mu;
        double v = d / (x + mu);
        double s1 = v * d;
        double s = Double.NaN;
        double ej = 2.0 * x * v;
        v = v * v;
        int j = 1;
        while (s1 != s) {
            s = s1;
            ej *= v;
            s1 = s + ej / ((j * 2) + 1);
            ++j;
        }
        ret = s1;
    } else {
        ret = x * FastMath.log(x / mu) + mu - x;
    }
    return ret;
}
 

/**
 * For this distribution, X, this method returns the critical point x, such
 * that P(X &lt; x) = <code>p</code>.
 * <p>
 * Returns 0 for p=0 and <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
 *
 * @param p the desired probability
 * @return x, such that P(X &lt; x) = <code>p</code>
 * @throws MathException if the inverse cumulative probability can not be
 *            computed due to convergence or other numerical errors.
 * @throws IllegalArgumentException if p < 0 or p > 1.
 */
@Override
public double inverseCumulativeProbability(double p) throws MathException {
    double ret;

    if (p < 0.0 || p > 1.0) {
        throw MathRuntimeException.createIllegalArgumentException(
              LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
    } else if (p == 1.0) {
        ret = Double.POSITIVE_INFINITY;
    } else {
        ret = -mean * FastMath.log(1.0 - p);
    }

    return ret;
}
 

/**
 * Return the probability density for a particular point.
 *
 * @param x The point at which the density should be computed.
 * @return The pdf at point x.
 * @since 2.1
 */
public double density(double x) {
    recomputeZ();
    if (x < 0 || x > 1) {
        return 0;
    } else if (x == 0) {
        if (alpha < 1) {
            throw MathRuntimeException.createIllegalArgumentException(
                    LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha);
        }
        return 0;
    } else if (x == 1) {
        if (beta < 1) {
            throw MathRuntimeException.createIllegalArgumentException(
                    LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta);
        }
        return 0;
    } else {
        double logX = FastMath.log(x);
        double log1mX = FastMath.log1p(-x);
        return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z);
    }
}
 

@Test
public void testGradientComponent5() {
    final double m = 1.2;
    final double k = 3.4;
    final double a = 2.3;
    final double q = 0.567;
    final double b = -FastMath.log(q);
    final double n = 3.4;

    final Logistic.Parametric f = new Logistic.Parametric();
    
    final double x = m - 1;
    final double qExp1 = 2;

    final double[] gf = f.gradient(x, new double[] {k, m, b, q, a, n});

    Assert.assertEquals((k - a) * FastMath.log(qExp1) / (n * n * FastMath.pow(qExp1, 1 / n)),
                        gf[5], EPS);
}
 

/**
 * Compute the PMF for a binomial distribution using the saddle point
 * expansion.
 *
 * @param x the value at which the probability is evaluated.
 * @param n the number of trials.
 * @param p the probability of success.
 * @param q the probability of failure (1 - p).
 * @return log(p(x)).
 */
static double logBinomialProbability(int x, int n, double p, double q) {
    double ret;
    if (x == 0) {
        if (p < 0.1) {
            ret = -getDeviancePart(n, n * q) - n * p;
        } else {
            ret = n * FastMath.log(q);
        }
    } else if (x == n) {
        if (q < 0.1) {
            ret = -getDeviancePart(n, n * p) - n * q;
        } else {
            ret = n * FastMath.log(p);
        }
    } else {
        ret = getStirlingError(n) - getStirlingError(x) -
              getStirlingError(n - x) - getDeviancePart(x, n * p) -
              getDeviancePart(n - x, n * q);
        double f = (MathUtils.TWO_PI * x * (n - x)) / n;
        ret = -0.5 * FastMath.log(f) + ret;
    }
    return ret;
}
 

static double getStirlingError(double z) {
    double ret;
    double z2;
    if (z < 15.0D) {
        z2 = 2.0D * z;
        if (FastMath.floor(z2) == z2) {
            ret = EXACT_STIRLING_ERRORS[(int)z2];
        } else {
            ret = Gamma.logGamma(z + 1.0D) - (z + 0.5D) * FastMath.log(z) + z - HALF_LOG_2_PI;
        }
    } else {
        z2 = z * z;
        ret = (0.08333333333333333D - (0.002777777777777778D - (7.936507936507937E-4D - (5.952380952380953E-4D - 8.417508417508417E-4D / z2) / z2) / z2) / z2) / z;
    }

    return ret;
}
 

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

/**
 * A part of the deviance portion of the saddle point approximation.
 * <p>
 * References:
 * <ol>
 * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial
 * Probabilities.". <a target="_blank"
 * href="http://www.herine.net/stat/papers/dbinom.pdf">
 * http://www.herine.net/stat/papers/dbinom.pdf</a></li>
 * </ol>
 * </p>
 *
 * @param x the x value.
 * @param mu the average.
 * @return a part of the deviance.
 */
static double getDeviancePart(double x, double mu) {
    double ret;
    if (FastMath.abs(x - mu) < 0.1 * (x + mu)) {
        double d = x - mu;
        double v = d / (x + mu);
        double s1 = v * d;
        double s = Double.NaN;
        double ej = 2.0 * x * v;
        v = v * v;
        int j = 1;
        while (s1 != s) {
            s = s1;
            ej *= v;
            s1 = s + ej / ((j * 2) + 1);
            ++j;
        }
        ret = s1;
    } else {
        ret = x * FastMath.log(x / mu) + mu - x;
    }
    return ret;
}
 

static double logBinomialProbability(int x, int n, double p, double q) {
    double ret;
    if (x == 0) {
        if (p < 0.1D) {
            ret = -getDeviancePart((double)n, (double)n * q) - (double)n * p;
        } else {
            ret = (double)n * FastMath.log(q);
        }
    } else if (x == n) {
        if (q < 0.1D) {
            ret = -getDeviancePart((double)n, (double)n * p) - (double)n * q;
        } else {
            ret = (double)n * FastMath.log(p);
        }
    } else {
        ret = getStirlingError((double)n) - getStirlingError((double)x) - getStirlingError((double)(n - x)) - getDeviancePart((double)x, (double)n * p) - getDeviancePart((double)(n - x), (double)n * q);
        double f = 6.283185307179586D * (double)x * (double)(n - x) / (double)n;
        ret += -0.5D * FastMath.log(f);
    }

    return ret;
}
 

@Test
public void testGradientComponent5() {
    final double m = 1.2;
    final double k = 3.4;
    final double a = 2.3;
    final double q = 0.567;
    final double b = -FastMath.log(q);
    final double n = 3.4;

    final Logistic.Parametric f = new Logistic.Parametric();
    
    final double x = m - 1;
    final double qExp1 = 2;

    final double[] gf = f.gradient(x, new double[] {k, m, b, q, a, n});

    Assert.assertEquals((k - a) * FastMath.log(qExp1) / (n * n * FastMath.pow(qExp1, 1 / n)),
                        gf[5], EPS);
}
 

/**
 * For this distribution, X, this method returns the critical point x, such
 * that {@code P(X < x) = p}.
 * It will return 0 when p = 0 and {@code Double.POSITIVE_INFINITY}
 * when p = 1.
 *
 * @param p Desired probability.
 * @return {@code x}, such that {@code P(X < x) = p}.
 * @throws MathException if the inverse cumulative probability can not be
 * computed due to convergence or other numerical errors.
 * @throws OutOfRangeException if {@code p < 0} or {@code p > 1}.
 */
@Override
public double inverseCumulativeProbability(double p) throws MathException {
    double ret;

    if (p < 0.0 || p > 1.0) {
        throw new OutOfRangeException(p, 0.0, 1.0);
    } else if (p == 1.0) {
        ret = Double.POSITIVE_INFINITY;
    } else {
        ret = -mean * FastMath.log(1.0 - p);
    }

    return ret;
}
 

/**
 * @param x Value at which to compute the logit.
 * @param lo Lower bound.
 * @param hi Higher bound.
 * @return the value of the logit function at {@code x}.
 */
private static double value(double x,
                            double lo,
                            double hi) {
    if (x < lo || x > hi) {
        throw new OutOfRangeException(x, lo, hi);
    }
    return FastMath.log((x - lo) / (hi - x));
}
 

/**
 * Returns the probability density for a particular point.
 *
 * @param x The point at which the density should be computed.
 * @return The pdf at point x.
 * @since 2.1
 */
@Override
public double density(double x) {
    final double nhalf = numeratorDegreesOfFreedom / 2;
    final double mhalf = denominatorDegreesOfFreedom / 2;
    final double logx = FastMath.log(x);
    final double logn = FastMath.log(numeratorDegreesOfFreedom);
    final double logm = FastMath.log(denominatorDegreesOfFreedom);
    final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x + denominatorDegreesOfFreedom);
    return FastMath.exp(nhalf*logn + nhalf*logx - logx + mhalf*logm - nhalf*lognxm -
           mhalf*lognxm - Beta.logBeta(nhalf, mhalf));
}
 

/** {@inheritDoc} */
@Override
public double value(double d) {
    return FastMath.log(d);
}
 

/** {@inheritDoc} */
public double value(double x) {
    return FastMath.log(x);
}
 

/**
 * {@inheritDoc}
 */
@Override
public void increment(final double d) {
    value += FastMath.log(d);
    n++;
}
 

/**
 * Returns the sum of the natural logs of the entries in the specified portion of
 * the input array, or <code>Double.NaN</code> if the designated subarray
 * is empty.
 * <p>
 * Throws <code>IllegalArgumentException</code> if the array is null.</p>
 * <p>
 * See {@link SumOfLogs}.</p>
 *
 * @param values the input array
 * @param begin index of the first array element to include
 * @param length the number of elements to include
 * @return the sum of the natural logs of the values or 0 if
 * length = 0
 * @throws IllegalArgumentException if the array is null or the array index
 *  parameters are not valid
 */
@Override
public double evaluate(final double[] values, final int begin, final int length) {
    double sumLog = Double.NaN;
    if (test(values, begin, length, true)) {
        sumLog = 0.0;
        for (int i = begin; i < begin + length; i++) {
            sumLog += FastMath.log(values[i]);
        }
    }
    return sumLog;
}
 

/**
 * Returns the sum of the natural logs of the entries in the specified portion of
 * the input array, or <code>Double.NaN</code> if the designated subarray
 * is empty.
 * <p>
 * Throws <code>IllegalArgumentException</code> if the array is null.</p>
 * <p>
 * See {@link SumOfLogs}.</p>
 *
 * @param values the input array
 * @param begin index of the first array element to include
 * @param length the number of elements to include
 * @return the sum of the natural logs of the values or Double.NaN if
 * length = 0
 * @throws IllegalArgumentException if the array is null or the array index
 *  parameters are not valid
 */
@Override
public double evaluate(final double[] values, final int begin, final int length) {
    double sumLog = Double.NaN;
    if (test(values, begin, length)) {
        sumLog = 0.0;
        for (int i = begin; i < begin + length; i++) {
            sumLog += FastMath.log(values[i]);
        }
    }
    return sumLog;
}
 

/**
 * Returns a random value from an Exponential distribution with the given
 * mean.
 * <p>
 * <strong>Algorithm Description</strong>: Uses the <a
 * href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html"> Inversion
 * Method</a> to generate exponentially distributed random values from
 * uniform deviates.
 * </p>
 *
 * @param mean the mean of the distribution
 * @return the random Exponential value
 * @throws NotStrictlyPositiveException if {@code mean <= 0}.
 */
public double nextExponential(double mean) {
    if (mean <= 0.0) {
        throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
    }
    final RandomGenerator generator = getRan();
    double unif = generator.nextDouble();
    while (unif == 0.0d) {
        unif = generator.nextDouble();
    }
    return -mean * FastMath.log(unif);
}