Starting point for fixing MATH-431. Still some validation tests and documentation are missing.
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1053836 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Mikkel Meyer Andersen
parent
392f47dfed
commit
794f28c200
|
|
@ -0,0 +1,96 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import org.apache.commons.math.MathException;
|
||||
|
||||
/**
|
||||
* An interface for Mann-Whitney U test (also called Wilcoxon rank-sum test).
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
public interface MannWhitneyUTest {
|
||||
|
||||
/**
|
||||
* Computes the <a
|
||||
* href="http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U"> Mann-Whitney
|
||||
* U statistic</a> comparing mean for two independent samples possibly of
|
||||
* different length.
|
||||
* <p>
|
||||
* This statistic can be used to perform a Mann-Whitney U test evaluating
|
||||
* the null hypothesis that the two independent samples has equal mean.
|
||||
* </p>
|
||||
* <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 would often have different length.
|
||||
* </p>
|
||||
* <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>
|
||||
* </p>
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @return mannWhitneyU statistic
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
*/
|
||||
double mannWhitneyU(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException;
|
||||
|
||||
/**
|
||||
* 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 statistic</a> comparing mean 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. Note that the
|
||||
* samples would often have different length.
|
||||
* </p>
|
||||
* <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>
|
||||
* </p>
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @param exactPValue
|
||||
* if the exact p-value is wanted (only works for x.length <= 30,
|
||||
* if true and x.length > 30, this is ignored because
|
||||
* calculations may take too long)
|
||||
* @return asymptotic p-value
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
* @throws MathException
|
||||
* if an error occurs computing the p-value
|
||||
*/
|
||||
double mannWhitneyUTest(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException, MathException;
|
||||
}
|
||||
|
|
@ -0,0 +1,199 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import org.apache.commons.math.MathException;
|
||||
import org.apache.commons.math.distribution.NormalDistributionImpl;
|
||||
import org.apache.commons.math.stat.ranking.NaNStrategy;
|
||||
import org.apache.commons.math.stat.ranking.NaturalRanking;
|
||||
import org.apache.commons.math.stat.ranking.TiesStrategy;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
/**
|
||||
* An implementation of the Mann-Whitney U test (also called Wilcoxon rank-sum
|
||||
* test).
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
public class MannWhitneyUTestImpl implements MannWhitneyUTest {
|
||||
private NaturalRanking naturalRanking;
|
||||
|
||||
/**
|
||||
* Create a test instance using where NaN's are left in place and ties get
|
||||
* the average of applicable ranks. Use this unless you are very sure of
|
||||
* what you are doing.
|
||||
*/
|
||||
public MannWhitneyUTestImpl() {
|
||||
naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
|
||||
TiesStrategy.AVERAGE);
|
||||
}
|
||||
|
||||
/**
|
||||
* Create a test instance using the given strategies for NaN's and ties.
|
||||
* Only use this if you are sure of what you are doing.
|
||||
*
|
||||
* @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 MannWhitneyUTestImpl(NaNStrategy nanStrategy,
|
||||
TiesStrategy tiesStrategy) {
|
||||
naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
|
||||
}
|
||||
|
||||
/**
|
||||
* Ensures that the provided arrays fulfills the assumptions.
|
||||
*
|
||||
* @param x
|
||||
* @param y
|
||||
* @throws IllegalArgumentException
|
||||
* if assumptions are not met
|
||||
*/
|
||||
private void ensureDataConformance(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException {
|
||||
if (x == null) {
|
||||
throw new IllegalArgumentException("x must not be null");
|
||||
}
|
||||
|
||||
if (y == null) {
|
||||
throw new IllegalArgumentException("y must not be null");
|
||||
}
|
||||
|
||||
if (x.length == 0) {
|
||||
throw new IllegalArgumentException(
|
||||
"x must contain at least one element");
|
||||
}
|
||||
|
||||
if (y.length == 0) {
|
||||
throw new IllegalArgumentException(
|
||||
"y must contain at least one element");
|
||||
}
|
||||
}
|
||||
|
||||
private double[] concatinateSamples(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;
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @return mannWhitneyU statistic U (maximum of U<sup>x</sup> and U<sup>y</sup>)
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
*/
|
||||
public double mannWhitneyU(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException {
|
||||
|
||||
ensureDataConformance(x, y);
|
||||
|
||||
final double[] z = concatinateSamples(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 - (x.length * (x.length + 1)) / 2;
|
||||
|
||||
/*
|
||||
* It can be shown that U1 + U2 = n1 * n2
|
||||
*/
|
||||
final double U2 = x.length * y.length - U1;
|
||||
|
||||
return FastMath.max(U1, U2);
|
||||
}
|
||||
|
||||
/**
|
||||
* @param Umin
|
||||
* smallest Mann-Whitney U value
|
||||
* @param N
|
||||
* number of subjects (corresponding to x.length)
|
||||
* @return two-sided asymptotic p-value
|
||||
* @throws MathException
|
||||
* if an error occurs computing the p-value
|
||||
*/
|
||||
private double calculateAsymptoticPValue(final double Umin, final int n1,
|
||||
final int n2) throws MathException {
|
||||
|
||||
final int n1n2prod = n1 * n2;
|
||||
|
||||
// http://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U#Normal_approximation
|
||||
final double EU = (double) n1n2prod / 2.0;
|
||||
final double VarU = (double) (n1n2prod * (n1 + n2 + 1)) / 12.0;
|
||||
|
||||
final double z = (Umin - EU) / FastMath.sqrt(VarU);
|
||||
|
||||
final NormalDistributionImpl standardNormal = new NormalDistributionImpl(
|
||||
0, 1);
|
||||
|
||||
return 2 * standardNormal.cumulativeProbability(z);
|
||||
}
|
||||
|
||||
/**
|
||||
* Ties give rise to biased variance at the moment. See e.g. <a
|
||||
* href="http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf"
|
||||
* >http://mlsc.lboro.ac.uk/resources/statistics/Mannwhitney.pdf</a>.
|
||||
*
|
||||
* {@inheritDoc}
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @param exactPValue
|
||||
* if the exact p-value is wanted (only for x.length <= 50)
|
||||
* @return asymptotic p-value (biased for samples with ties)
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
* @throws MathException
|
||||
* if an error occurs computing the p-value
|
||||
*/
|
||||
public double mannWhitneyUTest(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException, MathException {
|
||||
|
||||
ensureDataConformance(x, y);
|
||||
|
||||
final double Umax = mannWhitneyU(x, y);
|
||||
|
||||
/*
|
||||
* It can be shown that U1 + U2 = n1 * n2
|
||||
*/
|
||||
final double Umin = x.length * y.length - Umax;
|
||||
|
||||
return calculateAsymptoticPValue(Umin, x.length, y.length);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,107 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import org.apache.commons.math.MathException;
|
||||
|
||||
/**
|
||||
* An interface for Wilcoxon signed-rank test.
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
public interface WilcoxonSignedRankTest {
|
||||
|
||||
/**
|
||||
* Computes the <a
|
||||
* href="http://en.wikipedia.org/wiki/Wilcoxon_signed-rank_test">
|
||||
* Wilcoxon signed ranked statistic</a> comparing mean for two related
|
||||
* samples or repeated measurements on a single sample.
|
||||
* <p>
|
||||
* This statistic can be used to perform a Wilcoxon signed ranked test
|
||||
* evaluating the null hypothesis that the two related samples or repeated
|
||||
* measurements on a single sample has equal mean.
|
||||
* </p>
|
||||
* <p>
|
||||
* Let X<sub>i</sub> denote the i'th individual of the first sample and
|
||||
* Y<sub>i</sub> the related i'th individual in the second sample. Let
|
||||
* Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
|
||||
* </p>
|
||||
* <p>
|
||||
* <strong>Preconditions</strong>:
|
||||
* <ul>
|
||||
* <li>The differences Z<sub>i</sub> must be independent.</li>
|
||||
* <li>Each Z<sub>i</sub> comes from a continuous population (they must be
|
||||
* identical) and is symmetric about a common median.</li>
|
||||
* <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
|
||||
* ordered, so the comparisons greater than, less than, and equal to are
|
||||
* meaningful.</li>
|
||||
* </ul>
|
||||
* </p>
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @return wilcoxonSignedRank statistic
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
*/
|
||||
double wilcoxonSignedRank(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException;
|
||||
|
||||
/**
|
||||
* 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/Wilcoxon_signed-rank_test">
|
||||
* Wilcoxon signed ranked statistic</a> comparing mean for two related
|
||||
* samples or repeated measurements on a single sample.
|
||||
* <p>
|
||||
* Let X<sub>i</sub> denote the i'th individual of the first sample and
|
||||
* Y<sub>i</sub> the related i'th individual in the second sample. Let
|
||||
* Z<sub>i</sub> = Y<sub>i</sub> - X<sub>i</sub>.
|
||||
* </p>
|
||||
* <p>
|
||||
* <strong>Preconditions</strong>:
|
||||
* <ul>
|
||||
* <li>The differences Z<sub>i</sub> must be independent.</li>
|
||||
* <li>Each Z<sub>i</sub> comes from a continuous population (they must be
|
||||
* identical) and is symmetric about a common median.</li>
|
||||
* <li>The values that X<sub>i</sub> and Y<sub>i</sub> represent are
|
||||
* ordered, so the comparisons greater than, less than, and equal to are
|
||||
* meaningful.</li>
|
||||
* </ul>
|
||||
* </p>
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @param exactPValue
|
||||
* if the exact p-value is wanted (only works for x.length <= 30,
|
||||
* if true and x.length > 30, this is ignored because
|
||||
* calculations may take too long)
|
||||
* @return p-value
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
* @throws MathException
|
||||
* if an error occurs computing the p-value
|
||||
*/
|
||||
double wilcoxonSignedRankTest(final double[] x, final double[] y,
|
||||
boolean exactPValue) throws IllegalArgumentException,
|
||||
MathException;
|
||||
}
|
||||
|
|
@ -0,0 +1,271 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import org.apache.commons.math.MathException;
|
||||
import org.apache.commons.math.distribution.NormalDistributionImpl;
|
||||
import org.apache.commons.math.stat.ranking.NaNStrategy;
|
||||
import org.apache.commons.math.stat.ranking.NaturalRanking;
|
||||
import org.apache.commons.math.stat.ranking.TiesStrategy;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
|
||||
/**
|
||||
* An implementation of the Wilcoxon signed-rank test.
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
public class WilcoxonSignedRankTestImpl implements WilcoxonSignedRankTest {
|
||||
private NaturalRanking naturalRanking;
|
||||
|
||||
/**
|
||||
* Create a test instance where NaN's are left in place and ties get
|
||||
* the average of applicable ranks. Use this unless you are very sure
|
||||
* of what you are doing.
|
||||
*/
|
||||
public WilcoxonSignedRankTestImpl() {
|
||||
naturalRanking = new NaturalRanking(NaNStrategy.FIXED,
|
||||
TiesStrategy.AVERAGE);
|
||||
}
|
||||
|
||||
/**
|
||||
* Create a test instance using the given strategies for NaN's and ties.
|
||||
* Only use this if you are sure of what you are doing.
|
||||
*
|
||||
* @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 WilcoxonSignedRankTestImpl(NaNStrategy nanStrategy,
|
||||
TiesStrategy tiesStrategy) {
|
||||
naturalRanking = new NaturalRanking(nanStrategy, tiesStrategy);
|
||||
}
|
||||
|
||||
/**
|
||||
* Ensures that the provided arrays fulfills the assumptions.
|
||||
*
|
||||
* @param x
|
||||
* @param y
|
||||
* @throws IllegalArgumentException
|
||||
* if assumptions are not met
|
||||
*/
|
||||
private void ensureDataConformance(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException {
|
||||
if (x == null) {
|
||||
throw new IllegalArgumentException("x must not be null");
|
||||
}
|
||||
|
||||
if (y == null) {
|
||||
throw new IllegalArgumentException("y must not be null");
|
||||
}
|
||||
|
||||
if (x.length != y.length) {
|
||||
throw new IllegalArgumentException(
|
||||
"x and y must contain the same number of elements");
|
||||
}
|
||||
|
||||
if (x.length == 0) {
|
||||
throw new IllegalArgumentException(
|
||||
"x and y must contain at least one element");
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates y[i] - x[i] for all i
|
||||
*
|
||||
* @param x
|
||||
* @param y
|
||||
* @throws IllegalArgumentException
|
||||
* if assumptions are not met
|
||||
*/
|
||||
private double[] calculateDifferences(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException {
|
||||
|
||||
final double[] z = new double[x.length];
|
||||
|
||||
for (int i = 0; i < x.length; ++i) {
|
||||
z[i] = y[i] - x[i];
|
||||
}
|
||||
|
||||
return z;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculates |z[i]| for all i
|
||||
*
|
||||
* @param z
|
||||
* @throws IllegalArgumentException
|
||||
* if assumptions are not met
|
||||
*/
|
||||
private double[] calculateAbsoluteDifferences(final double[] z)
|
||||
throws IllegalArgumentException {
|
||||
if (z == null) {
|
||||
throw new IllegalArgumentException("z must not be null");
|
||||
}
|
||||
|
||||
if (z.length == 0) {
|
||||
throw new IllegalArgumentException(
|
||||
"z must contain at least one element");
|
||||
}
|
||||
|
||||
final double[] zAbs = new double[z.length];
|
||||
|
||||
for (int i = 0; i < z.length; ++i) {
|
||||
zAbs[i] = FastMath.abs(z[i]);
|
||||
}
|
||||
|
||||
return zAbs;
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @return wilcoxonSignedRank statistic (the larger of W+ and W-)
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met
|
||||
*/
|
||||
public double wilcoxonSignedRank(final double[] x, final double[] y)
|
||||
throws IllegalArgumentException {
|
||||
|
||||
ensureDataConformance(x, y);
|
||||
|
||||
// throws IllegalArgumentException if x and y are not correctly
|
||||
// specified
|
||||
final double[] z = calculateDifferences(x, y);
|
||||
final double[] zAbs = calculateAbsoluteDifferences(z);
|
||||
|
||||
final double[] ranks = naturalRanking.rank(zAbs);
|
||||
|
||||
double Wplus = 0;
|
||||
|
||||
for (int i = 0; i < z.length; ++i) {
|
||||
if (z[i] > 0) {
|
||||
Wplus += ranks[i];
|
||||
}
|
||||
}
|
||||
|
||||
final int N = x.length;
|
||||
final double Wminus = (((double) (N * (N + 1))) / 2.0) - Wplus;
|
||||
|
||||
return FastMath.max(Wplus, Wminus);
|
||||
}
|
||||
|
||||
/**
|
||||
* Algorithm inspired by
|
||||
* http://www.fon.hum.uva.nl/Service/Statistics/Signed_Rank_Algorihms.html#C
|
||||
* by Rob van Son, Institute of Phonetic Sciences & IFOTT,
|
||||
* University of Amsterdam
|
||||
*
|
||||
* @param Wmax largest Wilcoxon signed rank value
|
||||
* @param N number of subjects (corresponding to x.length)
|
||||
* @return two-sided exact p-value
|
||||
*/
|
||||
private double calculateExactPValue(final double Wmax, final int N) {
|
||||
|
||||
// Total number of outcomes (equal to 2^N but a lot faster)
|
||||
final int m = 1 << N;
|
||||
|
||||
int largerRankSums = 0;
|
||||
|
||||
for (int i = 0; i < m; ++i) {
|
||||
int rankSum = 0;
|
||||
|
||||
// Generate all possible rank sums
|
||||
for (int j = 0; j < N; ++j) {
|
||||
|
||||
// (i >> j) & 1 extract i's j-th bit from the right
|
||||
if (((i >> j) & 1) == 1) {
|
||||
rankSum += j + 1;
|
||||
}
|
||||
}
|
||||
|
||||
if (rankSum >= Wmax) {
|
||||
++largerRankSums;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* largerRankSums / m gives the one-sided p-value, so it's multiplied
|
||||
* with 2 to get the two-sided p-value
|
||||
*/
|
||||
return 2 * ((double) largerRankSums) / ((double) m);
|
||||
}
|
||||
|
||||
/**
|
||||
* @param Wmin smallest Wilcoxon signed rank value
|
||||
* @param N number of subjects (corresponding to x.length)
|
||||
* @return two-sided asymptotic p-value
|
||||
* @throws MathException if an error occurs computing the p-value
|
||||
*/
|
||||
private double calculateAsymptoticPValue(final double Wmin, final int N) throws MathException {
|
||||
|
||||
final double ES = (double) (N * (N + 1)) / 4.0;
|
||||
|
||||
/* Same as (but saves computations):
|
||||
* final double VarW = ((double) (N * (N + 1) * (2*N + 1))) / 24;
|
||||
*/
|
||||
final double VarS = ES * ((double) (2 * N + 1) / 6.0);
|
||||
|
||||
// - 0.5 is a continuity correction
|
||||
final double z = (Wmin - ES - 0.5) / FastMath.sqrt(VarS);
|
||||
|
||||
final NormalDistributionImpl standardNormal = new NormalDistributionImpl(0, 1);
|
||||
|
||||
return 2*standardNormal.cumulativeProbability(z);
|
||||
}
|
||||
|
||||
/**
|
||||
* {@inheritDoc}
|
||||
*
|
||||
* @param x
|
||||
* the first sample
|
||||
* @param y
|
||||
* the second sample
|
||||
* @param exactPValue
|
||||
* if the exact p-value is wanted (only for x.length <= 30)
|
||||
* @return p-value
|
||||
* @throws IllegalArgumentException
|
||||
* if preconditions are not met or exact p-value is wanted for
|
||||
* when x.length > 30
|
||||
* @throws MathException
|
||||
* if an error occurs computing the p-value
|
||||
*/
|
||||
public double wilcoxonSignedRankTest(final double[] x, final double[] y,
|
||||
boolean exactPValue) throws IllegalArgumentException,
|
||||
MathException {
|
||||
|
||||
ensureDataConformance(x, y);
|
||||
|
||||
final int N = x.length;
|
||||
final double Wmax = wilcoxonSignedRank(x, y);
|
||||
|
||||
if (exactPValue && N > 30) {
|
||||
throw new IllegalArgumentException("Exact test can only be made for N <= 30.");
|
||||
}
|
||||
|
||||
if (exactPValue) {
|
||||
return calculateExactPValue(Wmax, N);
|
||||
} else {
|
||||
final double Wmin = ( (double)(N*(N+1)) / 2.0 ) - Wmax;
|
||||
return calculateAsymptoticPValue(Wmin, N);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,101 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import junit.framework.TestCase;
|
||||
|
||||
/**
|
||||
* Test cases for the ChiSquareTestImpl class.
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
|
||||
public class MannWhitneyUTestTest extends TestCase {
|
||||
|
||||
protected MannWhitneyUTest testStatistic = new MannWhitneyUTestImpl();
|
||||
|
||||
public MannWhitneyUTestTest(String name) {
|
||||
super(name);
|
||||
}
|
||||
|
||||
public void testMannWhitneyUSimple() throws Exception {
|
||||
/* Target values computed using R version 2.11.1
|
||||
* x <- c(19, 22, 16, 29, 24)
|
||||
* y <- c(20, 11, 17, 12)
|
||||
* wilcox.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE, correct = FALSE)
|
||||
* W = 17, p-value = 0.08641
|
||||
*/
|
||||
final double x[] = {19, 22, 16, 29, 24};
|
||||
final double y[] = {20, 11, 17, 12};
|
||||
|
||||
assertEquals(17, testStatistic.mannWhitneyU(x, y), 1e-10);
|
||||
assertEquals(0.08641, testStatistic.mannWhitneyUTest(x, y), 1e-5);
|
||||
}
|
||||
|
||||
|
||||
public void testMannWhitneyUInputValidation() throws Exception {
|
||||
/* Samples must be present, i.e. length > 0
|
||||
*/
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(new double[] { }, new double[] { 1.0 });
|
||||
fail("x does not contain samples (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(new double[] { 1.0 }, new double[] { });
|
||||
fail("y does not contain samples (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/*
|
||||
* x and y is null
|
||||
*/
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(null, null);
|
||||
fail("x and y is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(null, null);
|
||||
fail("x and y is null (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/*
|
||||
* x or y is null
|
||||
*/
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(null, new double[] { 1.0 });
|
||||
fail("x is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.mannWhitneyUTest(new double[] { 1.0 }, null);
|
||||
fail("y is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,180 @@
|
|||
/*
|
||||
* 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
|
||||
*
|
||||
* http://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.
|
||||
*/
|
||||
package org.apache.commons.math.stat.inference;
|
||||
|
||||
import junit.framework.TestCase;
|
||||
|
||||
/**
|
||||
* Test cases for the ChiSquareTestImpl class.
|
||||
*
|
||||
* @version $Revision: $ $Date: $
|
||||
*/
|
||||
|
||||
public class WilcoxonSignedRankTestTest extends TestCase {
|
||||
|
||||
protected WilcoxonSignedRankTest testStatistic = new WilcoxonSignedRankTestImpl();
|
||||
|
||||
public WilcoxonSignedRankTestTest(String name) {
|
||||
super(name);
|
||||
}
|
||||
|
||||
public void testWilcoxonSignedRankSimple() throws Exception {
|
||||
/* Target values computed using R version 2.11.1
|
||||
* x <- c(1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
|
||||
* y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)
|
||||
*/
|
||||
final double x[] = {1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30};
|
||||
final double y[] = {0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29};
|
||||
|
||||
/* EXACT:
|
||||
* wilcox.test(x, y, alternative = "two.sided", mu = 0, paired = TRUE, exact = TRUE, correct = FALSE)
|
||||
* V = 40, p-value = 0.03906
|
||||
*
|
||||
* Corresponds to the value obtained in R.
|
||||
*/
|
||||
assertEquals(40, testStatistic.wilcoxonSignedRank(x, y), 1e-10);
|
||||
assertEquals(0.03906, testStatistic.wilcoxonSignedRankTest(x, y, true), 1e-5);
|
||||
|
||||
/* ASYMPTOTIC:
|
||||
* wilcox.test(x, y, alternative = "two.sided", mu = 0, paired = TRUE, exact = FALSE, correct = FALSE)
|
||||
* V = 40, p-value = 0.03815
|
||||
*
|
||||
* This is not entirely the same due to different corrects,
|
||||
* e.g. http://mlsc.lboro.ac.uk/resources/statistics/wsrt.pdf
|
||||
* and src/library/stats/R/wilcox.test.R in the R source
|
||||
*/
|
||||
assertEquals(40, testStatistic.wilcoxonSignedRank(x, y), 1e-10);
|
||||
assertEquals(0.0329693812, testStatistic.wilcoxonSignedRankTest(x, y, false), 1e-10);
|
||||
}
|
||||
|
||||
public void testWilcoxonSignedRankInputValidation() throws Exception {
|
||||
/*
|
||||
* Exact only for sample size <= 30
|
||||
*/
|
||||
final double[] x1 = new double[30];
|
||||
final double[] x2 = new double[31];
|
||||
final double[] y1 = new double[30];
|
||||
final double[] y2 = new double[31];
|
||||
for (int i = 0; i < 30; ++i) {
|
||||
x1[i] = x2[i] = y1[i] = y2[i] = i;
|
||||
}
|
||||
|
||||
// Exactly 30 is okay
|
||||
testStatistic.wilcoxonSignedRankTest(x1, y1, true);
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(x2, y2, true);
|
||||
fail("More than 30 samples and exact chosen, IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/* Samples must be present, i.e. length > 0
|
||||
*/
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { }, new double[] { 1.0 }, true);
|
||||
fail("x does not contain samples (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { }, new double[] { 1.0 }, false);
|
||||
fail("x does not contain samples (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0 }, new double[] { }, true);
|
||||
fail("y does not contain samples (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0 }, new double[] { }, false);
|
||||
fail("y does not contain samples (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/* Samples not same size, i.e. cannot be pairred
|
||||
*/
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0, 2.0 }, new double[] { 3.0 }, true);
|
||||
fail("x and y not same size (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0, 2.0 }, new double[] { 3.0 }, false);
|
||||
fail("x and y not same size (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/*
|
||||
* x and y is null
|
||||
*/
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(null, null, true);
|
||||
fail("x and y is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(null, null, false);
|
||||
fail("x and y is null (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
/*
|
||||
* x or y is null
|
||||
*/
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(null, new double[] { 1.0 }, true);
|
||||
fail("x is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(null, new double[] { 1.0 }, false);
|
||||
fail("x is null (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0 }, null, true);
|
||||
fail("y is null (exact), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
|
||||
try {
|
||||
testStatistic.wilcoxonSignedRankTest(new double[] { 1.0 }, null, false);
|
||||
fail("y is null (asymptotic), IllegalArgumentException expected");
|
||||
} catch (IllegalArgumentException ex) {
|
||||
// expected
|
||||
}
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue