d <- 8        # Sides of the die
throws <- 20  # Number of throws
no_cut <- 6   # Number at least on die
n_lower <- 12 # Lower cutoff
n_upper <- 16 # Upper cutoff

n_sim <- 1000000 # Number of simulations

res <- NULL

res <- replicate(n_sim, {
  x <- sample.int(d, size = throws, replace = TRUE)
  s <- sum(x > no_cut)
  (s >= n_lower && s <= n_upper)
})

sum(res)/n_sim # simulated probability