test_statistika = function(x, y) {
  m = length(x)
  n = length(y)
  return(sqrt((m*n)/(m+n))*(log(var(x)) - log(var(y))))
}

ocena_p_val = function(x, y) {
  Tmn = test_statistika(x, y)
  
  # prvi metod reuzorkovanja
  ts1 = rep(0, B)
  for (i in 1:B) {
    uz_b = sample(c(x, y), size = m + n, replace = T)
    ts1[i] = test_statistika(uz_b[1:m], uz_b[(m+1):(n+m)])
  }
  
  # drugi metod reuzorkovanja
  ts2 = rep(0, B)
  for (i in 1:B) {
    uz_b = sample(c(x - mean(x), y - mean(y)), size = m + n, replace = T)
    ts2[i] = test_statistika(uz_b[1:m], uz_b[(m+1):(n+m)])
  }
  
  # treci metod reuzorkovanja
  ts3 = rep(0, B)
  for (i in 1:B) {
    X = sample(x/sd(x), size = m, replace = T)
    Y = sample(y/sd(y), size = n, replace = T)
    ts3[i] = test_statistika(X, Y)
  }
  
  p_val = rep(0, 3)
  # odbacujemo H0 ako |T_{mn}|>= c
  # c = quantile(abs(niz_ts), 1 - 0.05)
  if (abs(Tmn) > quantile(abs(ts1), 0.95)) {
    p_val[1] = p_val[1] + 1
  }
  if (abs(Tmn) > quantile(abs(ts2), 0.95)) {
    p_val[2] = p_val[2] + 1
  }
  if (abs(Tmn) > quantile(abs(ts3), 0.95)) {
    p_val[3] = p_val[3] + 1
  }
  return(p_val)
}

n = 20; m = 25
B = 500
N = 1000

s = rep(0, 3)
for (i in 1:N) {
  x = rnorm(m)
  y = rnorm(n)
  s = s + ocena_p_val(x,y)
}
s/N

s = rep(0, 3)
for (i in 1:N) {
  x = rnorm(m, 2, 2)
  y = rnorm(n, 2, 2)
  s = s + ocena_p_val(x,y)
}
s/N

s = rep(0, 3)
for (i in 1:N) {
  x = rnorm(m, 10, 2)
  y = rnorm(n, 2, 2)
  s = s + ocena_p_val(x,y)
}
s/N

s = rep(0, 3)
for (i in 1:N) {
  x = rnorm(m, 2, 2)
  y = rnorm(n, 2, 5)
  s = s + ocena_p_val(x,y)
}
s/N

s = rep(0, 3)
for (i in 1:N) {
  x = rnorm(m, 10, 2)
  y = rnorm(n, 2, 5)
  s = s + ocena_p_val(x,y)
}
s/N