Simulate a tournament with N teams

Question

Based on a programming problem I saw recently, I am wondering if it is possible to simulate the following situation more efficiently:

• There are $N$ teams playing in a tournament.
• Each team plays every other team $K$ times.
• Each team has $1/2$ chance of winning every game it plays.
• There are no draws: in each game, one game wins and one game loses.
• All game outcomes are independent.
• After the tournament is complete, vector $S$ contains $N$ integers; $S[i]$ gives the number of wins by team $i$.

I am looking for a probability distribution with parameters $(N, K)$ with a typical realization $S$. I can simulate $S$ by drawing $K * N * (N - 1)$ Bernoulli random variables (see code below).

Question: Is there a known distribution (as a function of $(N, K)$) that I can sample from to get realizations $S$?.

Here is a sample Python code for $N = 8$ and $K = 3$.

from itertools import combinations
import random
random.seed(0)

n_teams = 8
k_playoffs = 3

def simulate_naive(n_teams=n_teams, k_playoffs=k_playoffs):
  # store score for each team
  # eventually, the score[i] contains the number of times
  # team i has won
  scores = [0] * n_teams
  # iterate over all pairs of teams
  for i, j in combinations(range(n_teams), 2):
    # each pair of teams plays k times
    for k in range(k_playoffs):
      # within each game, each team has 50/50 chance of winning
      if random.random() < 0.5: # condition is met 50% of the time
        scores[i] += 1
      else:
        scores[j] += 1
  return scores

example_scores = simulate_naive()
print(example_scores)
# [8, 9, 13, 4, 14, 11, 12, 13]

Answer 1

First thing to note is that there are $KN(N-1)/2$ Bernoulli trials to simulate, not $KN(N-1)$ since we are dealing with unique pairs of $N$ teams. For instance, if there were only 2 teams then clearly only $K$ games would be played, not $2K$.

Now, among each pair of teams, we can simulate all $K$ games by using a $Binomial(K,.5)$ distribution instead of performing $K$ separate $Bernoulli(.5)$ simulations. Unfortunately, I do not think we could do better due to the following reason. Let $S_{ij}$ denote the number of games that Team $i$ won against Team $j$. As mentioned before, we can simulate this by setting $S_{ij} \sim Binomial(K,.5)$. Now we do not need to simulate $S_{ji}$ after simulating $S_{ij}$ since the number of games Team $j$ won against Team $i$ is $K-S_{ij}$. Therefore, we could simulate $N(N-1)/2$ random variables ($\{S_{ij}: i<j; i,j \in \{1,\cdots,N\}\}$).

Now it may be asked if we can simulate less. For instance, when $N=3$ we could simulate $S_{12}, S_{13},$ and $S_{23}$. Let $S_i$ denote the number of games Team $i$ won in total. Now the number of games that Team 1 won is $S_1 = S_{12} + S_{13}$. Hence, one may be led to opine that it is possible to simulate $S_1$ as a $Binomial(2K,.5)$ distribution. Then we would need only simulate $S_{23}$, saving a single simulation call. Unfortunately, it would be impossible to solve for $S_2=(K-S_{12})+S_{23}$ and $S_3=(K-S_{13}) + (K-S_{23})$, since given $S_1 \ne 0$ we would not be able to recover $S_{12}$ and $S_{13}$.

I recommend initiating an $N\times N$ matrix with zeroes along the main diagonal. Call this $\boldsymbol{X} = (x_{ij})$, with $x_{ij}$ denoting the element in the $i$th row and $j$th column. Next simulate the $N(N-1)/2$ random variables $S_{ij} \sim Binomial(K,.5)$ for $i<j$ and set $x_{ij}=S_{ij}$. Thus you will have a strictly upper-triangular matrix. Next define a strictly lower triangular $N\times N$ matrix, $\boldsymbol{Y}$, where each element equals $K$. Then define the matrix $\boldsymbol{Z} = \boldsymbol{Y} - \boldsymbol{X}^{\prime} + \boldsymbol{X}$. The row sums of $\boldsymbol{Z}$ will give you $S_1, \cdots S_N$.

Here is some R code which performs quite quickly for $N$ as large as $10,000$. For larger $N$, there may be storage issues with the matrices and the $N(N-1)/2$ length vector of Binomial random variables.

library(Matrix)
set.seed(555)
K = 10
N = 3
X = diag(rep(0,N))
Y = X
S = rbinom(choose(N,2),K,.5)
up.index=which(upper.tri(X))
X[up.index] = S
low.index=which(lower.tri(X))
Y[low.index] = K
Z = Y-t(X)+X
team.points = rowSums(Z)