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

butstrep_moc = function(m, n, mu_x, sd_x, mu_y, sd_y, N = 1000, B = 500, alpha = 0.05){
  odb = rep(0, 3)
  
  for (i in 1:N) {
    x = rnorm(m, mu_x, sd_x)
    y = rnorm(n, mu_y, sd_y)
    ts = T_mn(x, y) 
  
   ts1 = rep(0, B)
    for(j in 1:B) {
      uzB = sample(c(x, y), size = m + n, replace = T)
      ts1[j] = T_mn(uzB[1:m], uzB[(m + 1):(m + n)])
    }
  
    ts2 = rep(0, B)
    for(j in 1:B) {
      uzB = sample(c(x - mean(x), y - mean(y)), size = m + n, replace = T)
      ts2[j] = T_mn(uzB[1:m], uzB[(m + 1):(m + n)])
    }
  
    ts3 = rep(0, B)
    for(j in 1:B) {
      xB = sample(x / sd(x), size = m, replace = TRUE)
      yB = sample(y / sd(y), size = n, replace = TRUE)
      ts3[j] = T_mn(xB, yB)
    }
  
    c1 = quantile(abs(ts1), 1 - alpha)
    c2 = quantile(abs(ts2), 1 - alpha)
    c3 = quantile(abs(ts3), 1 - alpha)
    
    odb[1] = odb[1] + I(abs(ts) >= c1)
    odb[2] = odb[2] + I(abs(ts) >= c2)
    odb[3] = odb[3] + I(abs(ts) >= c3)
  }
  return(odb / N)
}

butstrep_moc(m = 20, n = 25, mu_x = 0, sd_x = 2, mu_y = 0, sd_y = 2)
butstrep_moc(m = 20, n = 25, mu_x = 3, sd_x = 2, mu_y = 0, sd_y = 2)
butstrep_moc(m = 20, n = 25, mu_x = 0, sd_x = 2, mu_y = 0, sd_y = 3)
butstrep_moc(m = 20, n = 25, mu_x = 0, sd_x = 2, mu_y = 0, sd_y = 10)