x = c(82, 79, 81, 79, 77, 79, 79, 78, 79, 82, 76, 73, 64)
y = c(84, 86, 85, 82, 77, 76, 77, 80, 83, 81, 78, 78, 78)

Tn = function(x, y) {
  mean(x) - mean(y)
}

butstrep_p_vrednost = function(x, y, B = 5000){
  m = length(x); n = length(y)
  T_obs = Tn(x, y)
  
  T_B = rep(0, B)
  for(i in 1:B){
    xB = sample(x - mean(x) + mean(c(x, y)), size = m, replace = T)
    yB = sample(y - mean(y) + mean(c(x, y)), size = n, replace = T)
    T_B[i] = Tn(xB, yB)
  }
  
  # odbacujemo ako Tn <= c
  p_vrednost = mean(T_B <= T_obs)
  # p_vrednost = ecdf(T_B)(T_obs)
  return(p_vrednost)
}

butstrep_p_vrednost(x, y)