MannWhitneyUTest.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/*
* This is not the original file distributed by the Apache Software Foundation
* It has been modified by the Hipparchus project
*/
package org.hipparchus.stat.inference;
import java.util.Map;
import java.util.TreeMap;
import java.util.stream.LongStream;
import org.hipparchus.distribution.continuous.NormalDistribution;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.exception.NullArgumentException;
import org.hipparchus.stat.LocalizedStatFormats;
import org.hipparchus.stat.ranking.NaNStrategy;
import org.hipparchus.stat.ranking.NaturalRanking;
import org.hipparchus.stat.ranking.TiesStrategy;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
/**
* An implementation of the Mann-Whitney U test.
* <p>
* The definitions and computing formulas used in this implementation follow
* those in the article,
* <a href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney U
* Test</a>
* <p>
* In general, results correspond to (and have been tested against) the R
* wilcox.test function, with {@code exact} meaning the same thing in both APIs
* and {@code CORRECT} uniformly true in this implementation. For example,
* wilcox.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE
* correct = TRUE) will return the same p-value as mannWhitneyUTest(x, y,
* false). The minimum of the W value returned by R for wilcox.test(x, y...) and
* wilcox.test(y, x...) should equal mannWhitneyU(x, y...).
*/
public class MannWhitneyUTest { // NOPMD - this is not a Junit test class, PMD false positive here
/**
* If the combined dataset contains no more values than this, test defaults to
* exact test.
*/
private static final int SMALL_SAMPLE_SIZE = 50;
/** Ranking algorithm. */
private final NaturalRanking naturalRanking;
/** Normal distribution */
private final NormalDistribution standardNormal;
/**
* Create a test instance using where NaN's are left in place and ties get
* the average of applicable ranks.
*/
public MannWhitneyUTest() {
naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
TiesStrategy.AVERAGE);
standardNormal = new NormalDistribution(0, 1);
}
/**
* Create a test instance using the given strategies for NaN's and ties.
*
* @param nanStrategy specifies the strategy that should be used for
* Double.NaN's
* @param tiesStrategy specifies the strategy that should be used for ties
*/
public MannWhitneyUTest(final NaNStrategy nanStrategy,
final TiesStrategy tiesStrategy) {
naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
standardNormal = new NormalDistribution(0, 1);
}
/**
* Computes the
* <a href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U">
* Mann-Whitney U statistic</a> comparing means for two independent samples
* possibly of different lengths.
* <p>
* This statistic can be used to perform a Mann-Whitney U test evaluating
* the null hypothesis that the two independent samples have equal mean.
* <p>
* Let X<sub>i</sub> denote the i'th individual of the first sample and
* Y<sub>j</sub> the j'th individual in the second sample. Note that the
* samples can have different lengths.
* <p>
* <strong>Preconditions</strong>:
* <ul>
* <li>All observations in the two samples are independent.</li>
* <li>The observations are at least ordinal (continuous are also
* ordinal).</li>
* </ul>
*
* @param x the first sample
* @param y the second sample
* @return Mann-Whitney U statistic (minimum of U<sup>x</sup> and
* U<sup>y</sup>)
* @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
* @throws MathIllegalArgumentException if {@code x} or {@code y} are
* zero-length.
*/
public double mannWhitneyU(final double[] x, final double[] y)
throws MathIllegalArgumentException, NullArgumentException {
ensureDataConformance(x, y);
final double[] z = concatenateSamples(x, y);
final double[] ranks = naturalRanking.rank(z);
double sumRankX = 0;
/*
* The ranks for x is in the first x.length entries in ranks because x
* is in the first x.length entries in z
*/
for (int i = 0; i < x.length; ++i) {
sumRankX += ranks[i];
}
/*
* U1 = R1 - (n1 * (n1 + 1)) / 2 where R1 is sum of ranks for sample 1,
* e.g. x, n1 is the number of observations in sample 1.
*/
final double U1 = sumRankX - ((long) x.length * (x.length + 1)) / 2;
/*
* U1 + U2 = n1 * n2
*/
final double U2 = (long) x.length * y.length - U1;
return FastMath.min(U1, U2);
}
/**
* Concatenate the samples into one array.
*
* @param x first sample
* @param y second sample
* @return concatenated array
*/
private double[] concatenateSamples(final double[] x, final double[] y) {
final double[] z = new double[x.length + y.length];
System.arraycopy(x, 0, z, 0, x.length);
System.arraycopy(y, 0, z, x.length, y.length);
return z;
}
/**
* Returns the asymptotic <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/Mann%E2%80%93Whitney_U">Mann-Whitney U
* Test</a> comparing means for two independent samples.
* <p>
* Let X<sub>i</sub> denote the i'th individual of the first sample and
* Y<sub>j</sub> the j'th individual in the second sample.
* <p>
* <strong>Preconditions</strong>:
* <ul>
* <li>All observations in the two samples are independent.</li>
* <li>The observations are at least ordinal.</li>
* </ul>
* <p>
* If there are no ties in the data and both samples are small (less than or
* equal to 50 values in the combined dataset), an exact test is performed;
* otherwise the test uses the normal approximation (with continuity
* correction).
* <p>
* If the combined dataset contains ties, the variance used in the normal
* approximation is bias-adjusted using the formula in the reference above.
*
* @param x the first sample
* @param y the second sample
* @return approximate 2-sized p-value
* @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
* @throws MathIllegalArgumentException if {@code x} or {@code y} are
* zero-length
*/
public double mannWhitneyUTest(final double[] x, final double[] y)
throws MathIllegalArgumentException, NullArgumentException {
ensureDataConformance(x, y);
// If samples are both small and there are no ties, perform exact test
if (x.length + y.length <= SMALL_SAMPLE_SIZE &&
tiesMap(x, y).isEmpty()) {
return mannWhitneyUTest(x, y, true);
} else { // Normal approximation
return mannWhitneyUTest(x, y, false);
}
}
/**
* Returns the asymptotic <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/Mann%E2%80%93Whitney_U">Mann-Whitney U
* Test</a> comparing means for two independent samples.
* <p>
* Let X<sub>i</sub> denote the i'th individual of the first sample and
* Y<sub>j</sub> the j'th individual in the second sample.
* <p>
* <strong>Preconditions</strong>:
* <ul>
* <li>All observations in the two samples are independent.</li>
* <li>The observations are at least ordinal.</li>
* </ul>
* <p>
* If {@code exact} is {@code true}, the p-value reported is exact, computed
* using the exact distribution of the U statistic. The computation in this
* case requires storage on the order of the product of the two sample
* sizes, so this should not be used for large samples.
* <p>
* If {@code exact} is {@code false}, the normal approximation is used to
* estimate the p-value.
* <p>
* If the combined dataset contains ties and {@code exact} is {@code true},
* MathIllegalArgumentException is thrown. If {@code exact} is {@code false}
* and the ties are present, the variance used to compute the approximate
* p-value in the normal approximation is bias-adjusted using the formula in
* the reference above.
*
* @param x the first sample
* @param y the second sample
* @param exact true means compute the p-value exactly, false means use the
* normal approximation
* @return approximate 2-sided p-value
* @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
* @throws MathIllegalArgumentException if {@code x} or {@code y} are
* zero-length or if {@code exact} is {@code true} and ties are
* present in the data
*/
public double mannWhitneyUTest(final double[] x, final double[] y,
final boolean exact)
throws MathIllegalArgumentException, NullArgumentException {
ensureDataConformance(x, y);
final Map<Double, Integer> tiesMap = tiesMap(x, y);
final double u = mannWhitneyU(x, y);
if (exact) {
if (!tiesMap.isEmpty()) {
throw new MathIllegalArgumentException(LocalizedStatFormats.TIES_ARE_NOT_ALLOWED);
}
return exactP(x.length, y.length, u);
}
return approximateP(u, x.length, y.length,
varU(x.length, y.length, tiesMap));
}
/**
* Ensures that the provided arrays fulfills the assumptions.
*
* @param x first sample
* @param y second sample
* @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
* @throws MathIllegalArgumentException if {@code x} or {@code y} are
* zero-length.
*/
private void ensureDataConformance(final double[] x, final double[] y)
throws MathIllegalArgumentException, NullArgumentException {
if (x == null || y == null) {
throw new NullArgumentException();
}
if (x.length == 0 || y.length == 0) {
throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
}
}
/**
* Estimates the 2-sided p-value associated with a Mann-Whitney U statistic
* value using the normal approximation.
* <p>
* The variance passed in is assumed to be corrected for ties. Continuity
* correction is applied to the normal approximation.
*
* @param u Mann-Whitney U statistic
* @param n1 number of subjects in first sample
* @param n2 number of subjects in second sample
* @param varU variance of U (corrected for ties if these exist)
* @return two-sided asymptotic p-value
* @throws MathIllegalStateException if the p-value can not be computed due
* to a convergence error
* @throws MathIllegalStateException if the maximum number of iterations is
* exceeded
*/
private double approximateP(final double u, final int n1, final int n2,
final double varU)
throws MathIllegalStateException {
final double mu = (long) n1 * n2 / 2.0;
// If u == mu, return 1
if (Precision.equals(mu, u)) {
return 1;
}
// Force z <= 0 so we get tail probability. Also apply continuity
// correction
final double z = -Math.abs((u - mu) + 0.5) / FastMath.sqrt(varU);
return 2 * standardNormal.cumulativeProbability(z);
}
/**
* Calculates the (2-sided) p-value associated with a Mann-Whitney U
* statistic.
* <p>
* To compute the p-value, the probability densities for each value of U up
* to and including u are summed and the resulting tail probability is
* multiplied by 2.
* <p>
* The result of this computation is only valid when the combined n + m
* sample has no tied values.
* <p>
* This method should not be used for large values of n or m as it maintains
* work arrays of size n*m.
*
* @param u Mann-Whitney U statistic value
* @param n first sample size
* @param m second sample size
* @return two-sided exact p-value
*/
private double exactP(final int n, final int m, final double u) {
final double nm = m * n;
if (u > nm) { // Quick exit if u is out of range
return 1;
}
// Need to convert u to a mean deviation, so cumulative probability is
// tail probability
final double crit = u < nm / 2 ? u : nm / 2 - u;
double cum = 0d;
for (int ct = 0; ct <= crit; ct++) {
cum += uDensity(n, m, ct);
}
return 2 * cum;
}
/**
* Computes the probability density function for the Mann-Whitney U
* statistic.
* <p>
* This method should not be used for large values of n or m as it maintains
* work arrays of size n*m.
*
* @param n first sample size
* @param m second sample size
* @param u U-statistic value
* @return the probability that a U statistic derived from random samples of
* size n and m (containing no ties) equals u
*/
private double uDensity(final int n, final int m, double u) {
if (u < 0 || u > m * n) {
return 0;
}
final long[] freq = uFrequencies(n, m);
return freq[(int) FastMath.round(u + 1)] /
(double) LongStream.of(freq).sum();
}
/**
* Computes frequency counts for values of the Mann-Whitney U statistc. If
* freq[] is the returned array, freq[u + 1] counts the frequency of U = u
* among all possible n-m orderings. Therefore, P(u = U) = freq[u + 1] / sum
* where sum is the sum of the values in the returned array.
* <p>
* Implements the algorithm presented in "Algorithm AS 62: A Generator for
* the Sampling Distribution of the Mann-Whitney U Statistic", L. C. Dinneen
* and B. C. Blakesley Journal of the Royal Statistical Society. Series C
* (Applied Statistics) Vol. 22, No. 2 (1973), pp. 269-273.
*
* @param n first sample size
* @param m second sample size
* @return array of U statistic value frequencies
*/
private long[] uFrequencies(final int n, final int m) {
final int max = FastMath.max(m, n);
if (max > 100) {
throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
max, 100);
}
final int min = FastMath.min(m, n);
final long[] out = new long[n * m + 2];
final long[] work = new long[n * m + 2];
for (int i = 1; i < out.length; i++) {
out[i] = (i <= (max + 1)) ? 1 : 0;
}
work[1] = 0;
int in = max;
for (int i = 2; i <= min; i++) {
work[i] = 0;
in = in + max;
int n1 = in + 2;
long l = 1 + in / 2;
int k = i;
for (int j = 1; j <= l; j++) {
k++;
n1 = n1 - 1;
final long sum = out[j] + work[j];
out[j] = sum;
work[k] = sum - out[n1];
out[n1] = sum;
}
}
return out;
}
/**
* Computes the variance for a U-statistic associated with samples of
* sizes{@code n} and {@code m} and ties described by {@code tiesMap}. If
* {@code tiesMap} is non-empty, the multiplicity counts in its values set
* are used to adjust the variance.
*
* @param n first sample size
* @param m second sample size
* @param tiesMap map of <value, multiplicity>
* @return ties-adjusted variance
*/
private double varU(final int n, final int m,
Map<Double, Integer> tiesMap) {
final double nm = (long) n * m;
if (tiesMap.isEmpty()) {
return nm * (n + m + 1) / 12.0;
}
final long tSum = tiesMap.entrySet().stream()
.mapToLong(e -> e.getValue() * e.getValue() * e.getValue() -
e.getValue())
.sum();
final double totalN = n + m;
return (nm / 12) * (totalN + 1 - tSum / (totalN * (totalN - 1)));
}
/**
* Creates a map whose keys are values occurring more than once in the
* combined dataset formed from x and y. Map entry values are the number of
* occurrences. The returned map is empty iff there are no ties in the data.
*
* @param x first dataset
* @param y second dataset
* @return map of <value, number of times it occurs> for values occurring
* more than once or an empty map if there are no ties (the returned
* map is <em>not</em> thread-safe, which is OK in the context of the callers)
*/
private Map<Double, Integer> tiesMap(final double[] x, final double[] y) {
final Map<Double, Integer> tiesMap = new TreeMap<>(); // NOPMD - no concurrent access in the callers context
for (int i = 0; i < x.length; i++) {
tiesMap.merge(x[i], 1, Integer::sum);
}
for (int i = 0; i < y.length; i++) {
tiesMap.merge(y[i], 1, Integer::sum);
}
tiesMap.entrySet().removeIf(e -> e.getValue() == 1);
return tiesMap;
}
}