Example of permutation test in R

permutation-test-example-360892.R

# Code for:
# https://stats.stackexchange.com/questions/360863/probability-based-on-observed-data/360892#360892

set.seed(123)

A <- c(34.5, 34.2, 35, 35.1, 36, 35.2, 35.7, 34.8, 34.9, 34.4)
B <- c(35, 34.2, 36, 35, 34, 34.2, 34.2, 34.5, 34.4)



AB <- c(A, B)

k <- length(A)
n <- length(AB)

Tstat <- function(x) {
  A <- x[1:k]
  B <- x[(k+1):n]
  mean(outer(A,B,">"))
}

Tstat_sim <- replicate(50000, Tstat(sample(AB)))

mean(Tstat(AB) <= Tstat_sim)



sim_dens <- function(x, alpha=0.5) {
  A <- x[1:k]
  B <- x[(k+1):n]
  lines(density(A), col=add.alpha("blue", alpha))
  lines(density(B), col=add.alpha("red", alpha))
}


add.alpha <- function(col, alpha = 1){
  out <- apply(sapply(col, col2rgb)/255, 2L,
               function(x)
                 rgb(x[1L], x[2L], x[3L], alpha=alpha))
  out[is.na(col)] <- NA
  out
}




op <- par(mfrow=c(1,2))

plot(0, 0, col=NA, main="Distributions of both groups\n vs simulated data", ylim=c(0, 1.8), xlim=c(32.5, 37), xlab="", ylab="")
replicate(100, sim_dens(sample(AB), alpha=0.1))
lines(density(A), col="blue", lwd=2)
lines(density(B), col="red", lwd=2)

plot(density(Tstat_sim), main="Null distribution of the test statistic\n vs the observed result")
abline(v=Tstat(AB), lwd=2, col="red")

par(op)