So I'm trying to simulate one winner for each of the 4 NFL divisions, assuming that all teams are of equal ability. It's clear that the probability of each team being the division winner is $1/4$.

I wonder if the following simulation strategy works. Instead of simulating each division, I randomly rank all 16 teams from 1 to 16, then compare them in group of 4. I did run the code, and the probability does converge to .25, but it takes 100,000 simulations to even get kinda close.

Thus, I'm unsure whether the proposed simulation strategy is valid or not?

one_simu <- function() {
  rank <- sample(1:16, replace=F)
  groups <- split(rank, ceiling(seq_along(rank)/4)) # Split into 4 groups of 4
  result <- rep(0, 16)
  for (i in seq_along(rank)) {
    group_of_i <- groups[[ceiling(i/4)] ] # The group that team 1 belongs to
    if (rank[i] == max(group_of_i)) {
      result[i] <- result[i] + 1  
    }
  } 
  return(result)
}

apply(replicate(1000, one_simu()), 1, mean)
apply(replicate(100000, one_simu()), 1, mean)