Couples in a Line and the Lonesome King

The Riddler for November 18, 2016. The Express is a tidy inclusion-exclusion count; the Classic is an elimination game whose answer refuses to settle down.

Riddler Express

Five couples line up at random for a photo, all $10!$ orderings equally likely. What is the probability that no person stands next to his or her partner?

The Riddler, FiveThirtyEight, November 18, 2016(original post)

Solution

It is easier to count the lineups we do not want, the ones with at least one couple together, and subtract. The trouble is that “at least one” double-counts lineups where several couples are together, which is exactly what inclusion-exclusion is built to fix.

If we insist a particular set of $k$ couples each stand as a pair, glue each chosen pair into a single block. That leaves $10-k$ items to arrange in $(10-k)!$ ways, and each glued pair can be in either internal order, a factor $2^k$. There are $\binom{5}{k}$ ways to choose which $k$ couples are glued. Summing with alternating signs over how many couples we force together gives the probability that none are: $$P = \sum_{k=0}^{5} (-1)^k \binom{5}{k}\, 2^k\, \frac{(10-k)!}{10!}.$$ Evaluating the six terms, $$P = 1 - 1 + \tfrac{4}{9} - \tfrac19 + \tfrac{1}{63} - \tfrac{1}{945} = \boxed{\dfrac{47}{135}} \approx 34.8\%.$$

The computation

Encode it two ways: the exact inclusion-exclusion sum with rational arithmetic, and a direct simulation of random lineups that checks the no-adjacent-partner condition.

from fractions import Fraction as F
from math import comb, factorial as fac
import random

P = sum(F((-1)**k * comb(5, k) * 2**k * fac(10 - k), fac(10)) for k in range(6))
print("exact:", P, "=", round(float(P), 4))

rng = random.Random(0); good = 0; T = 2_000_000
partner = lambda x: x ^ 1                      # couples (0,1),(2,3),(4,5),(6,7),(8,9)
for _ in range(T):
    p = list(range(10)); rng.shuffle(p)
    if all(partner(p[i]) != p[i+1] for i in range(9)): good += 1
print("simulated:", round(good / T, 4))
# exact: 47/135 = 0.3481
# simulated: 0.3479

Riddler Classic

The lonely King of Solitaria picks a prince or princess by a game. On the village green, every remaining subject simultaneously points at a random other remaining subject; everyone who is pointed at is eliminated. This repeats. If exactly one subject is ever left, that subject wins; if a round wipes out everyone still standing, nobody wins and the king stays lonely. With $56{,}000$ subjects, is a winner crowned more often than not?

Extra credit: how does the answer behave for a kingdom of arbitrary size?

The Riddler, FiveThirtyEight, November 18, 2016(original post)

Solution

Let $f(n)$ be the probability that, starting a round with $n$ subjects, the game ends with exactly one winner. A round turns $n$ subjects into some number of survivors, the subjects nobody pointed at, and then the same question recurs on the survivors, so $$f(n) = \sum_{s} \Pr(S = s \mid n)\, f(s), \qquad f(1) = 1,\ f(0) = 0,$$ where $S$ is the number of survivors. A subject survives a round only if none of the other $n-1$ subjects points at them, and since each points at a uniform one of the $n-1$ others, the chance a given subject survives is $\bigl(\tfrac{n-2}{n-1}\bigr)^{n-1}$, which tends to $1/e$. So a round keeps roughly a fraction $1/e$ of the players, and the whole game is decided by the luck of the final few rounds, where it can land on $1$ (a winner) or crash to $0$ (nobody).

There is no clean closed form, and the recurrence’s survivor distribution must be evaluated numerically. Doing so, for $56{,}000$ subjects the chance a winner emerges is $$\boxed{\approx 0.48},$$ just under one half, so it is slightly more likely that nobody wins and the king stays lonely.

Extra Credit

Arbitrary kingdom size.

Solution

The genuinely strange part is what happens as $n$ grows: the win probability does not converge. It hovers near $\tfrac12$ but oscillates forever, a little above for some population sizes and a little below for others, with no tidy limit. The shrink-by-$1/e$ dynamics mean the endgame depends on the fractional part of $\log n$ in a way that never settles, behaviour studied under the name “group Russian roulette”. The honest answer is that the value is computed per $n$, and the only universal statement is the non-convergent oscillation around $\tfrac12$.

The computation

Two encodings. For small $n$, evaluate the recurrence exactly by enumerating every choice profile to get the survivor distribution, which validates the model. For $n = 56{,}000$, the state is far too large to enumerate, so simulate the actual rounds and measure the win frequency.

import numpy as np
def simulate(n, trials, seed=0):
    rng = np.random.default_rng(seed); wins = 0
    for _ in range(trials):
        m = n
        while m > 1:
            picks = rng.integers(0, m - 1, size=m); idx = np.arange(m)
            picks = np.where(picks >= idx, picks + 1, picks)   # point at a non-self subject
            m -= len(np.unique(picks))                          # survivors = those unchosen
        wins += (m == 1)
    return wins / trials

from fractions import Fraction as F
from functools import lru_cache
import itertools
def survivor_dist(n):                       # exact P(S=s), tiny n only
    others = [[j for j in range(n) if j != i] for i in range(n)]
    counts = {}
    for prof in itertools.product(*others):
        s = n - len(set(prof)); counts[s] = counts.get(s, 0) + 1
    tot = (n - 1) ** n
    return {s: F(c, tot) for s, c in counts.items()}
@lru_cache(None)
def f(n):
    if n <= 1: return F(n)
    return sum(p * f(s) for s, p in survivor_dist(n).items())

for n in (3, 4, 5, 6):
    print(f"n={n}: exact f={float(f(n)):.4f}  sim={simulate(n, 200000):.4f}")
print("n=56000:", round(simulate(56000, 8000, seed=12345), 3))
# n=3: exact f=0.7500  sim=0.7494
# n=4: exact f=0.5926  sim=0.5923
# n=5: exact f=0.4688  sim=0.4682
# n=6: exact f=0.4161  sim=0.4177
# n=56000: 0.481