Java Gamma类代码示例

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
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));
    }
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/**
 * Creates a t distribution.
 *
 * @param rng Random number generator.
 * @param degreesOfFreedom Degrees of freedom.
 * @param inverseCumAccuracy the maximum absolute error in inverse
 * cumulative probability estimates
 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
 * @throws NotStrictlyPositiveException if {@code degreesOfFreedom <= 0}
 * @since 3.1
 */
public TDistribution(RandomGenerator rng,
                     double degreesOfFreedom,
                     double inverseCumAccuracy)
    throws NotStrictlyPositiveException {
    super(rng);

    if (degreesOfFreedom <= 0) {
        throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
                                               degreesOfFreedom);
    }
    this.degreesOfFreedom = degreesOfFreedom;
    solverAbsoluteAccuracy = inverseCumAccuracy;

    final double n = degreesOfFreedom;
    final double nPlus1Over2 = (n + 1) / 2;
    factor = Gamma.logGamma(nPlus1Over2) -
             0.5 * (FastMath.log(FastMath.PI) + FastMath.log(n)) -
             Gamma.logGamma(n / 2);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/** {@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);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
@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);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
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;
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/**
 * Computes log likelihood for word-topic vectors (P(w|z)) according to Eq. [2] in the reference.
 * <ul>
 *   <li>T: {@code numTopics}</li>
 *   <li>W: {@code numVocabs}</li>
 *   <li>n(j, w): <i>j</i>th topic's number of assignments to <i>w</i>th vocabulary</li>
 * </ul>

 * @return a portion of log likelihood computed from the given word-topic vectors
 */
double computeWordLLH(final Collection<int[]> wordTopicCounts, final int[] wordTopicCountsSummary) {
  double result = numTopics * (Gamma.logGamma(numVocabs * beta) - numVocabs * Gamma.logGamma(beta));
  for (final int[] wordTopicCount : wordTopicCounts) {
    // For computing log-likelihood, we need only the values. Please refer to SparseArrayCodec.
    for (int j = 1; j < wordTopicCount.length; j += 2) {
      result += Gamma.logGamma(wordTopicCount[j] + beta);
    }
    // handle the case of zero values separately
    result += logGammaBeta * (numTopics - wordTopicCount.length / 2);
  }
  for (int j = 1; j < wordTopicCountsSummary.length; j += 2) {
    result -= Gamma.logGamma(wordTopicCountsSummary[j] + numVocabs * beta);
  }
  // handle the case of zero values separately
  result -= Gamma.logGamma(numVocabs * beta) * (numTopics - wordTopicCountsSummary.length / 2);
  return result;
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
protected double[] calcEffectivePhis(final double E_alpha, final double[] responsibilitiesByRho) {

        final double sumResponsibilities = MathUtils.sum(responsibilitiesByRho);

        final double[] result = new double[responsibilitiesByRho.length];
        final int k = responsibilitiesByRho.length;

        // Initialize all pseudocounts to 1, except for index 0, which is 20;
        //  This artificially increases the odds for a rho = 0.
        final RealVector pseudocounts = new ArrayRealVector(responsibilitiesByRho.length);
        pseudocounts.set(1.0);
        pseudocounts.setEntry(0, 20.0);

        final double sumPseudocounts = MathUtils.sum(pseudocounts.toArray());
        final double term2 = Gamma.digamma(E_alpha + sumPseudocounts + sumResponsibilities);
        for (int i=0; i < result.length; i++) {
            final double term1 = Gamma.digamma(E_alpha/k + pseudocounts.getEntry(i) + responsibilitiesByRho[i]);
            result[i] = Math.exp(term1 - term2);
        }

        return result;
    }
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
@Test
public void testHetLogLikelihoodMinorFractionNearZero() {
    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(1e-6, 1e-7, 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, 2, 3)) {  //alt count
                    for (final int r : Arrays.asList(50, 100, 200)) { //ref count
                        final AllelicCount count = new AllelicCount(DUMMY, r, a);
                        final double actual = AlleleFractionLikelihoods.hetLogLikelihood(parameters, f, count, AlleleFractionIndicator.ALT_MINOR);
                        final double expected = a * log(beta) + Gamma.logGamma(alpha - a) - Gamma.logGamma(alpha)
                                + log((1 - pi) / 2) + a * log(f / (1 - f));
                        Assert.assertEquals(actual, expected, 1e-3);
                    }
                }
            }
        }
    }
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
@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);
                    }
                }
            }
        }
    }
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
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);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/**
 * 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);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/**
 * {@inheritDoc}
 *
 * The implementation of this method is based on:
 * <ul>
 * <li>
 * <a href="http://mathworld.wolfram.com/Chi-SquaredDistribution.html"> Chi-Squared Distribution</a>, equation
 * (9).</li>
 * <li>Casella, G., & Berger, R. (1990). <i>Statistical Inference</i>. Belmont, CA: Duxbury Press.</li>
 * </ul>
 */
public double cumulativeProbability(double x)
{
	double ret;

	if (x <= 0)
	{
		ret = 0;
	}
	else
	{
		ret = Gamma.regularizedGammaP(shape, x / scale);
	}

	return ret;
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
/**
 * Frequency Test with a Block
 * <p>
 * The focus of the test is the proportion of ones within M-bit
 * blocks. The purpose of this test is to determine whether the
 * frequency of ones in an M-bit block is approximately block_size/2,
 * as would be expected under an assumption of randomness. For block
 * size = 1, this test degenerates to test 1, frequency test for RNG.
 */
private void testBlockFrequency(final int[] epsilon, final int M) {
    final int n = epsilon.length;
    final int N = n / M;

    double sum = 0.0;

    for (int i = 0; i < N; i++) {
        int blockSum = 0;
        for (int j = 0; j < M; j++) {
            blockSum += epsilon[j + i * M];
        }
        double pi = (double) blockSum / (double) M;
        double v = pi - 0.5;
        sum += v * v;
    }
    double chi_squared = 4.0 * M * sum;
    double p_value = Gamma.regularizedGammaQ(N / 2, chi_squared / 2);

    assertTrue("RNG test failed, test block frequency.", p_value >= 0.01);
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
private double pr(int u, double eta) {
    int l;
    double sum, p;

    if (u == 0)
        p = Math.exp(-eta);
    else {
        sum = 0.0;
        for (l = 1; l <= u; l++)
            sum += Math.exp(-eta - u * Math.log(2)
                    + l * Math.log(eta) - Gamma.logGamma(l + 1)
                    + Gamma.logGamma(u) - Gamma.logGamma(l)
                    - Gamma.logGamma(u - l + 1));
        p = sum;
    }
    return p;
}
 

import org.apache.commons.math3.special.Gamma; //导入依赖的package包/类
public Map<String, Double> MStep(List<Double> expectation){
		Map<String, Double> parameters=new HashMap<String, Double>();

		double pi=EMUtils.getPi(expectation);
		parameters.put("pi", pi);

		double mu=EMUtils.getMu(sample, expectation);		
		double logmu=EMUtils.getLogMu(sample, expectation);
//		double sd=EMUtils.getSD(sample, expectation,mu);
//		double n=EMUtils.getN(expectation);
		double s=Math.log(mu)-logmu;

//		double scale=Math.pow(sd,2)/mu;
//		double shape=mu/scale;
		
//		double shape=1/Math.pow(sd/mu,2)-1d/n;
//		double scale=mu/shape;
		
		double shape=Double.POSITIVE_INFINITY;
		double shape2=(3-s+Math.pow(Math.pow(s-3,2)+24*s,0.5))/(12*s);
		int it=0;
		while(it<500&&Math.abs(shape-shape2)>1E-13){
			shape=shape2;
			shape2=shape-(Math.log(shape)-Gamma.digamma(shape)-s)/(1/shape-Gamma.trigamma(shape));
			it++;
		}
		shape=shape2;
		double scale=mu/shape;//wikipedia
		if(Double.isInfinite(shape)){
			System.out.println();
		}
		
		parameters.put("scale", scale);
		parameters.put("shape",shape);
		
		return parameters;
	}