Python中排列检验的实现，输入为排列数

Question

R语言有一个很有名的库，专门用来做排列检验，叫做 perm。我感兴趣的例子是这样的：

 x <- c(12.6, 11.4, 13.2, 11.2, 9.4, 12.0)
 y <- c(16.4, 14.1, 13.4, 15.4, 14.0, 11.3)
 permTS(x,y, alternative="two.sided", method="exact.mc", control=permControl(nmc=30000))$p.value

这个例子会输出一个结果，p值是：0.01999933。注意，这里使用的函数permTS允许我们输入的排列次数是30000。请问在Python中有没有类似的实现呢？

我查看了Python的 perm_stat，但这并不是我想要的，而且似乎还有一些bug。

Answer 1

这是一个使用蒙特卡洛方法进行排列检验的可能实现：

def exact_mc_perm_test(xs, ys, nmc):
    n, k = len(xs), 0
    diff = np.abs(np.mean(xs) - np.mean(ys))
    zs = np.concatenate([xs, ys])
    for j in range(nmc):
        np.random.shuffle(zs)
        k += diff < np.abs(np.mean(zs[:n]) - np.mean(zs[n:]))
    return k / nmc

请注意，由于这个算法是蒙特卡洛性质的，所以每次运行时你得到的结果可能会有所不同，不能保证每次都是一样的数字：

>>> xs = np.array([12.6, 11.4, 13.2, 11.2, 9.4, 12.0])
>>> ys = np.array([16.4, 14.1, 13.4, 15.4, 14.0, 11.3])
>>> exact_mc_perm_test(xs, ys, 30000)
0.019466666666666667