秩和检验(rank sum test)及其Python实现
不是成对的符号秩和检验
秩和检验
假设我们有两组数据A和B,每组数据的样本数为
n
A
n_A
nA何
n
B
n_B
nB,我们想要检验的空假设是A与B具有相同的分布。
首先我们给出秩的定义,如果将一组数据按从小到达排列,那么某观测值的顺序就是该观察值的秩。如假设有一组数据:1, 21, 9, 4, 3,将其按照从小到大的顺序排列,结果为1, 3, 4, 9, 21。那么观测值1的秩就是1,观测值9的秩就是4,以此类推。Wilcoxon 秩和检验的统计量就是一组数据所有观测值的秩的和。这里,我们定义
ω
A
\omega_A
ω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
则可以得到统计量
ω
A
=
3
+
6
+
11
+
12
+
13
+
14
+
16
=
75
\omega_A= 3 + 6 + 11 + 12 + 13 + 14 + 16 = 75
ωA=3+6+11+12+13+14+16=75
接下来的问题是如何得到统计量对应的p值。如果空假设
H
0
H_0
H0成立,即A与B的数据来自同一分布,那么A中数据的秩应该均匀的,我们可以通过秩和检验表获得
W
A
W_A
WA,下面根据备选假设的讨论如下:
- 如果备选假设是 H 1 : A > B H_1:A>B H1:A>B,那么如果备选假设成立, ω A > W A \omega_A>W_A ωA>WA,因而对应的p值为 p ( ω A < = W A ) p(\omega_A<=W_A) p(ωA<=WA)
- 如果备选假设是 H 1 : A < B H_1:A<B H1:A<B,那么如果备选假设成立, ω A < W A \omega_A<W_A ωA<WA,因而对应的p值为 p ( ω A > = W A ) p(\omega_A>=W_A) p(ωA>=WA)
- 对于备选假设 H 1 : A ≠ B H_1: A\neq B H1:A=B,应该进行双边检验,如果 ω A \omega_A ωA偏小,对应的p值应该检验为 2 ∗ p ( ω A > = W A ) 2*p(\omega_A>=W_A) 2∗p(ωA>=WA),如果 ω A \omega_A ωA偏大,对应的p值应该检验为 2 ∗ p ( ω A < = W A ) 2*p(\omega_A<=W_A) 2∗p(ωA<=WA)
其中 W A W_A WA通过查找秩和表得到。
Python实现
为了更加直观的理解秩和检验的过程,我们不适用查表的方法来得到 W A W_A WA,而是通过仿真的方法来获得近似的 W A W_A WA的分布。仿真的代码如下:
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
这里我们观测到的
ω
A
\omega_A
ωA为75,进行双边检验得到的p值为:0.1178。我们与scipy.stats
库中的秩和检验进行对比,代码如下:
t, p = stats.ranksums(a, b) # p=>0.1
注意,这里标准库的算法假设数据量足够大,利用
ω
A
\omega_A
ωA近似服从正态分布得到的结果。均值和方差如下: