    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     */
    017    package org.apache.commons.math3.stat.inference;
    018
    019    import org.apache.commons.math3.distribution.NormalDistribution;
    020    import org.apache.commons.math3.exception.ConvergenceException;
    021    import org.apache.commons.math3.exception.DimensionMismatchException;
    022    import org.apache.commons.math3.exception.MaxCountExceededException;
    023    import org.apache.commons.math3.exception.NoDataException;
    024    import org.apache.commons.math3.exception.NullArgumentException;
    025    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    026    import org.apache.commons.math3.stat.ranking.NaNStrategy;
    027    import org.apache.commons.math3.stat.ranking.NaturalRanking;
    028    import org.apache.commons.math3.stat.ranking.TiesStrategy;
    029    import org.apache.commons.math3.util.FastMath;
    030
    031    /**
    032     * An implementation of the Wilcoxon signed-rank test.
    033     *
    034     * @version $Id: WilcoxonSignedRankTest.java 1416643 2012-12-03 19:37:14Z tn $
    035     */
    036    public class WilcoxonSignedRankTest {
    037
    038        /** Ranking algorithm. */
    039        private NaturalRanking naturalRanking;
    040
    041        /**
    042         * Create a test instance where NaN's are left in place and ties get
    043         * the average of applicable ranks. Use this unless you are very sure
    044         * of what you are doing.
    045         */
    046        public WilcoxonSignedRankTest() {
    047            naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
    048                    TiesStrategy.AVERAGE);
    049        }
    050
    051        /**
    052         * Create a test instance using the given strategies for NaN's and ties.
    053         * Only use this if you are sure of what you are doing.
    054         *
    055         * @param nanStrategy
    056         *            specifies the strategy that should be used for Double.NaN's
    057         * @param tiesStrategy
    058         *            specifies the strategy that should be used for ties
    059         */
    060        public WilcoxonSignedRankTest(final NaNStrategy nanStrategy,
    061                                      final TiesStrategy tiesStrategy) {
    062            naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
    063        }
    064
    065        /**
    066         * Ensures that the provided arrays fulfills the assumptions.
    067         *
    068         * @param x first sample
    069         * @param y second sample
    070         * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
    071         * @throws NoDataException if {@code x} or {@code y} are zero-length.
    072         * @throws DimensionMismatchException if {@code x} and {@code y} do not
    073         * have the same length.
    074         */
    075        private void ensureDataConformance(final double[] x, final double[] y)
    076            throws NullArgumentException, NoDataException, DimensionMismatchException {
    077
    078            if (x == null ||
    079                y == null) {
    080                    throw new NullArgumentException();
    081            }
    082            if (x.length == 0 ||
    083                y.length == 0) {
    084                throw new NoDataException();
    085            }
    086            if (y.length != x.length) {
    087                throw new DimensionMismatchException(y.length, x.length);
    088            }
    089        }
    090
    091        /**
    092         * Calculates y[i] - x[i] for all i
    093         *
    094         * @param x first sample
    095         * @param y second sample
    096         * @return z = y - x
    097         */
    098        private double[] calculateDifferences(final double[] x, final double[] y) {
    099
    100            final double[] z = new double[x.length];
    101
    102            for (int i = 0; i < x.length; ++i) {
    103                z[i] = y[i] - x[i];
    104            }
    105
    106            return z;
    107        }
    108
    109        /**
    110         * Calculates |z[i]| for all i
    111         *
    112         * @param z sample
    113         * @return |z|
    114         * @throws NullArgumentException if {@code z} is {@code null}
    115         * @throws NoDataException if {@code z} is zero-length.
    116         */
    117        private double[] calculateAbsoluteDifferences(final double[] z)
    118            throws NullArgumentException, NoDataException {
    119
    120            if (z == null) {
    121                throw new NullArgumentException();
    122            }
    123
    124            if (z.length == 0) {
    125                throw new NoDataException();
    126            }
    127
    128            final double[] zAbs = new double[z.length];
    129
    130            for (int i = 0; i < z.length; ++i) {
    131                zAbs[i] = FastMath.abs(z[i]);
    132            }
    133
    134            return zAbs;
    135        }
    136
    137        /**
    138         * Computes the <a
    139         * href="http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test">
    140         * Wilcoxon signed ranked statistic</a> comparing mean for two related
    141         * samples or repeated measurements on a single sample.
    142         * <p>
    143         * This statistic can be used to perform a Wilcoxon signed ranked test
    144         * evaluating the null hypothesis that the two related samples or repeated
    145         * measurements on a single sample has equal mean.
    146         * </p>
    147         * <p>
    148         * Let X<sub>i</sub> denote the i'th individual of the first sample and
    149         * Y<sub>i</sub> the related i'th individual in the second sample. Let
    150         * Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
    151         * </p>
    152         * <p>
    153         * <strong>Preconditions</strong>:
    154         * <ul>
    155         * <li>The differences Z<sub>i</sub> must be independent.</li>
    156         * <li>Each Z<sub>i</sub> comes from a continuous population (they must be
    157         * identical) and is symmetric about a common median.</li>
    158         * <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
    159         * ordered, so the comparisons greater than, less than, and equal to are
    160         * meaningful.</li>
    161         * </ul>
    162         * </p>
    163         *
    164         * @param x the first sample
    165         * @param y the second sample
    166         * @return wilcoxonSignedRank statistic (the larger of W+ and W-)
    167         * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
    168         * @throws NoDataException if {@code x} or {@code y} are zero-length.
    169         * @throws DimensionMismatchException if {@code x} and {@code y} do not
    170         * have the same length.
    171         */
    172        public double wilcoxonSignedRank(final double[] x, final double[] y)
    173            throws NullArgumentException, NoDataException, DimensionMismatchException {
    174
    175            ensureDataConformance(x, y);
    176
    177            // throws IllegalArgumentException if x and y are not correctly
    178            // specified
    179            final double[] z = calculateDifferences(x, y);
    180            final double[] zAbs = calculateAbsoluteDifferences(z);
    181
    182            final double[] ranks = naturalRanking.rank(zAbs);
    183
    184            double Wplus = 0;
    185
    186            for (int i = 0; i < z.length; ++i) {
    187                if (z[i] > 0) {
    188                    Wplus += ranks[i];
    189                }
    190            }
    191
    192            final int N = x.length;
    193            final double Wminus = (((double) (N * (N + 1))) / 2.0) - Wplus;
    194
    195            return FastMath.max(Wplus, Wminus);
    196        }
    197
    198        /**
    199         * Algorithm inspired by
    200         * http://www.fon.hum.uva.nl/Service/Statistics/Signed_Rank_Algorihms.html#C
    201         * by Rob van Son, Institute of Phonetic Sciences & IFOTT,
    202         * University of Amsterdam
    203         *
    204         * @param Wmax largest Wilcoxon signed rank value
    205         * @param N number of subjects (corresponding to x.length)
    206         * @return two-sided exact p-value
    207         */
    208        private double calculateExactPValue(final double Wmax, final int N) {
    209
    210            // Total number of outcomes (equal to 2^N but a lot faster)
    211            final int m = 1 << N;
    212
    213            int largerRankSums = 0;
    214
    215            for (int i = 0; i < m; ++i) {
    216                int rankSum = 0;
    217
    218                // Generate all possible rank sums
    219                for (int j = 0; j < N; ++j) {
    220
    221                    // (i >> j) & 1 extract i's j-th bit from the right
    222                    if (((i >> j) & 1) == 1) {
    223                        rankSum += j + 1;
    224                    }
    225                }
    226
    227                if (rankSum >= Wmax) {
    228                    ++largerRankSums;
    229                }
    230            }
    231
    232            /*
    233             * largerRankSums / m gives the one-sided p-value, so it's multiplied
    234             * with 2 to get the two-sided p-value
    235             */
    236            return 2 * ((double) largerRankSums) / ((double) m);
    237        }
    238
    239        /**
    240         * @param Wmin smallest Wilcoxon signed rank value
    241         * @param N number of subjects (corresponding to x.length)
    242         * @return two-sided asymptotic p-value
    243         */
    244        private double calculateAsymptoticPValue(final double Wmin, final int N) {
    245
    246            final double ES = (double) (N * (N + 1)) / 4.0;
    247
    248            /* Same as (but saves computations):
    249             * final double VarW = ((double) (N * (N + 1) * (2*N + 1))) / 24;
    250             */
    251            final double VarS = ES * ((double) (2 * N + 1) / 6.0);
    252
    253            // - 0.5 is a continuity correction
    254            final double z = (Wmin - ES - 0.5) / FastMath.sqrt(VarS);
    255
    256            // No try-catch or advertised exception because args are valid
    257            final NormalDistribution standardNormal = new NormalDistribution(0, 1);
    258
    259            return 2*standardNormal.cumulativeProbability(z);
    260        }
    261
    262        /**
    263         * Returns the <i>observed significance level</i>, or <a href=
    264         * "http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">
    265         * p-value</a>, associated with a <a
    266         * href="http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test">
    267         * Wilcoxon signed ranked statistic</a> comparing mean for two related
    268         * samples or repeated measurements on a single sample.
    269         * <p>
    270         * Let X<sub>i</sub> denote the i'th individual of the first sample and
    271         * Y<sub>i</sub> the related i'th individual in the second sample. Let
    272         * Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
    273         * </p>
    274         * <p>
    275         * <strong>Preconditions</strong>:
    276         * <ul>
    277         * <li>The differences Z<sub>i</sub> must be independent.</li>
    278         * <li>Each Z<sub>i</sub> comes from a continuous population (they must be
    279         * identical) and is symmetric about a common median.</li>
    280         * <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
    281         * ordered, so the comparisons greater than, less than, and equal to are
    282         * meaningful.</li>
    283         * </ul>
    284         * </p>
    285         *
    286         * @param x the first sample
    287         * @param y the second sample
    288         * @param exactPValue
    289         *            if the exact p-value is wanted (only works for x.length <= 30,
    290         *            if true and x.length > 30, this is ignored because
    291         *            calculations may take too long)
    292         * @return p-value
    293         * @throws NullArgumentException if {@code x} or {@code y} are {@code null}.
    294         * @throws NoDataException if {@code x} or {@code y} are zero-length.
    295         * @throws DimensionMismatchException if {@code x} and {@code y} do not
    296         * have the same length.
    297         * @throws NumberIsTooLargeException if {@code exactPValue} is {@code true}
    298         * and {@code x.length} > 30
    299         * @throws ConvergenceException if the p-value can not be computed due to
    300         * a convergence error
    301         * @throws MaxCountExceededException if the maximum number of iterations
    302         * is exceeded
    303         */
    304        public double wilcoxonSignedRankTest(final double[] x, final double[] y,
    305                                             final boolean exactPValue)
    306            throws NullArgumentException, NoDataException, DimensionMismatchException,
    307            NumberIsTooLargeException, ConvergenceException, MaxCountExceededException {
    308
    309            ensureDataConformance(x, y);
    310
    311            final int N = x.length;
    312            final double Wmax = wilcoxonSignedRank(x, y);
    313
    314            if (exactPValue && N > 30) {
    315                throw new NumberIsTooLargeException(N, 30, true);
    316            }
    317
    318            if (exactPValue) {
    319                return calculateExactPValue(Wmax, N);
    320            } else {
    321                final double Wmin = ( (double)(N*(N+1)) / 2.0 ) - Wmax;
    322                return calculateAsymptoticPValue(Wmin, N);
    323            }
    324        }
    325    }