Java Code Examples for org.apache.commons.math3.special.Gamma#logGamma()

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

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

    double sumOfU_i = this.naturalParameters.sum() ;

    double sumLogGammaOfU_i = this.getSparseNaturalParameters().sumApply(new Function<Double, Double>() {
        @Override
        public Double apply(Double aDouble) {
            return Gamma.logGamma(aDouble);
        }
    });

    return sumLogGammaOfU_i - Gamma.logGamma(sumOfU_i);

}
 

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

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

@Test
public void testFactorialLog() {
    // Cache size allows some of the factorials to be cached and some
    // to be under the precomputed factorials.
    FactorialLog factorialLog = FactorialLog.create().withCache(MAX_REPRESENTABLE / 2);
    Assert.assertEquals(0, factorialLog.value(0), 1e-10);
    for (int n = 1; n <= MAX_REPRESENTABLE + 5; n++) {
        // Use Commons math to compute logGamma(1 + n);
        double expected = Gamma.logGamma(1 + n);
        Assert.assertEquals(expected, factorialLog.value(n), 1e-10);
    }
}
 

/**
 *  this implements the isActive() algorithm described in docs/mutect/mutect.pdf
 *  the multiplicative factor is for the special case where we pass a singleton list
 *  of alt quals and want to duplicate that alt qual over multiple reads
 * @param nRef          ref read count
 * @param altQuals      Phred-scaled qualities of alt-supporting reads
 * @param repeatFactor  Number of times each alt qual is duplicated
 * @param afPrior       Beta prior on alt allele fraction
 * @return
 */
public static double logLikelihoodRatio(final int nRef, final List<Byte> altQuals, final int repeatFactor, final Optional<BetaDistributionShape> afPrior) {
    final int nAlt = repeatFactor * altQuals.size();
    final int n = nRef + nAlt;

    final double fTildeRatio = FastMath.exp(MathUtils.digamma(nRef + 1) - MathUtils.digamma(nAlt + 1));


    double readSum = 0;
    for (final byte qual : altQuals) {
        final double epsilon = QualityUtils.qualToErrorProb(qual);
        final double zBarAlt = (1 - epsilon) / (1 - epsilon + epsilon * fTildeRatio);
        final double logEpsilon = NaturalLogUtils.qualToLogErrorProb(qual);
        final double logOneMinusEpsilon = NaturalLogUtils.qualToLogProb(qual);
        readSum += zBarAlt * (logOneMinusEpsilon - logEpsilon) + MathUtils.fastBernoulliEntropy(zBarAlt);
    }

    final double betaEntropy;
    if (afPrior.isPresent()) {
        final double alpha = afPrior.get().getAlpha();
        final double beta = afPrior.get().getBeta();
        betaEntropy = Gamma.logGamma(alpha + beta) - Gamma.logGamma(alpha) - Gamma.logGamma(beta)
                - Gamma.logGamma(alpha + beta + n) + Gamma.logGamma(alpha + nAlt) + Gamma.logGamma(beta + nRef);
    } else {
        betaEntropy = MathUtils.log10ToLog(-MathUtils.log10Factorial(n + 1) + MathUtils.log10Factorial(nAlt) + MathUtils.log10Factorial(nRef));
    }
    return betaEntropy + readSum * repeatFactor;
}
 

/**
 * {@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}
 */
@Override
public double computeLogNormalizer() {
    double alpha = this.naturalParameters.get(EF_JointNormalGamma.LOGGAMMA) + 0.5 ;
    double precision =  this.naturalParameters.get(EF_JointNormalGamma.MUSQUARE_GAMMA)/(-0.5);
    double mean =  this.naturalParameters.get(EF_JointNormalGamma.MU_GAMMA)/precision;
    double beta = -this.naturalParameters.get(EF_JointNormalGamma.GAMMA) -0.5*precision*mean*mean;

    return -0.5*Math.log(precision) - alpha*Math.log(beta) + Gamma.logGamma(alpha);
}
 

private static double betabinPMFG(int n, int k, double a, double b) {
    double b1 = Gamma.logGamma(n + 1);
    double b2 = Gamma.logGamma(k + 1) + Gamma.logGamma(n - k + 1);
    double b3 = Gamma.logGamma(k + a) + Gamma.logGamma(n - k + b);
    double b4 = Gamma.logGamma(n + a + b);
    double b5 = Gamma.logGamma(a + b);
    double b6 = Gamma.logGamma(a) + Gamma.logGamma(b);
    double v = b1 - b2 + b3 - b4 + b5 - b6;
    return Math.exp(v);
}
 

/**
 * Creates an instance.
 * 
 * @param degrees The number of degrees of freedom, not negative or zero
 * @param nonCentrality The non-centrality parameter, not negative
 */
public NonCentralChiSquaredDistribution(double degrees, double nonCentrality) {
  ArgChecker.isTrue(degrees > 0, "degrees of freedom must be > 0, have " + degrees);
  ArgChecker.isTrue(nonCentrality >= 0, "non-centrality must be >= 0, have " + nonCentrality);
  _dofOverTwo = degrees / 2.0;
  _lambdaOverTwo = nonCentrality / 2.0;
  _k = (int) Math.round(_lambdaOverTwo);

  if (_lambdaOverTwo == 0) {
    _pStart = 0.0;
  } else {
    double logP = -_lambdaOverTwo + _k * Math.log(_lambdaOverTwo) - Gamma.logGamma(_k + 1);
    _pStart = Math.exp(logP);
  }
}
 

/**
 * {@inheritDoc}
 */
@Override
public double getCDF(Double x) {
  ArgChecker.notNull(x, "x");
  if (x < 0) {
    return 0.0;
  }

  if ((_dofOverTwo + _lambdaOverTwo) > 1000) {
    return getFraserApproxCDF(x);
  }

  double regGammaStart = 0;
  double halfX = x / 2.0;
  double logX = Math.log(halfX);
  try {
    regGammaStart = Gamma.regularizedGammaP(_dofOverTwo + _k, halfX);
  } catch (MaxCountExceededException ex) {
    throw new MathException(ex);
  }

  double sum = _pStart * regGammaStart;
  double oldSum = Double.NEGATIVE_INFINITY;
  double p = _pStart;
  double regGamma = regGammaStart;
  double temp;
  int i = _k;

  // first add terms below _k
  while (i > 0 && Math.abs(sum - oldSum) / sum > _eps) {
    i--;
    p *= (i + 1) / _lambdaOverTwo;
    temp = (_dofOverTwo + i) * logX - halfX - Gamma.logGamma(_dofOverTwo + i + 1);
    regGamma += Math.exp(temp);
    oldSum = sum;
    sum += p * regGamma;
  }

  p = _pStart;
  regGamma = regGammaStart;
  oldSum = Double.NEGATIVE_INFINITY;
  i = _k;
  while (Math.abs(sum - oldSum) / sum > _eps) {
    i++;
    p *= _lambdaOverTwo / i;
    temp = (_dofOverTwo + i - 1) * logX - halfX - Gamma.logGamma(_dofOverTwo + i);
    regGamma -= Math.exp(temp);
    oldSum = sum;
    sum += p * regGamma;
  }

  return sum;
}
 

/**
 * Optimize beta using Minka's Fixed Point Iteration
 */
public void optimizeBetaMinka() {

	// Optimize beta
	final int nTopics = Topic.nTopics;
	final int nBackTopics = Topic.nBackTopics;
	final double[] residual = new double[nTopics];
	final double W = (double) nTokensCorpus;
	do {
		// Get necessary sums
		final double[] topSum = new double[nTopics];
		final double[] bottomSum = new double[nTopics];

		// minus parts
		for (int k = 0; k < nTopics; k++) {
			topSum[k] -= W * digamma(beta[k]);
			bottomSum[k] -= digamma(W * beta[k]);
		}

		// Get summed topic counts
		@SuppressWarnings("unchecked")
		final HashMultiset<Integer>[] sums = (HashMultiset<Integer>[]) new HashMultiset[nTopics];
		for (int k = 0; k < nTopics; k++)
			sums[k] = HashMultiset.create();

		for (int b = 0; b < nBackTopics; b++)
			sums[Topic.BACKGROUND[b]].addAll(btopic[b].getTokenMultiSet());
		for (int ci = 0; ci < nclusters; ci++) {
			sums[Topic.CONTENT].addAll(ctopic[ci].getTokenMultiSet());
			for (int di = 0; di < corpus.getCluster(ci).ndocs(); di++) {
				sums[Topic.DOCUMENT].addAll(dtopic[ci].get(di)
						.getTokenMultiSet());
				// for (int si = 0; si < corpus.getCluster(ci).getDoc(di)
				// .nsents(); si++)
				// sums[Topic.SENTENCE].addAll(stopic[ci].get(di)[si]
				// .getTokenMultiSet());
			}
		}

		// topics top sum and total summed topic counts
		final double[] sumTotals = new double[nTopics];
		for (int k = 0; k < nTopics; k++) {

			for (final Integer ti : sums[k].elementSet()) {
				topSum[k] += Gamma.logGamma(sums[k].count(ti) + beta[k]);
				sumTotals[k] += sums[k].count(ti);
			}
		}

		// topics bottom sum
		for (int k = 0; k < nTopics; k++)
			bottomSum[k] += Gamma.logGamma(sumTotals[k] + W * beta[k]);

		// Estimate beta_k
		for (int k = 0; k < nTopics; k++) {
			residual[k] = beta[k];
			beta[k] = (beta[k] * topSum[k]) / (W * bottomSum[k]);
			residual[k] -= beta[k];
		}
	} while (StatsUtil.norm(residual) > HYPER_BETA_TOL);
}
 

/** Recompute the normalization factor. */
private void recomputeZ() {
    if (Double.isNaN(z)) {
        z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
    }
}
 

/** Recompute the normalization factor. */
private void recomputeZ() {
    if (Double.isNaN(z)) {
        z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
    }
}
 

public static double logDirichletNormalization(final double... dirichletParams) {
    final double logNumerator = Gamma.logGamma(MathUtils.sum(dirichletParams));
    final double logDenominator = MathUtils.sum(MathUtils.applyToArray(dirichletParams, Gamma::logGamma));
    return logNumerator - logDenominator;
}
 

/** Recompute the normalization factor. */
private void recomputeZ() {
    if (Double.isNaN(z)) {
        z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
    }
}
 

public static double[] computeLocalKLDirichlet(EF_Dirichlet dirichletQ, EF_Dirichlet dirichletP) {

        double[] localKL = new double[dirichletP.sizeOfSufficientStatistics()];

        int nstates = localKL.length;
        double sumCountsQ = dirichletQ.getNaturalParameters().sum();
        double sumCountsP = dirichletP.getNaturalParameters().sum();

        for (int i = 0; i < localKL.length; i++) {
            double alphaQ = dirichletQ.getNaturalParameters().get(i);
            double alphaP = dirichletP.getNaturalParameters().get(i);

            double kl =  alphaQ - alphaP;
            kl*=dirichletQ.getMomentParameters().get(i);

            kl-=(Gamma.logGamma(alphaQ) - Gamma.logGamma(sumCountsQ)/nstates);
            kl+=(Gamma.logGamma(alphaP) - Gamma.logGamma(sumCountsP)/nstates);

            localKL[i] = kl;
        }


        return localKL;
    }
 

/**
 * Calculate lgamma of the input value.
 */
static double lgamma(double x) {
    return Gamma.logGamma(x);
}
 

/** Recompute the normalization factor. */
private void recomputeZ() {
    if (Double.isNaN(z)) {
        z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
    }
}
 

/**
 * Calculates the log of the integral of the allelic-bias posterior (approximated as a Gamma distribution
 * with mode and curvature given by {@link AlleleFractionLikelihoods#biasPosteriorMode}
 * and {@link AlleleFractionLikelihoods#biasPosteriorCurvature}, respectively)
 * at given values of the hyperparameters for the allelic-bias Gamma-distribution prior,
 * the minor-allele fraction parameter, and the observed counts at a site.
 * See docs/CNVs/CNV-methods.pdf (where this quantity is referred to as phi) for details.
 * @param alpha alpha hyperparameter for allelic-bias Gamma-distribution prior
 * @param beta  beta hyperparameter for allelic-bias Gamma-distribution prior
 * @param f     minor-allele fraction
 * @param a     alt read count
 * @param r     ref read count
 */
protected static double logIntegralOverAllelicBias(final double alpha, final double beta, final double f, final int a, final int r) {
    final double lambda0 = biasPosteriorMode(alpha, beta, f, a, r);
    final int n = a + r;
    final double kappa = biasPosteriorCurvature(alpha, f, r, n, lambda0);
    final double rho = biasPosteriorEffectiveAlpha(lambda0, kappa);
    final double tau = biasPosteriorEffectiveBeta(lambda0, kappa);
    final double logc = alpha*log(beta) - Gamma.logGamma(alpha) + a * log(f) + r * log(1 - f)
            + (r + alpha - rho) * log(lambda0) + (tau - beta) * lambda0 - n * log(f + (1 - f) * lambda0);
    return logc + Gamma.logGamma(rho) - rho * log(tau);
}