What's the correct way to simulate a division winner?

Question

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)

Answer 1

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.

EDIT, added on request: I think the above symmetry argument is just fine, but if we want to get our hands dirty with numbers and stuff, we can do as follows.

All the permutations (rankings) are equally likely to occur. There is a total of $16!$. We can find the probability of an event by simply counting how many permutations are in this event. Let's say that the first number in a permutation corresponds to the placement of team $a$ of Division 1. We want to count how many rankings will have team $a$ as the winner of Division 1. To count the amount of rankings fulfilling this criterion, we can sum over the placement of team $a$. Thus, fix the ranking of team $a$, $i$, where $i \in \{1,...,16\}$. With $i$ fixed we have to have that the three other teams of Division 1 have a smaller ranking than team $a$. These rankings can be chosen in $(i-1), (i-2)$ and $(i-3)$ ways, respectively. The last $12$ rankings do not affect whether or not team $a$ will be the Division 1-winner, thus we can distribute them in $12!$ different ways with the last $12$ remaining placements. This gives us, summing over $i$,

$ \sum_{i=4}^{16} 12!(i-1)(i-2)(i-3) = 15!4 $.

Naturally, team $a$ cannot win with a placement strictly less than fourth, and these cases are omitted. As we have $16!$ different rankings in total, this gives us a probability of $1/4$ of team $a$ winning its division by simply dividing the number of rankings with team $a$ as the winner with the total number of rankings.