org.apache.commons.math3.special.Gamma Java Examples

@Test
public void testMoments() {
    final double tol = 1e-9;
    WeibullDistribution dist;

    dist = new WeibullDistribution(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);

    dist = new WeibullDistribution(10.4, 2.222);
    Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

/**
 * Calculate the posterior probability that a single biallelic genotype is non-ref
 *
 * The nth genotype (n runs from 0 to the sample ploidy, inclusive) contains n copies of the alt allele
 * @param log10GenotypeLikelihoods
 * @return
 */
public double calculateSingleSampleBiallelicNonRefPosterior(final double[] log10GenotypeLikelihoods, final boolean returnZeroIfRefIsMax) {
    Utils.nonNull(log10GenotypeLikelihoods);

    if (returnZeroIfRefIsMax && MathUtils.maxElementIndex(log10GenotypeLikelihoods) == 0) {
        return 0;
    }

    final int ploidy = log10GenotypeLikelihoods.length - 1;

    final double[] log10UnnormalizedPosteriors = new IndexRange(0, ploidy + 1)
            .mapToDouble(n -> log10GenotypeLikelihoods[n] + MathUtils.log10BinomialCoefficient(ploidy, n)
                    + MathUtils.logToLog10(Gamma.logGamma(n + snpPseudocount ) + Gamma.logGamma(ploidy - n + refPseudocount)));

    return (returnZeroIfRefIsMax && MathUtils.maxElementIndex(log10UnnormalizedPosteriors) == 0) ? 0.0 :
            1 - MathUtils.normalizeFromLog10ToLinearSpace(log10UnnormalizedPosteriors)[0];
}
 

@Test
public void testMoments() {
    final double tol = 1e-9;
    WeibullDistribution dist;

    dist = new WeibullDistribution(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);

    dist = new WeibullDistribution(10.4, 2.222);
    Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

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

@Test
public void testMoments() {
    final double tol = 1e-9;
    WeibullDistribution dist;

    dist = new WeibullDistribution(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);

    dist = new WeibullDistribution(10.4, 2.222);
    Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

public void precomputeDenominatorForVariationalBayes( final double sumHyperParameterLambda ) {

        // Variational Bayes calculations from Bishop
        precomputeInverse();
        cachedSigmaInverse.timesEquals( hyperParameter_a );
        double sum = 0.0;
        for(int jjj = 1; jjj <= mu.length; jjj++) {
            sum += Gamma.digamma( (hyperParameter_a + 1.0 - jjj) / 2.0 );
        }
        sum -= Math.log( sigma.det() );
        sum += Math.log(2.0) * mu.length;
        final double lambda = 0.5 * sum;
        final double pi = Gamma.digamma( hyperParameter_lambda ) - Gamma.digamma( sumHyperParameterLambda );
        final double beta = (-1.0 * mu.length) / (2.0 * hyperParameter_b);
        cachedDenomLog10 = (pi / Math.log(10.0)) + (lambda / Math.log(10.0)) + (beta / Math.log(10.0));
    }
 

@Test
public void testHetLogLikelihoodMinorFractionNearOne() {
    final double pi = 0.01; //pi is just a prefactor so we don't need to test it thoroughly here
    for (final double f : Arrays.asList(1 - 1e-6, 1 - 1e-7, 1 - 1e-8)) {
        for (final double mean : Arrays.asList(0.9, 1.0, 1.1)) {
            for (final double variance : Arrays.asList(0.01, 0.005, 0.001)) {
                final double alpha = mean * mean / variance;
                final double beta = mean / variance;
                final AlleleFractionGlobalParameters parameters = new AlleleFractionGlobalParameters(mean, variance, pi);
                for (final int a : Arrays.asList(1, 10, 20)) {  //alt count
                    for (final int r : Arrays.asList(1, 10, 20)) { //ref count
                        final AllelicCount count = new AllelicCount(DUMMY, r, a);
                        final double actual = AlleleFractionLikelihoods.hetLogLikelihood(parameters, f, count, AlleleFractionIndicator.ALT_MINOR);
                        final double expected = -r * log(beta) + Gamma.logGamma(alpha + r) - Gamma.logGamma(alpha)
                                + log((1 - pi) / 2) - r * log(f / (1 - f));
                        Assert.assertEquals(actual, expected,1e-4);
                    }
                }
            }
        }
    }
}
 

public MultivariateTDistribution(RealVector mean, RealMatrix covarianceMatrix, double degreesOfFreedom) {
    this.mean = mean;
    if (mean.getDimension() > 1) {
        this.precisionMatrix = MatrixUtils.blockInverse(covarianceMatrix, (-1 + covarianceMatrix.getColumnDimension()) / 2);
    } else {
        this.precisionMatrix = MatrixUtils.createRealIdentityMatrix(1).scalarMultiply(1. / covarianceMatrix.getEntry(0, 0));
    }
    this.dof = degreesOfFreedom;

    this.D = mean.getDimension();

    double determinant = new LUDecomposition(covarianceMatrix).getDeterminant();

    this.multiplier = Math.exp(Gamma.logGamma(0.5 * (D + dof)) - Gamma.logGamma(0.5 * dof)) /
            Math.pow(Math.PI * dof, 0.5 * D) /
            Math.pow(determinant, 0.5);
}
 

@Test
public void testFactorialLogCacheExpansion() {
    // There is no way to determine if the cache values were reused but this test
    // exercises the method to ensure it does not error.
    final FactorialLog factorialLog = FactorialLog.create()
                                                  // Edge case where cache should not be copied (<2)
                                                  .withCache(1)
                                                  // Expand
                                                  .withCache(5)
                                                  // Expand more
                                                  .withCache(10)
                                                  // Contract
                                                  .withCache(5);
    for (int n = 1; n <= 5; n++) {
        // Use Commons math to compute logGamma(1 + n);
        double expected = Gamma.logGamma(1 + n);
        Assert.assertEquals(expected, factorialLog.value(n), 1e-10);
    }
}
 

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

@Override
public void learn(final List<Datum> data, final double[] responsibilities) {
    double alpha = betaDistributionShape.getAlpha();
    double beta = betaDistributionShape.getBeta();

    for (int epoch = 0; epoch < NUM_EPOCHS; epoch++) {
        for (int n = 0; n < data.size(); n++) {
            final Datum datum = data.get(n);
            final int alt = datum.getAltCount();
            final int ref = datum.getTotalCount() - alt;

            final double digammaOfTotalPlusAlphaPlusBeta = Gamma.digamma(datum.getTotalCount() + alpha + beta);
            final double digammaOfAlphaPlusBeta = Gamma.digamma(alpha + beta);
            final double alphaGradient = Gamma.digamma(alpha + alt) - digammaOfTotalPlusAlphaPlusBeta - Gamma.digamma(alpha) + digammaOfAlphaPlusBeta;
            final double betaGradient = Gamma.digamma(beta + ref) - digammaOfTotalPlusAlphaPlusBeta - Gamma.digamma(beta) + digammaOfAlphaPlusBeta;

            alpha = Math.max(alpha + RATE * alphaGradient * responsibilities[n], 1.0);
            beta = Math.max(beta + RATE * betaGradient * responsibilities[n], 0.5);
        }
    }

    betaDistributionShape = new BetaDistributionShape(alpha, beta);

}
 

@Test
public void testMoments() {
    final double tol = 1e-9;
    WeibullDistribution dist;

    dist = new WeibullDistribution(2.5, 3.5);
    // In R: 3.5*gamma(1+(1/2.5)) (or emperically: mean(rweibull(10000, 2.5, 3.5)))
    Assert.assertEquals(dist.getNumericalMean(), 3.5 * FastMath.exp(Gamma.logGamma(1 + (1 / 2.5))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (3.5 * 3.5) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 2.5))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);

    dist = new WeibullDistribution(10.4, 2.222);
    Assert.assertEquals(dist.getNumericalMean(), 2.222 * FastMath.exp(Gamma.logGamma(1 + (1 / 10.4))), tol);
    Assert.assertEquals(dist.getNumericalVariance(), (2.222 * 2.222) *
            FastMath.exp(Gamma.logGamma(1 + (2 / 10.4))) -
            (dist.getNumericalMean() * dist.getNumericalMean()), tol);
}
 

@Test
// TODO 通过https://www.wolframalpha.com/input验证,Apache比较准确.
public void testTrigamma() {
    // trigamma遇到负整数会变为NaN或者无穷.
    Assert.assertTrue(Float.isNaN(GammaUtility.trigamma(0F)) == Double.isNaN(Gamma.trigamma(0D)));
    Assert.assertTrue(Float.isNaN(GammaUtility.trigamma(-0F)) == Double.isNaN(Gamma.trigamma(-0D)));
    Assert.assertTrue(Float.isNaN(GammaUtility.trigamma(-1F)) == Double.isNaN(Gamma.trigamma(-1D)));

    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(0.1F), (float) Gamma.trigamma(0.1D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(0.5F), (float) Gamma.trigamma(0.5D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(1F), (float) Gamma.trigamma(1D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(1.5F), (float) Gamma.trigamma(1.5D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(8F), (float) Gamma.trigamma(8D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(8.1F), (float) Gamma.trigamma(8.1D)));

    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(-0.1F), (float) Gamma.trigamma(-0.1D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(-0.5F), (float) Gamma.trigamma(-0.5D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(-1.5F), (float) Gamma.trigamma(-1.5D)));
    Assert.assertTrue(MathUtility.equal(GammaUtility.trigamma(-8.1F), (float) Gamma.trigamma(-8.1D)));
}
 

/**
 * {@inheritDoc}
 */
@Override
public SufficientStatistics createInitSufficientStatistics() {

    ArrayVector vector = new ArrayVector(this.sizeOfSufficientStatistics());

    double alpha = 1;
    double beta = 1;
    double mean = 0;
    double precision = 1e-10;

    this.momentParameters.set(EF_JointNormalGamma.MU_GAMMA, mean*alpha/beta);
    this.momentParameters.set(EF_JointNormalGamma.MUSQUARE_GAMMA, 1/precision + mean*mean*alpha/beta);
    this.momentParameters.set(EF_JointNormalGamma.GAMMA, alpha/beta);
    this.momentParameters.set(EF_JointNormalGamma.LOGGAMMA, Gamma.digamma(alpha) - Math.log(beta));

    return vector;
}
 

/** {@inheritDoc} */
public double density(double x) {
    final double n = degreesOfFreedom;
    final double nPlus1Over2 = (n + 1) / 2;
    return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
                        0.5 * (FastMath.log(FastMath.PI) +
                               FastMath.log(n)) -
                        Gamma.logGamma(n / 2) -
                        nPlus1Over2 * FastMath.log(1 + x * x / n));
}
 

/**
 * {@inheritDoc}
 */
@Override
public SufficientStatistics createInitSufficientStatistics() {

    ArrayVector vector = new ArrayVector(this.sizeOfSufficientStatistics());

    double alpha = 1.0;
    double beta = 1.0;
    vector.set(0, Math.log(beta) - Gamma.digamma(alpha));
    vector.set(1, alpha / beta);

    return vector;
}
 

/**
 * Overloaded digamma function: returns 0 if argument is zero
 *
 * @param x
 * @return digamma(x) if x nonzero, otherwise zero
 */
private double digamma(final double x) {
	if (Math.abs(x) < 1e-15)
		return 0.0;
	else
		return Gamma.digamma(x);
}
 

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

/** {@inheritDoc} */
public double density(double x) {
    final double n = degreesOfFreedom;
    final double nPlus1Over2 = (n + 1) / 2;
    return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
                        0.5 * (FastMath.log(FastMath.PI) +
                               FastMath.log(n)) -
                        Gamma.logGamma(n / 2) -
                        nPlus1Over2 * FastMath.log(1 + x * x / n));
}
 

/**
 * used by {@link #getNumericalMean()}
 *
 * @return the mean of this distribution
 */
protected double calculateNumericalMean() {
    final double sh = getShape();
    final double sc = getScale();

    return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
}
 

/** {@inheritDoc} */
public double density(double x) {
    final double n = degreesOfFreedom;
    final double nPlus1Over2 = (n + 1) / 2;
    return FastMath.exp(Gamma.logGamma(nPlus1Over2) -
                        0.5 * (FastMath.log(FastMath.PI) +
                               FastMath.log(n)) -
                        Gamma.logGamma(n / 2) -
                        nPlus1Over2 * FastMath.log(1 + x * x / n));
}
 

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

/**
 * used by {@link #getNumericalMean()}
 *
 * @return the mean of this distribution
 */
protected double calculateNumericalMean() {
    final double sh = getShape();
    final double sc = getScale();

    return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
}
 

/** {@inheritDoc} */
public double cumulativeProbability(int x) {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
                                   maxIterations);
}
 

/**
 * used by {@link #getNumericalMean()}
 *
 * @return the mean of this distribution
 */
protected double calculateNumericalMean() {
    final double sh = getShape();
    final double sc = getScale();

    return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
}
 

/**
 * used by {@link #getNumericalVariance()}
 *
 * @return the variance of this distribution
 */
protected double calculateNumericalVariance() {
    final double sh = getShape();
    final double sc = getScale();
    final double mn = getNumericalMean();

    return (sc * sc) * FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
           (mn * mn);
}
 

/**
 * used by {@link #getNumericalMean()}
 *
 * @return the mean of this distribution
 */
protected double calculateNumericalMean() {
    final double sh = getShape();
    final double sc = getScale();

    return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
}
 

/** {@inheritDoc} */
public double density(double x) {
    if (x <= 0) {
        return 0.0;
    }
    return 2.0 * FastMath.pow(mu, mu) / (Gamma.gamma(mu) * FastMath.pow(omega, mu)) *
                 FastMath.pow(x, 2 * mu - 1) * FastMath.exp(-mu * x * x / omega);
}
 

/**
 * {@inheritDoc}
 */
@Override
public double computeLogNormalizer() {

    double sumOfU_i = 0;
    double sumLogGammaOfU_i = 0;
    for (int i = 0; i < nOfStates; i++) {
        sumOfU_i += naturalParameters.get(i);
        sumLogGammaOfU_i += Gamma.logGamma(naturalParameters.get(i));
    }

    return sumLogGammaOfU_i - Gamma.logGamma(sumOfU_i);
}
 

/** {@inheritDoc} */
public double cumulativeProbability(int x) {
    if (x < 0) {
        return 0;
    }
    if (x == Integer.MAX_VALUE) {
        return 1;
    }
    return Gamma.regularizedGammaQ((double) x + 1, mean, epsilon,
                                   maxIterations);
}