Python's implementation of Permutation Test with permutation number as input

Question

R well-known library for permutation test i.e. perm. The example I'm interested in is this:

 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

Which prints result with p-value: 0.01999933. Note there the function permTS allows us to input number of permutation = 30000. Is there such similar implmentation in Python?

I was looking at Python's perm_stat, but it's not what I'm looking for and seems to be buggy.

Answer 1

This is a possible implementation of permutation test using monte-carlo method:

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

note that given the monte-carlo nature of the algorithm you will not get exact same number on each run:

>>> 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