Java Code Examples for org.apache.commons.math3.distribution.BinomialDistribution#probability()

double ploidyLikelihood(double ploidy, @NotNull final WeightedPloidy weighted) {
    final String binomialKey = weighted.alleleReadCount() + ":" + weighted.totalReadCount();
    final BinomialDistribution binomialDistribution = binomialDistributionMap.computeIfAbsent(binomialKey,
            s -> new BinomialDistribution(weighted.totalReadCount(), weighted.alleleFrequency()));

    double lowerBoundAlleleReadCount = Math.max(0, ploidy - modelWidth / 2d) / weighted.ploidy() * weighted.alleleReadCount();
    int lowerBoundAlleleReadCountRounded = (int) Math.round(lowerBoundAlleleReadCount);
    double lowerBoundAddition = lowerBoundAlleleReadCountRounded + 0.5 - lowerBoundAlleleReadCount;

    double upperBoundAlleleReadCount = Math.max(0, ploidy + modelWidth / 2d) / weighted.ploidy() * weighted.alleleReadCount();
    int upperBoundAlleleReadCountRounded = (int) Math.round(upperBoundAlleleReadCount);
    double upperBoundSubtraction = upperBoundAlleleReadCountRounded + 0.5 - upperBoundAlleleReadCount;

    double rawResult =
            binomialDistribution.cumulativeProbability(upperBoundAlleleReadCountRounded) - binomialDistribution.cumulativeProbability(
                    lowerBoundAlleleReadCountRounded) + lowerBoundAddition * binomialDistribution.probability(
                    lowerBoundAlleleReadCountRounded) - upperBoundSubtraction * binomialDistribution.probability(
                    upperBoundAlleleReadCountRounded);

    return Math.round(rawResult * 100) / 100d;
}
 

/**
 * Test sampling from a binomial distribution.
 */
@Test
public void testBinomialSamples() {
    final int trials = 67;
    final double probabilityOfSuccess = 0.345;
    final BinomialDistribution dist = new BinomialDistribution(null, trials, probabilityOfSuccess);
    final double[] expected = new double[trials + 1];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected, 1.0);
}
 

/**
 * Test sampling from a binomial distribution.
 */
@Test
public void testBinomialSamples() {
    final int trials = 67;
    final double probabilityOfSuccess = 0.345;
    final BinomialDistribution dist = new BinomialDistribution(trials, probabilityOfSuccess);
    final double[] expected = new double[trials + 1];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected);
}
 

/**
 * Returns the <i>observed significance level</i>, or
 * <a href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">p-value</a>,
 * associated with a <a href="http://en.wikipedia.org/wiki/Binomial_test"> Binomial test</a>.
 * <p>
 * The number returned is the smallest significance level at which one can reject the null hypothesis.
 * The form of the hypothesis depends on {@code alternativeHypothesis}.</p>
 * <p>
 * The p-Value represents the likelihood of getting a result at least as extreme as the sample,
 * given the provided {@code probability} of success on a single trial. For single-sided tests,
 * this value can be directly derived from the Binomial distribution. For the two-sided test,
 * the implementation works as follows: we start by looking at the most extreme cases
 * (0 success and n success where n is the number of trials from the sample) and determine their likelihood.
 * The lower value is added to the p-Value (if both values are equal, both are added). Then we continue with
 * the next extreme value, until we added the value for the actual observed sample.</p>
 * <p>
 * <strong>Preconditions</strong>:
 * <ul>
 * <li>Number of trials must be &ge; 0.</li>
 * <li>Number of successes must be &ge; 0.</li>
 * <li>Number of successes must be &le; number of trials.</li>
 * <li>Probability must be &ge; 0 and &le; 1.</li>
 * </ul></p>
 *
 * @param numberOfTrials number of trials performed
 * @param numberOfSuccesses number of successes observed
 * @param probability assumed probability of a single trial under the null hypothesis
 * @param alternativeHypothesis type of hypothesis being evaluated (one- or two-sided)
 * @return p-value
 * @throws NotPositiveException if {@code numberOfTrials} or {@code numberOfSuccesses} is negative
 * @throws OutOfRangeException if {@code probability} is not between 0 and 1
 * @throws MathIllegalArgumentException if {@code numberOfTrials} &lt; {@code numberOfSuccesses} or
 * if {@code alternateHypothesis} is null.
 * @see AlternativeHypothesis
 */
public double binomialTest(int numberOfTrials, int numberOfSuccesses, double probability,
                           AlternativeHypothesis alternativeHypothesis) {
    if (numberOfTrials < 0) {
        throw new NotPositiveException(numberOfTrials);
    }
    if (numberOfSuccesses < 0) {
        throw new NotPositiveException(numberOfSuccesses);
    }
    if (probability < 0 || probability > 1) {
        throw new OutOfRangeException(probability, 0, 1);
    }
    if (numberOfTrials < numberOfSuccesses) {
        throw new MathIllegalArgumentException(
            LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
            numberOfTrials, numberOfSuccesses);
    }
    if (alternativeHypothesis == null) {
        throw new NullArgumentException();
    }

    // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
    final BinomialDistribution distribution = new BinomialDistribution(null, numberOfTrials, probability);
    switch (alternativeHypothesis) {
    case GREATER_THAN:
        return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
    case LESS_THAN:
        return distribution.cumulativeProbability(numberOfSuccesses);
    case TWO_SIDED:
        int criticalValueLow = 0;
        int criticalValueHigh = numberOfTrials;
        double pTotal = 0;

        while (true) {
            double pLow = distribution.probability(criticalValueLow);
            double pHigh = distribution.probability(criticalValueHigh);

            if (pLow == pHigh) {
                pTotal += 2 * pLow;
                criticalValueLow++;
                criticalValueHigh--;
            } else if (pLow < pHigh) {
                pTotal += pLow;
                criticalValueLow++;
            } else {
                pTotal += pHigh;
                criticalValueHigh--;
            }

            if (criticalValueLow > numberOfSuccesses || criticalValueHigh < numberOfSuccesses) {
                break;
            }
        }
        return pTotal;
    default:
        throw new MathInternalError(LocalizedFormats. OUT_OF_RANGE_SIMPLE, alternativeHypothesis,
                  AlternativeHypothesis.TWO_SIDED, AlternativeHypothesis.LESS_THAN);
    }
}
 

/**
 * Returns the <i>observed significance level</i>, or
 * <a href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">p-value</a>,
 * associated with a <a href="http://en.wikipedia.org/wiki/Binomial_test"> Binomial test</a>.
 * <p>
 * The number returned is the smallest significance level at which one can reject the null hypothesis.
 * The form of the hypothesis depends on {@code alternativeHypothesis}.</p>
 * <p>
 * The p-Value represents the likelihood of getting a result at least as extreme as the sample,
 * given the provided {@code probability} of success on a single trial. For single-sided tests,
 * this value can be directly derived from the Binomial distribution. For the two-sided test,
 * the implementation works as follows: we start by looking at the most extreme cases
 * (0 success and n success where n is the number of trials from the sample) and determine their likelihood.
 * The lower value is added to the p-Value (if both values are equal, both are added). Then we continue with
 * the next extreme value, until we added the value for the actual observed sample.</p>
 * <p>
 * <strong>Preconditions</strong>:
 * <ul>
 * <li>Number of trials must be &ge; 0.</li>
 * <li>Number of successes must be &ge; 0.</li>
 * <li>Number of successes must be &le; number of trials.</li>
 * <li>Probability must be &ge; 0 and &le; 1.</li>
 * </ul></p>
 *
 * @param numberOfTrials number of trials performed
 * @param numberOfSuccesses number of successes observed
 * @param probability assumed probability of a single trial under the null hypothesis
 * @param alternativeHypothesis type of hypothesis being evaluated (one- or two-sided)
 * @return p-value
 * @throws NotPositiveException if {@code numberOfTrials} or {@code numberOfSuccesses} is negative
 * @throws OutOfRangeException if {@code probability} is not between 0 and 1
 * @throws MathIllegalArgumentException if {@code numberOfTrials} &lt; {@code numberOfSuccesses} or
 * if {@code alternateHypothesis} is null.
 * @see AlternativeHypothesis
 */
public double binomialTest(int numberOfTrials, int numberOfSuccesses, double probability,
                           AlternativeHypothesis alternativeHypothesis) {
    if (numberOfTrials < 0) {
        throw new NotPositiveException(numberOfTrials);
    }
    if (numberOfSuccesses < 0) {
        throw new NotPositiveException(numberOfSuccesses);
    }
    if (probability < 0 || probability > 1) {
        throw new OutOfRangeException(probability, 0, 1);
    }
    if (numberOfTrials < numberOfSuccesses) {
        throw new MathIllegalArgumentException(
            LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
            numberOfTrials, numberOfSuccesses);
    }
    if (alternativeHypothesis == null) {
        throw new NullArgumentException();
    }

    // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
    final BinomialDistribution distribution = new BinomialDistribution(null, numberOfTrials, probability);
    switch (alternativeHypothesis) {
    case GREATER_THAN:
        return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
    case LESS_THAN:
        return distribution.cumulativeProbability(numberOfSuccesses);
    case TWO_SIDED:
        int criticalValueLow = 0;
        int criticalValueHigh = numberOfTrials;
        double pTotal = 0;

        while (true) {
            double pLow = distribution.probability(criticalValueLow);
            double pHigh = distribution.probability(criticalValueHigh);

            if (pLow == pHigh) {
                pTotal += 2 * pLow;
                criticalValueLow++;
                criticalValueHigh--;
            } else if (pLow < pHigh) {
                pTotal += pLow;
                criticalValueLow++;
            } else {
                pTotal += pHigh;
                criticalValueHigh--;
            }

            if (criticalValueLow > numberOfSuccesses || criticalValueHigh < numberOfSuccesses) {
                break;
            }
        }
        return pTotal;
    default:
        throw new MathInternalError(LocalizedFormats. OUT_OF_RANGE_SIMPLE, alternativeHypothesis,
                  AlternativeHypothesis.TWO_SIDED, AlternativeHypothesis.LESS_THAN);
    }
}
 

/**
 * Returns the <i>observed significance level</i>, or
 * <a
 * href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">p-value</a>,
 * associated with a <a href="http://en.wikipedia.org/wiki/Binomial_test">
 * Binomial test</a>.
 * <p>
 * The number returned is the smallest significance level at which one can
 * reject the null hypothesis. The form of the hypothesis depends on
 * {@code alternativeHypothesis}.</p>
 * <p>
 * The p-Value represents the likelihood of getting a result at least as
 * extreme as the sample, given the provided {@code probability} of success
 * on a single trial. For single-sided tests, this value can be directly
 * derived from the Binomial distribution. For the two-sided test, the
 * implementation works as follows: we start by looking at the most extreme
 * cases (0 success and n success where n is the number of trials from the
 * sample) and determine their likelihood. The lower value is added to the
 * p-Value (if both values are equal, both are added). Then we continue with
 * the next extreme value, until we added the value for the actual observed
 * sample.</p>
 * <p>
 * <strong>Preconditions</strong>:
 * <ul>
 * <li>Number of trials must be &ge; 0.</li>
 * <li>Number of successes must be &ge; 0.</li>
 * <li>Number of successes must be &le; number of trials.</li>
 * <li>Probability must be &ge; 0 and &le; 1.</li>
 * </ul></p>
 *
 * @param numberOfTrials number of trials performed
 * @param numberOfSuccesses number of successes observed
 * @param probability assumed probability of a single trial under the null
 * hypothesis
 * @param alternativeHypothesis type of hypothesis being evaluated (one- or
 * two-sided)
 * @return p-value
 * @throws NotPositiveException if {@code numberOfTrials} or
 * {@code numberOfSuccesses} is negative
 * @throws OutOfRangeException if {@code probability} is not between 0 and 1
 * @throws MathIllegalArgumentException if {@code numberOfTrials} &lt;
 * {@code numberOfSuccesses} or if {@code alternateHypothesis} is null.
 * @see AlternativeHypothesis
 */
public double binomialTest(int numberOfTrials, int numberOfSuccesses, double probability,
		AlternativeHypothesis alternativeHypothesis) {
	if (numberOfTrials < 0) {
		throw new NotPositiveException(numberOfTrials);
	}
	if (numberOfSuccesses < 0) {
		throw new NotPositiveException(numberOfSuccesses);
	}
	if (probability < 0 || probability > 1) {
		throw new OutOfRangeException(probability, 0, 1);
	}
	if (numberOfTrials < numberOfSuccesses) {
		throw new MathIllegalArgumentException(
				LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
				numberOfTrials, numberOfSuccesses);
	}
	if (alternativeHypothesis == null) {
		throw new NullArgumentException();
	}

	final BinomialDistribution distribution = new BinomialDistribution(numberOfTrials, probability);
	switch (alternativeHypothesis) {
		case GREATER_THAN:
			return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
		case LESS_THAN:
			return distribution.cumulativeProbability(numberOfSuccesses);
		case TWO_SIDED:
			int criticalValueLow = 0;
			int criticalValueHigh = numberOfTrials;
			double pTotal = 0;

			while (true) {
				double pLow = distribution.probability(criticalValueLow);
				double pHigh = distribution.probability(criticalValueHigh);

				if (pLow == pHigh) {
					pTotal += 2 * pLow;
					criticalValueLow++;
					criticalValueHigh--;
				} else if (pLow < pHigh) {
					pTotal += pLow;
					criticalValueLow++;
				} else {
					pTotal += pHigh;
					criticalValueHigh--;
				}

				if (criticalValueLow > numberOfSuccesses || criticalValueHigh < numberOfSuccesses) {
					break;
				}
			}
			return pTotal;
		default:
			throw new MathInternalError(LocalizedFormats.OUT_OF_RANGE_SIMPLE, alternativeHypothesis,
					AlternativeHypothesis.TWO_SIDED, AlternativeHypothesis.LESS_THAN);
	}
}
 

/**
public static double getTwoSidedPvalueColt(int numTrials, int numSuccesses, double probability) {
	double result =0;
	//  R code to replicate.
	//twoTailedUnbiasedTest<-function (numSuccesses, numTrials) {
	//	x=c(numSuccesses, numTrials-numSuccesses)
	//	dbinom(min(x), sum(x), 0.5) + ifelse( min(x)==0, 0, 2*sum(dbinom(0:(min(x)-1), sum(x), 0.5)))
	//


	// make numSuccesses less than 1/2 of the num successes.
	numSuccesses=Math.min(numSuccesses, numTrials-numSuccesses);

	//result=2*sum(dbinom(0:(min(x)), sum(x), 0.5)
	Binomial b = new Binomial(numTrials, probability, rand);

	for (int i=0; i<numSuccesses; i++) {
		result+=b.pdf(i);
	}

	result=result*2;
	result+=b.pdf(numSuccesses);

	return (result);
}
*/

public static double getTwoSidedPvalue(final int numTrials, int numSuccesses, final double probability) {
	double result =0;

	/*  R code to replicate.
	twoTailedUnbiasedTest<-function (numSuccesses, numTrials) {
		x=c(numSuccesses, numTrials-numSuccesses)
		dbinom(min(x), sum(x), 0.5) + ifelse( min(x)==0, 0, 2*sum(dbinom(0:(min(x)-1), sum(x), 0.5)))
		#2*sum(dbinom(x=0:min(x[1], x[2]), size=x[1]+x[2], prob=0.5))
	}
	*/

	// make numSuccesses less than 1/2 of the num successes.
	numSuccesses=Math.min(numSuccesses, numTrials-numSuccesses);

	//result=2*sum(dbinom(0:(min(x)), sum(x), 0.5)
	BinomialDistribution bd = new BinomialDistribution(numTrials, probability);
	// double two = bd.cumulativeProbability(numberOfSuccesses);



	for (int i=0; i<numSuccesses; i++)
		result+=bd.probability(i);

	result=result*2;
	result+=bd.probability(numSuccesses);

	return (result);
}