001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.stat.inference;
018
019import org.apache.commons.statistics.distribution.NormalDistribution;
020import org.apache.commons.math4.legacy.exception.ConvergenceException;
021import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
022import org.apache.commons.math4.legacy.exception.NoDataException;
023import org.apache.commons.math4.legacy.exception.NullArgumentException;
024import org.apache.commons.math4.legacy.stat.ranking.NaNStrategy;
025import org.apache.commons.math4.legacy.stat.ranking.NaturalRanking;
026import org.apache.commons.math4.legacy.stat.ranking.TiesStrategy;
027import org.apache.commons.math4.core.jdkmath.JdkMath;
028
029import java.util.stream.IntStream;
030
031/**
032 * An implementation of the Mann-Whitney U test (also called Wilcoxon rank-sum test).
033 *
034 */
035public class MannWhitneyUTest {
036
037    /** Ranking algorithm. */
038    private NaturalRanking naturalRanking;
039
040    /**
041     * Create a test instance using where NaN's are left in place and ties get
042     * the average of applicable ranks. Use this unless you are very sure of
043     * what you are doing.
044     */
045    public MannWhitneyUTest() {
046        naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
047                TiesStrategy.AVERAGE);
048    }
049
050    /**
051     * Create a test instance using the given strategies for NaN's and ties.
052     * Only use this if you are sure of what you are doing.
053     *
054     * @param nanStrategy
055     *            specifies the strategy that should be used for Double.NaN's
056     * @param tiesStrategy
057     *            specifies the strategy that should be used for ties
058     */
059    public MannWhitneyUTest(final NaNStrategy nanStrategy,
060                            final TiesStrategy tiesStrategy) {
061        naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
062    }
063
064    /**
065     * Ensures that the provided arrays fulfills the assumptions.
066     *
067     * @param x first sample
068     * @param y second sample
069     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
070     * @throws NoDataException if {@code x} or {@code y} are zero-length.
071     */
072    private void ensureDataConformance(final double[] x, final double[] y)
073        throws NullArgumentException, NoDataException {
074
075        if (x == null ||
076            y == null) {
077            throw new NullArgumentException();
078        }
079        if (x.length == 0 ||
080            y.length == 0) {
081            throw new NoDataException();
082        }
083    }
084
085    /** Concatenate the samples into one array.
086     * @param x first sample
087     * @param y second sample
088     * @return concatenated array
089     */
090    private double[] concatenateSamples(final double[] x, final double[] y) {
091        final double[] z = new double[x.length + y.length];
092
093        System.arraycopy(x, 0, z, 0, x.length);
094        System.arraycopy(y, 0, z, x.length, y.length);
095
096        return z;
097    }
098
099    /**
100     * Computes the <a
101     * href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney
102     * U statistic</a> comparing mean for two independent samples possibly of
103     * different length.
104     * <p>
105     * This statistic can be used to perform a Mann-Whitney U test evaluating
106     * the null hypothesis that the two independent samples has equal mean.
107     * </p>
108     * <p>
109     * Let X<sub>i</sub> denote the i'th individual of the first sample and
110     * Y<sub>j</sub> the j'th individual in the second sample. Note that the
111     * samples would often have different length.
112     * </p>
113     * <p>
114     * <strong>Preconditions</strong>:
115     * <ul>
116     * <li>All observations in the two samples are independent.</li>
117     * <li>The observations are at least ordinal (continuous are also ordinal).</li>
118     * </ul>
119     *
120     * @param x the first sample
121     * @param y the second sample
122     * @return Mann-Whitney U statistic (minimum of U<sup>x</sup> and U<sup>y</sup>)
123     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
124     * @throws NoDataException if {@code x} or {@code y} are zero-length.
125     */
126    public double mannWhitneyU(final double[] x, final double[] y)
127        throws NullArgumentException, NoDataException {
128
129        ensureDataConformance(x, y);
130
131        final double[] z = concatenateSamples(x, y);
132        final double[] ranks = naturalRanking.rank(z);
133
134        double sumRankX = 0;
135
136        /*
137         * The ranks for x is in the first x.length entries in ranks because x
138         * is in the first x.length entries in z
139         */
140        sumRankX = IntStream.range(0, x.length).mapToDouble(i -> ranks[i]).sum();
141
142        /*
143         * U1 = R1 - (n1 * (n1 + 1)) / 2 where R1 is sum of ranks for sample 1,
144         * e.g. x, n1 is the number of observations in sample 1.
145         */
146        final double u1 = sumRankX - ((long) x.length * (x.length + 1)) / 2;
147
148        /*
149         * It can be shown that U1 + U2 = n1 * n2
150         */
151        final double u2 = (long) x.length * y.length - u1;
152
153        return JdkMath.min(u1, u2);
154    }
155
156    /**
157     * @param umin smallest Mann-Whitney U value
158     * @param n1 number of subjects in first sample
159     * @param n2 number of subjects in second sample
160     * @return two-sided asymptotic p-value
161     * @throws ConvergenceException if the p-value can not be computed
162     * due to a convergence error
163     * @throws MaxCountExceededException if the maximum number of
164     * iterations is exceeded
165     */
166    private double calculateAsymptoticPValue(final double umin,
167                                             final int n1,
168                                             final int n2)
169        throws ConvergenceException, MaxCountExceededException {
170
171        /* long multiplication to avoid overflow (double not used due to efficiency
172         * and to avoid precision loss)
173         */
174        final long n1n2prod = (long) n1 * n2;
175
176        // http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U#Normal_approximation
177        final double eU = n1n2prod / 2.0;
178        final double varU = n1n2prod * (n1 + n2 + 1) / 12.0;
179
180        final double z = (umin - eU) / JdkMath.sqrt(varU);
181
182        // No try-catch or advertised exception because args are valid
183        // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
184        final NormalDistribution standardNormal = NormalDistribution.of(0, 1);
185
186        return 2 * standardNormal.cumulativeProbability(z);
187    }
188
189    /**
190     * Returns the asymptotic <i>observed significance level</i>, or <a href=
191     * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
192     * p-value</a>, associated with a <a
193     * href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney
194     * U statistic</a> comparing mean for two independent samples.
195     * <p>
196     * Let X<sub>i</sub> denote the i'th individual of the first sample and
197     * Y<sub>j</sub> the j'th individual in the second sample. Note that the
198     * samples would often have different length.
199     * </p>
200     * <p>
201     * <strong>Preconditions</strong>:
202     * <ul>
203     * <li>All observations in the two samples are independent.</li>
204     * <li>The observations are at least ordinal (continuous are also ordinal).</li>
205     * </ul><p>
206     * Ties give rise to biased variance at the moment. See e.g. <a
207     * href="http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf"
208     * >http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf</a>.</p>
209     *
210     * @param x the first sample
211     * @param y the second sample
212     * @return asymptotic p-value
213     * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
214     * @throws NoDataException if {@code x} or {@code y} are zero-length.
215     * @throws ConvergenceException if the p-value can not be computed due to a
216     * convergence error
217     * @throws MaxCountExceededException if the maximum number of iterations
218     * is exceeded
219     */
220    public double mannWhitneyUTest(final double[] x, final double[] y)
221        throws NullArgumentException, NoDataException,
222        ConvergenceException, MaxCountExceededException {
223
224        ensureDataConformance(x, y);
225
226        final double uMin = mannWhitneyU(x, y);
227
228        return calculateAsymptoticPValue(uMin, x.length, y.length);
229    }
230}