permutation_test.py

def permutation_simulation(higher, lower, n=100_000):
    xs = higher
    ys = lower
    m = np.mean(xs) - np.mean(ys)
    zs = np.concatenate((xs, ys))
    z_perms = np.random.choice(ys, size=(n, len(zs)))
    x_perms = z_perms[:, 0:len(xs)]
    y_perms = z_perms[:, len(xs)+1:]
    x_perm_means = np.mean(x_perms, 1)
    y_perm_means = np.mean(y_perms, 1)
    return np.sum(x_perm_means - y_perm_means > m)

def permutation_pvalue(higher, lower, n=100_000):
    return permutation_simulation(higher, lower, n) / n