秩和检验(rank sum test)及其Python实现

秩和检验(rank sum test)及其Python实现

不是成对的符号秩和检验

秩和检验

假设我们有两组数据A和B，每组数据的样本数为$n_A$何$n_B$，我们想要检验的空假设是A与B具有相同的分布。
首先我们给出秩的定义，如果将一组数据按从小到达排列，那么某观测值的顺序就是该观察值的秩。如假设有一组数据：1, 21, 9, 4, 3，将其按照从小到大的顺序排列，结果为1, 3, 4, 9, 21。那么观测值1的秩就是1，观测值9的秩就是4，以此类推。Wilcoxon 秩和检验的统计量就是一组数据所有观测值的秩的和。这里，我们定义${\omega }_A$为数据集A中的秩和。
例，收集到数据集如下：
A: 8.50,9.48,8.65,8.16,8.83,7.76,8.63
B: 8.27,8.20,8.25,8.14,9.00,8.10,7.20,8.32,7.70

将两组数据混合在一起进行排序，得到序列：7.2 , 7.7 , 7.76, 8.1 , 8.14, 8.16, 8.2 , 8.25, 8.27, 8.32, 8.5 , 8.63, 8.65, 8.83, 9. , 9.48
则可以得到统计量 ${\omega }_A=3+6+11+12+13+14+16=75$
接下来的问题是如何得到统计量对应的p值。如果空假设$H_0$成立，即A与B的数据来自同一分布，那么A中数据的秩应该均匀的，我们可以通过秩和检验表获得$W_A$，下面根据备选假设的讨论如下：

1. 如果备选假设是$H_1：A>B$，那么如果备选假设成立，${\omega }_A>W_A$，因而对应的p值为 $p({\omega }_A<=W_A)$
2. 如果备选假设是$H_1：A<B$，那么如果备选假设成立，${\omega }_A<W_A$，因而对应的p值为 $p({\omega }_A>=W_A)$
3. 对于备选假设$H_1:A\ne B$，应该进行双边检验，如果${\omega }_A$偏小，对应的p值应该检验为$2\ast p({\omega }_A>=W_A)$，如果${\omega }_A$偏大，对应的p值应该检验为$2\ast p({\omega }_A<=W_A)$

其中$W_A$通过查找秩和表得到。

Python实现

为了更加直观的理解秩和检验的过程，我们不适用查表的方法来得到$W_A$，而是通过仿真的方法来获得近似的$W_A$的分布。仿真的代码如下：

import numpy as np
import matplotlib.pyplot as plt

a=np.array([8.50,9.48,8.65,8.16,8.83,7.76,8.63])
b = np.array([8.27,8.20,8.25,8.14,9.00,8.10,7.20,8.32,7.70])
c=np.hstack((a,b))
c.sort()

na = a.shape[0]
nb = b.shape[0]

# 仿真计算H0成立的前提下，WA的分布
N = 100000

randrk = np.zeros((N,))

for i in range(N):
    tmp = np.random.permutation(na+nb) + 1 #出现顺序随机，注意这里需要加1
    randrk[i] = tmp[:na].sum()    
plt.hist(randrk, bins=10)

p = 2*randrk[randrk>=75].shape[0]/N # 双边检验的t值
 #=> 0.1178

这里我们观测到的${\omega }_A$为75，进行双边检验得到的p值为：0.1178。我们与scipy.stats库中的秩和检验进行对比，代码如下：

t, p = stats.ranksums(a, b) # p=>0.1

注意，这里标准库的算法假设数据量足够大，利用${\omega }_A$近似服从正态分布得到的结果。均值和方差如下：