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

Introduction

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

Prototype

public static double logGamma(double x)

Source Link

Document

Returns the value of log Γ(x) for x > 0.

Usage

From source file:beast.math.distribution.GammaDistributionTest.java

/**
 * This test stochastically draws gamma// w  w  w.j a v  a  2s  . c om
 * variates and compares the coded pdf 
 * with the actual pdf.  
 * The tolerance is required to be at most 1e-10.
 */

static double mypdf(double value, double shape, double scale) {
    return Math.exp(
            (shape - 1) * Math.log(value) - value / scale - Gamma.logGamma(shape) - shape * Math.log(scale));
}

From source file:com.opengamma.strata.math.impl.function.special.NaturalLogGammaFunction.java

@Override
public Double apply(Double x) {
    ArgChecker.isTrue(x > 0, "x must be greater than zero");
    return Gamma.logGamma(x);
}

From source file:edu.byu.nlp.math.GammaFunctions.java

/**
 * Computes ln Gamma(numerator) - ln Gamma(denominator).
 *//*ww  w .ja  v  a 2s .c  om*/
public static double logRatioOfGammas(double numerator, double denominator) {
    // There are potential speed ups when the difference is integral
    double diff = numerator - denominator;
    if (Math2.isIntegral(diff, IS_INTEGER_EPS)) {
        // When n = numerator - denominator is an integer, then 
        // Gamma(numerator) / Gamma(denominator) = x * (x + 1) * ... * (x + n - 1)
        // Which is, by definition the rising factorial of x^(n)
        return logRisingFactorial(denominator, (int) (Math.round(diff)));
    }
    return Gamma.logGamma(numerator) - Gamma.logGamma(denominator);
}

From source file:gedi.core.data.reads.functions.ReadDirichletLikelihoodRatioTest.java

public static double logProbability(double[] alpha1, double[] p) {
    double re = 0;
    double asum = 0;
    double gsum = 0;
    for (int i = 0; i < p.length; i++) {
        re += Math.log(p[i]) * (alpha1[i]);
        asum += alpha1[i] + 1;/* w  w w .  j  a v  a 2s . c o  m*/
        gsum += Gamma.logGamma(alpha1[i] + 1);
    }
    return re + Gamma.logGamma(asum) - gsum;
}

From source file:com.opengamma.strata.math.impl.function.special.GammaFunction.java

@Override
public double applyAsDouble(double x) {
    if (x > 0.0) {
        return Math.exp(Gamma.logGamma(x));
    }//from   w w w  .  ja  v  a2 s.  c o  m
    return Math.PI / Math.sin(Math.PI * x) / applyAsDouble(1 - x);
}

From source file:bide.prior.PriorBeta.java

private static double recomputeZ(double alpha, double beta) {

    return (Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta));

}

From source file:edu.byu.nlp.stats.RandomGeneratorsTest.java

/**
 * Dirichlet Compounad Multinomial where \alpha_k = 1 for all k. 
 * // www. j ava 2  s  .  c om
 * NOTE: this DCM is based on iid samples from a Categorical distribution and hence lacks the
 * multinomial factorial.
 */
// TODO(rhaertel): Implement a DCM class and just call that
private double dcmLogDensityWithUniformPrior(int[] observations) {
    double sumAlpha = observations.length;
    double logDensity = Gamma.logGamma(sumAlpha) - Gamma.logGamma(IntArrays.sum(observations) + sumAlpha);
    for (int i = 0; i < observations.length; i++) {
        logDensity += Gamma.logGamma(observations[i] + 1.0) - Gamma.logGamma(1.0);
    }
    return logDensity;
}

From source file:edu.byu.nlp.stats.SymmetricDirichletMLEOptimizableTest.java

/**
 * Test method for {@link edu.byu.nlp.stats.SymmetricDirichletMLENROptimizable#computeNext(java.lang.Double)}.
 *///from w  w w  .  j a  v  a 2 s.co  m
@Test
public void testComputeNext() {
    final double alpha = 2.13;

    final double[][] data = DirichletTestUtils.sampleDataset();
    SymmetricDirichletMLENROptimizable o = SymmetricDirichletMLENROptimizable.newOptimizable(data);

    ValueAndObject<Double> vao = o.computeNext(alpha);

    // Compute the update that should have happened
    final int N = data.length;
    final int K = data[0].length;
    double sumOfLogThetas = DoubleArrays.sum(Matrices.sumOverFirst(data));
    double value = N * Gamma.logGamma(K * alpha) - N * K * Gamma.logGamma(alpha) + (alpha - 1) * sumOfLogThetas;
    double firstDeriv = N * Gamma.digamma(K * alpha) - N * K * Gamma.digamma(alpha) + sumOfLogThetas;
    double secondDeriv = N * Gamma.trigamma(K * alpha) - N * K * Gamma.trigamma(alpha);
    double nextAlpha = alpha - firstDeriv / secondDeriv;

    // Ensure that the update was computed correctly
    assertThat(vao.getValue()).isEqualTo(value, delta);
    assertThat(vao.getObject().doubleValue()).isEqualTo(nextAlpha, delta);
}

From source file:edu.byu.nlp.math.GammaFunctions.java

public static double logRisingFactorial(double x, int diff) {
    Preconditions.checkArgument(diff >= 0);
    Preconditions.checkArgument(x >= 0.0);
    //    if (diff < -0.0) {
    //      throw new IllegalArgumentException("Argument must be positive; was " + diff);
    //    }// w  w  w  .j a  v  a2s.  c om

    // We can use a falling factorial when the difference is small.
    if (diff < MANUALLY_COMPUTE_RISING_FACTORIAL_THRESHOLD) {
        return logRisingFactorialManual(x, diff);
    }
    return Gamma.logGamma(x + diff) - Gamma.logGamma(x);
}

From source file:edu.byu.nlp.stats.SymmetricDirichletOptimizationHelper.java

public double valueAt(double alpha) {
    Preconditions.checkArgument(alpha > 0.0, "alpha must be strictly greater than zero; was %s", alpha);
    double value = N * (Gamma.logGamma(K * alpha) - K * Gamma.logGamma(alpha)) + (alpha - 1) * sumOfLogThetas;
    logger.debug(String.format("valueAt(%f) = %f", alpha, value));
    return value;
}