It is valid. As I understand your question, we say that all the teams of a division are equally likely to win their division. We want to know if the probability of each team to be the division winner of a simulation is .25. You assign the numbers 1 through 16 to the teams in a uniform manner. Each team will be the division winner if and only if that team has the highest number within the division. By symmetry, this probability is .25, as every team must have the same probability of getting the highest number within the division (this comes from way the numbers are handed out in the first place). Your ${\tt R}$-code could be slimmed down somewhat. ${\tt ceiling(seq\_along(rank)/4)) } $will give the same vector for each run, for instance. You also store a lot of information in each run, as the only thing you really need is the winning team for each division. Finally, you say that the convergence seems kind of slow. The pace will be the same if you do the simulation more straight-forwardly (just simulating for one division at a time): the (true) probabilities are certainly the same and the division results are also for your method independent of other divisions, even within a simulation. A heuristic argument of this is to say that no matter what you know about the other divisions, within one division four numbers have been distributed and only their order count. Thus, information on Divisions 1, 2 and 3 doesn't change anything for Division 4. In the end, the two approaches for simulation boils down to the same thing.