In [1]:
xsample <- c(10, 15, 21)
ysample <- c(20, 25, 30)

In [2]:
mean_dist <- function(x, y) {
    return(abs(mean(x) - mean(y))) 
}

In [3]:
perm_test <- function(xspl, yspl, N=2000) {
    nx <- length(xspl)
    ny <- length(yspl)
    n <- nx + ny

    t0 <- mean_dist(xspl, yspl)
    
    combined <- c(xspl, yspl)
    
    z <- array(,N)
    for (i in 1:N) {
        spl <- sample(combined, n, replace=FALSE)
        x <- spl[1:nx]
        y <- spl[(nx+1):n]
        z[i] <- mean_dist(x, y)
    }
    pval <- (sum(z >= t0))/N
    return(pval)
}

In [4]:
perm_test(xsample, ysample)


Out[4]:
0.2025

In [ ]: