Esercizio - Tre persone con lo stesso compleanno e due scelte nel bilanciamento

Teoria: Palline e urne - compleanno, carico massimo e approssimazione di PoissonModello palline e urne (m palline, n urne, scelte indipendenti e uniformi): paradosso del compleanno con soglia sqrt(2 n ln 2), urne vuote attese n e^(-m/n), carico massimo ln n / ln ln n con dimostrazione (first moment e union bound), approssimazione di Poisson con il fattore e sqrt(m), potere delle due scelte (ln ln n / ln 2).Palline e urne - compleanno, carico massimo e approssimazione di Poisson → (soglia del compleanno per le terne, potere delle due scelte), Linearità del valore atteso e variabili indicatriciLa linearità $E[X+Y]=E[X]+E[Y]$ vale sempre, anche senza indipendenza; con le variabili indicatrici ($E[I_A]=P(A)$) trasforma un conteggio in una somma di probabilità. Applicazioni: attesa del primo successo $1/p$, collezionista di figurine $nH_n$, punti fissi di una permutazione (media $1$), record (media $H_n$), inversioni ($n(n-1)/4$), attesa condizionata, identità di Wald $E[S_N]=E[N]E[X]$.Linearità del valore atteso e variabili indicatrici →.

Testo.

1. Con $n=365$ giorni, quante persone servono perché il numero atteso di terne con lo stesso compleanno sia almeno $\frac12$? E almeno $1$? Confrontare con la formula $6^{1/3}n^{2/3}$.
2. Calcolare la probabilità esatta che tra $m$ persone almeno tre abbiano lo stesso compleanno e trovare il minimo $m$ per cui supera $\frac12$. (Suggerimento: contare le assegnazioni in cui ogni giorno ha al più $2$ persone con la funzione generatrice esponenziale $(1+x+\frac{x^2}2)^{365}$.)
3. Confrontare il carico massimo di $n$ palline in $n$ urne con una scelta e con due scelte (la pallina va nell'urna meno carica fra due estratte a caso) per $n=10^3,10^4,10^5$, con una simulazione.

Soluzione

1. Terne attese. Le terne di persone sono $\binom m3$, e ognuna ha lo stesso compleanno con probabilità $\frac1{n^2}$ (la seconda e la terza persona devono coincidere con la prima: $\frac1n\cdot\frac1n$). Per la linearità $E[\#\text{terne}]=\binom m3\frac1{365^2}$:

Il valore $\frac12$ si raggiunge a $\boxed{m=75}$ ($\binom{74}3=64\,824$ dà $0{,}487$, $\binom{75}3=67\,525$ dà $0{,}507$; $365^2=133\,225$). Il valore $1$ richiede $\binom m3\ge133\,225$, cioè $m=94$ ($\binom{93}3=129\,766$ dà $0{,}974$, $\binom{94}3=134\,044$ dà $1{,}006$). La formula approssimata $\binom m3\approx\frac{m^3}6=n^2\Rightarrow m=6^{1/3}n^{2/3}=1{,}817\cdot51{,}09=92{,}8$ per il valore $1$.

2. Probabilità esatta. Se $m$ persone si distribuiscono sui $n$ giorni, il numero di assegnazioni (persone distinguibili) in cui ogni giorno ha $0$, $1$ o $2$ persone è $m!\,[x^m]\bigl(1+x+\frac{x^2}2\bigr)^n$: ogni giorno contribuisce $\sum_{j\le2}\frac{x^j}{j!}$ (funzione generatrice esponenziale). Quindi $$P(\text{nessuna terna})=\frac{m!\,[x^m](1+x+x^2/2)^{365}}{365^m}.$$ Il coefficiente si calcola moltiplicando il polinomio per se stesso $365$ volte, troncando al grado $m$. Risultati:

Il minimo è $\boxed{m=88}$ (con $m=87$ si ha $0{,}49945$, appena sotto $\frac12$; l'enumerazione con frazioni esatte per $m=88$ dà $0{,}51107$). Quindi per le terne servono $88$ persone, contro le $23$ delle coppie: la soglia cresce come $n^{2/3}$ anziché $n^{1/2}$. Si noti che a $m=88$ il numero atteso di terne è $0{,}82$, non $1$: la probabilità di almeno una terna ($0{,}51$) è molto inferiore al valore atteso, perché le terne sono molto sovrapposte.

3. Una scelta contro due scelte. Frequenze osservate del carico massimo (100 prove per $n=10^3,10^4$, 30 per $n=10^5$):

Con una scelta il massimo cresce come $\ln n/\ln\ln n$ ($3{,}6;\ 4{,}1;\ 4{,}7$ nei tre casi, con costanti più grandi); con due scelte resta fra $3$ e $4$, in accordo con $\frac{\ln\ln n}{\ln2}+O(1)$ ($2{,}8;\ 3{,}2;\ 3{,}5$ più una costante). Il costo è una seconda scelta casuale e un confronto per ogni pallina.

Codice di verifica

import random
import numpy as np
from math import comb, lgamma, log, exp

def p_terna(m, n=365):
    """P(almeno 3 con lo stesso giorno) = 1 - m! [x^m](1+x+x^2/2)^n / n^m."""
    poly = np.zeros(m + 1); poly[0] = 1.0
    for _ in range(n):                    # moltiplica per (1 + x + x^2/2), troncando al grado m
        nuovo = poly.copy()
        nuovo[1:] += poly[:-1]
        nuovo[2:] += 0.5 * poly[:-2]
        poly = nuovo
    return 1 - exp(log(poly[m]) + lgamma(m + 1) - m * log(n))

print([round(comb(m, 3) / 365 ** 2, 3) for m in (75, 80, 85, 88, 90, 94)])
print([round(p_terna(m), 4) for m in (80, 85, 87, 88, 89, 90)])

def massimo(n, d):
    c = [0] * n
    for _ in range(n):
        i = min((random.randrange(n) for _ in range(d)), key=lambda j: c[j])
        c[i] += 1                         # nell'urna meno carica fra le d estratte
    return max(c)

random.seed(3)
for n in (10 ** 3, 10 ** 4):
    print(n, sorted(set(massimo(n, 1) for _ in range(20))), sorted(set(massimo(n, 2) for _ in range(20))))

Risultato: [0.507, 0.617, 0.741, 0.824, 0.882, 1.006], poi [0.4182, 0.4762, 0.4995, 0.5111, 0.5226, 0.5342] e, per $n=10^3$ e $n=10^4$, due righe con i massimi osservati in 20 prove (cambiano a ogni esecuzione: con una scelta tra $5$ e $8$, con due scelte tra $3$ e $4$), ad esempio 1000 [5, 6] [3, 4] e 10000 [6, 7, 8] [3, 4].