in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/FisherExactTest.java [171:213]
private static double twoSidedTest(int k, HypergeometricDistribution distribution) {
// Find all i where Pr(X = i) <= Pr(X = k) and sum them.
// Exploit the known unimodal distribution to increase the
// search speed. Note the search depends only on magnitude differences.
// The current HypergeometricDistribution is faster using log probability
// as it omits a call to Math.exp.
// Use the mode as the point of largest probability.
// The lower or upper mode is important for the search below.
final int nn = distribution.getPopulationSize();
final int kk = distribution.getNumberOfSuccesses();
final int n = distribution.getSampleSize();
final double v = ((double) n + 1) * ((double) kk + 1) / (nn + 2.0);
final int m1 = (int) Math.ceil(v) - 1;
final int m2 = (int) Math.floor(v);
if (k < m1) {
final double pk = distribution.logProbability(k);
// Lower half = cdf(k)
// Find upper half. As k < lower mode i should never
// reach the lower mode based on the probability alone.
// Bracket with the upper mode.
final int i = Searches.searchDescending(m2, distribution.getSupportUpperBound(), pk,
distribution::logProbability);
return distribution.cumulativeProbability(k) +
distribution.survivalProbability(i - 1);
} else if (k > m2) {
final double pk = distribution.logProbability(k);
// Upper half = sf(k - 1)
// Find lower half. As k > upper mode i should never
// reach the upper mode based on the probability alone.
// Bracket with the lower mode.
final int i = Searches.searchAscending(distribution.getSupportLowerBound(), m1, pk,
distribution::logProbability);
return distribution.cumulativeProbability(i) +
distribution.survivalProbability(k - 1);
}
// k == mode
// Edge case where the sum of probabilities will be either
// 1 or 1 - Pr(X = mode) where mode != k
final double pk = distribution.probability(k);
final double pm = distribution.probability(k == m1 ? m2 : m1);
return pm > pk ? 1 - pm : 1;
}