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    }