02 - Variabili aleatorie discrete: Bernoulli, Binomiale, Poisson

Fai un quiz di 20 domande a risposta multipla, ciascuna con quattro opzioni, e rispondi tirando a indovinare. Quante risposte esatte ti aspetti? E con quale probabilità ne azzecchi almeno dieci, abbastanza per essere promosso? Sono domande sul conteggio di successi in prove ripetute, e la loro risposta ha un nome: la distribuzione binomiale.

Qui la ricavi, la calcoli in modo esatto (la media è 5, la varianza 15/4, la probabilità di arrivare a dieci è circa 1,4%) e la confermi con un milione di quiz simulati. Poi vedi come la binomiale diventa, quando il successo è raro, una distribuzione di Poisson, quanto è grande l'errore che si commette, e come la distribuzione geometrica descrive l'attesa del primo successo. Ti serve aver letto il primo articolo della serie sugli assiomi, sull'indipendenza e sulla probabilità condizionata; il resto lo definisco qui.

Dare un numero a ogni esito

Una variabile aleatoria è una regola che associa un numero a ogni esito dell'esperimento: nel quiz, il numero di risposte esatte. Si dice discreta se i valori possibili si possono elencare (0, 1, 2, …). La sua descrizione completa è la funzione di massa di probabilità (in inglese PMF), p(k) = P(X = k): un numero non negativo per ogni valore, con somma 1. Dalla PMF si ricavano due riassunti.

Il valore atteso E[X] = Σ k p(k) è la media dei valori pesata dalle probabilità, il «baricentro» della distribuzione. La varianza Var(X) = E[(X − E[X])²] = E[X²] − (E[X])² misura quanto i valori si disperdono attorno al valore atteso; la sua radice quadrata è la deviazione standard, nella stessa unità di misura di X.

Due proprietà fanno quasi tutto il lavoro nei calcoli che seguono. La prima è la linearità del valore atteso: E[X + Y] = E[X] + E[Y] sempre, anche se X e Y dipendono l'una dall'altra.

La seconda riguarda la varianza, e qui serve cautela: Var(X + Y) = Var(X) + Var(Y) vale quando X e Y sono indipendenti, non in generale (nell'articolo sulla covarianza vedi il termine che manca quando non lo sono). La differenza fra le due regole è la ragione per cui la media di una somma si calcola sempre, mentre la sua varianza richiede di sapere come si relazionano i pezzi.

Bernoulli e binomiale: contare i successi

La variabile più semplice è la Bernoulli di parametro p: vale 1 («successo») con probabilità p e 0 con probabilità q = 1 − p. Il suo valore atteso è p e la sua varianza pq, perché X² = X e quindi E[X²] = p. Una domanda del quiz, risposta a caso, è una Bernoulli con p = 1/4. Se sommi n Bernoulli indipendenti con lo stesso p ottieni il numero di successi Y, che ha distribuzione binomiale Bin(n, p).

La PMF si ricava contando: una particolare sequenza con k successi e n − k fallimenti ha probabilità p^kq^n−k per l'indipendenza, e le sequenze con k successi sono C(n, k) = n! / (k!(n − k)!), da cui P(Y = k) = C(n, k) p^kq^n−k.

La somma di n Bernoulli dà subito i momenti, senza calcolare serie: per la linearità E[Y] = np, e per l'indipendenza Var(Y) = npq.

Il programma verifica queste formule nel modo più diretto: costruisce la PMF con le frazioni esatte (per il quiz, n = 20 e p = 1/4), controlla che sommi a 1, calcola media e varianza dalla definizione, e confronta la probabilità di arrivare alla soglia con scipy.stats.binom e con un milione di quiz simulati. In coda cerca anche la soglia oltre la quale indovinare costa meno di una probabilità su mille.

from fractions import Fraction as F
from math import comb
import numpy as np
from scipy.stats import binom

# Un quiz di 20 domande a risposta multipla con 4 opzioni, risposte a caso: X = numero di risposte esatte
n, p = 20, F(1, 4)
q = 1 - p
pmf = [comb(n, k) * p**k * q**(n - k) for k in range(n + 1)]     # frazioni esatte
assert sum(pmf) == 1

media = sum(k * pk for k, pk in enumerate(pmf))
secondo = sum(k * k * pk for k, pk in enumerate(pmf))
var = secondo - media**2
print("E[X] =", media, "(np =", n * p, ")  Var(X) =", var, "(npq =", n * p * q, ")")
print("moda:", max(range(n + 1), key=lambda k: pmf[k]), " P(X=5) =", f"{float(pmf[5]):.4f}")

soglia = 10                                  # 'sufficiente' se almeno 10 su 20
esatto = sum(pmf[soglia:])
print(f"P(X >= {soglia}) esatto = {float(esatto):.6f}")
print(f"P(X >= {soglia}) scipy  = {binom.sf(soglia - 1, n, 0.25):.6f}")

N = 10**6
rng = np.random.default_rng(11)
campione = rng.binomial(n, 0.25, size=N)     # N studenti che tirano a indovinare
freq = np.mean(campione >= soglia)
pe = float(esatto)
tol = 3 * (pe * (1 - pe) / N) ** 0.5
print(f"Monte Carlo N={N}: {freq:.6f}; tolleranza 3 sigma {tol:.6f}; scarto {abs(freq - pe):.6f}")
print(f"media campionaria {campione.mean():.4f}, varianza campionaria {campione.var():.4f}")
assert abs(freq - pe) < tol

# Soglia sotto la quale indovinare "per caso" ha probabilita' inferiore a 1 su 1000
k_min = next(k for k in range(n + 1) if sum(pmf[k:]) < F(1, 1000))
print(f"soglia: P(X >= {k_min}) = {float(sum(pmf[k_min:])):.5f} < 0.001, mentre P(X >= {k_min - 1}) = {float(sum(pmf[k_min - 1:])):.5f}")

# Output:
# E[X] = 5 (np = 5 )  Var(X) = 15/4 (npq = 15/4 )
# moda: 5  P(X=5) = 0.2023
# P(X >= 10) esatto = 0.013864
# P(X >= 10) scipy  = 0.013864
# Monte Carlo N=1000000: 0.013781; tolleranza 3 sigma 0.000351; scarto 0.000083
# media campionaria 5.0024, varianza campionaria 3.7576
# soglia: P(X >= 12) = 0.00094 < 0.001, mentre P(X >= 11) = 0.00394

La media 5 e la varianza 15/4 escono esatte dalla definizione e coincidono con np e npq, quindi la dimostrazione per somma di Bernoulli e il calcolo diretto dicono la stessa cosa. La probabilità di prendere almeno 10 risposte esatte è 0,013864 per entrambi i calcoli, esatto e di SciPy, e la simulazione dà 0,013781, con scarto 0,000083 contro una tolleranza di 0,000351 (tre errori standard).

Il risultato pratico è netto: tirare a indovinare porta la sufficienza in meno di un caso su settanta, e con 12 risposte esatte la probabilità scende sotto l'uno per mille (0,00094). Nota anche la moda, il valore più probabile, che è 5 e vale solo il 20%: la distribuzione è molto più larga di un picco, la deviazione standard è circa 1,94. La varianza campionaria 3,758 è vicina a 15/4 = 3,75 come previsto.

Poisson: la binomiale quando il successo è raro

Quando n è grande e p è piccolo, calcolare C(n, k) diventa scomodo e, per molti fenomeni, nemmeno n e p separati sono noti: si conosce solo il numero medio di eventi, λ = np. La distribuzione di Poisson di parametro λ ha PMF P(Z = k) = e^−λ λ^k / k! per k = 0, 1, 2, …, e ha media e varianza entrambe uguali a λ.

È il limite della Bin(n, λ/n) per n che tende a infinito: fissato k, C(n, k) (λ/n)^k tende a λ^k/k!, e (1 − λ/n)^n−k tende a e^−λ. L'uguaglianza fra media e varianza è anche un test diagnostico: dati di conteggio con varianza molto maggiore della media non sono Poisson, e serve un modello diverso.

«Tende a» non dice quanto vicina sia la Poisson per un n concreto. La misura naturale è la distanza di variazione totale fra due distribuzioni, metà della somma delle differenze assolute delle PMF: è la massima differenza che le due leggi possono dare alla probabilità di un evento. Per una somma di Bernoulli indipendenti con probabilità p_i, la disuguaglianza di Le Cam (1960) la maggiora con Σ p_i², che per n prove uguali con p = λ/n vale λ²/n.

Il programma calcola la distanza esatta con λ = 3 per quattro valori di n e verifica la maggiorazione; poi applica la Poisson a un caso d'uso ipotetico, costruito per l'esempio: una compagnia vende 105 biglietti per 100 posti, e ogni passeggero si presenta con probabilità 0,95 indipendentemente dagli altri.

import numpy as np
from scipy.stats import binom, poisson

lam = 3.0
print(" n      TV(Bin(n,lam/n), Poisson(lam))   maggiorante lam^2/n   rapporto")
for n in (10, 100, 1000, 10000):
    k = np.arange(0, 200)
    tv = 0.5 * np.abs(binom.pmf(k, n, lam / n) - poisson.pmf(k, lam)).sum()
    bound = lam**2 / n                               # disuguaglianza di Le Cam: sum p_i^2 = n (lam/n)^2
    assert tv <= bound
    print(f"{n:6d}  {tv:.6f}                       {bound:.6f}            {tv / bound:.3f}")

# Un caso d'uso (ipotetico): 105 biglietti venduti per 100 posti, ogni passeggero si presenta con prob. 0,95
vendute, posti, p_assente = 105, 100, 0.05
assenti = binom(vendute, p_assente)
esatto = assenti.cdf(4)                              # overbooking <=> assenti <= 4  (presenti >= 101)
approx = poisson.cdf(4, vendute * p_assente)         # Poisson con lam = 5.25 sul numero di ASSENTI
print(f"P(overbooking) esatto = {esatto:.5f}  Poisson(5,25) = {approx:.5f}  differenza {abs(esatto - approx):.5f}")
print(f"maggiorante di Le Cam sulla distanza: {vendute * p_assente**2:.5f}")

rng = np.random.default_rng(5)
N = 10**6
presenti = rng.binomial(vendute, 1 - p_assente, size=N)
freq = np.mean(presenti > posti)
tol = 3 * (esatto * (1 - esatto) / N) ** 0.5
print(f"Monte Carlo N={N}: {freq:.5f}; tolleranza {tol:.5f}; scarto {abs(freq - esatto):.5f}")
assert abs(freq - esatto) < tol

# Output:
#  n      TV(Bin(n,lam/n), Poisson(lam))   maggiorante lam^2/n   rapporto
#     10  0.086409                       0.900000            0.096
#    100  0.007607                       0.090000            0.085
#   1000  0.000752                       0.009000            0.084
#  10000  0.000075                       0.000900            0.083
# P(overbooking) esatto = 0.39243  Poisson(5,25) = 0.39777  differenza 0.00534
# maggiorante di Le Cam sulla distanza: 0.26250
# Monte Carlo N=1000000: 0.39187; tolleranza 0.00146; scarto 0.00057

La distanza esatta scende come 1/n, da 0,086 con n = 10 a 0,000075 con n = 10.000, ed è sempre sotto il maggiorante: nei quattro casi misurati la maggiorazione di Le Cam è tra dieci e dodici volte più larga del valore vero (rapporto 0,083-0,096). Serve quindi a garantire un errore massimo, non a stimarlo.

Nel caso dei biglietti si applica la Poisson al numero di assenti, non a quello dei presenti: gli assenti sono l'evento raro, con media 105 · 0,05 = 5,25, e c'è sovraprenotazione quando gli assenti sono al più 4. La probabilità esatta è 0,39243, la Poisson(5,25) dà 0,39777, e la simulazione 0,39187. L'approssimazione sbaglia di 0,005 mentre il maggiorante di Le Cam, 0,2625, è cinquanta volte più largo dello scarto osservato sulla probabilità richiesta: la garanzia è prudente, non precisa. Il numero è ipotetico e illustra il metodo, non descrive una compagnia reale.

Geometrica: aspettare il primo successo

Se ripeti prove indipendenti finché non ottieni il primo successo, il numero X di prove ha distribuzione geometrica: P(X = k) = q^k−1p, con media 1/p e varianza q/p². Qui «tentativi» conta anche l'ultimo, il successo; altri testi contano solo i fallimenti e partono da zero, per cui la media diventa q/p: controlla sempre quale convenzione usa il testo o la libreria.

La proprietà caratteristica è l'assenza di memoria: P(X > m + n | X > m) = P(X > n), perché la probabilità di superare m prove senza successo è q^m e il rapporto q^m+n/q^m è qⁿ. L'esempio è il dado: lanci finché non esce un 6, quindi p = 1/6. Il programma verifica l'uguaglianza con le frazioni per tre coppie (m, n), somma la serie della media, e simula un milione di sequenze.

from fractions import Fraction as F
import numpy as np

# Lancio un dado finche' non esce il 6: X = numero di lanci (X = 1, 2, 3, ...)
p = F(1, 6); q = 1 - p
coda = lambda m: q**m                                # P(X > m): i primi m lanci falliscono
for m, n in ((2, 3), (5, 4), (10, 1)):
    cond = coda(m + n) / coda(m)                     # P(X > m+n | X > m)
    print(f"m={m:2d} n={n}: P(X>m+n | X>m) = {cond}  =  P(X>n) = {coda(n)}  -> {cond == coda(n)}")

media = 1 / p; var = q / p**2
somma = sum(k * q**(k - 1) * p for k in range(1, 400))   # somma parziale della serie E[X]
print(f"E[X] = 1/p = {media}, Var(X) = q/p^2 = {var}; serie troncata a 399 termini: {float(somma):.6f}")

rng = np.random.default_rng(3)
N = 10**6
x = rng.geometric(1/6, size=N)                       # numpy conta i tentativi, incluso il successo
sigma_media = float(var) ** 0.5 / N ** 0.5
print(f"Monte Carlo: media {x.mean():.4f} (atteso {float(media):.0f} +/- {3*sigma_media:.4f}), varianza {x.var():.2f} (atteso {float(var):.0f})")
assert abs(x.mean() - float(media)) < 3 * sigma_media
# Condizionata: tra le prove con X > 5, la distribuzione di X-5 e' di nuovo geometrica
resto = x[x > 5] - 5
print(f"dopo 5 fallimenti: {len(resto)} prove, media dei lanci ancora necessari {resto.mean():.4f}")

# Output:
# m= 2 n=3: P(X>m+n | X>m) = 125/216  =  P(X>n) = 125/216  -> True
# m= 5 n=4: P(X>m+n | X>m) = 625/1296  =  P(X>n) = 625/1296  -> True
# m=10 n=1: P(X>m+n | X>m) = 5/6  =  P(X>n) = 5/6  -> True
# E[X] = 1/p = 6, Var(X) = q/p^2 = 30; serie troncata a 399 termini: 6.000000
# Monte Carlo: media 6.0037 (atteso 6 +/- 0.0164), varianza 29.97 (atteso 30)
# dopo 5 fallimenti: 402794 prove, media dei lanci ancora necessari 5.9900

Le tre uguaglianze sono esatte (per esempio 625/1296 da entrambe le parti), la serie troncata a 399 termini dà 6,000000, e la simulazione trova media 6,0037 contro un atteso di 6 ± 0,0164 e varianza 29,97 contro 30.

L'ultima riga mostra la memoria che manca: tra le sequenze in cui i primi cinque lanci hanno fallito, servono ancora in media 5,99 lanci, quasi lo stesso 6 di una sequenza nuova, e non uno in meno come suggerirebbe l'idea che il 6 sia «ormai dovuto». È la correzione al ragionamento del giocatore d'azzardo: il dado non ricorda i lanci passati. L'analogo continuo di questa proprietà è la distribuzione esponenziale, che incontri nel prossimo articolo.

Quando il modello non regge

Le formule di questo articolo richiedono ipotesi, e le ipotesi sono il punto debole di ogni applicazione. La binomiale presuppone prove indipendenti e con la stessa probabilità: se i passeggeri viaggiano in gruppo, o se i quesiti del quiz non hanno la stessa difficoltà, la varianza reale è diversa da npq e la distribuzione non è più binomiale. La Poisson presuppone eventi rari e indipendenti: conteggi con varianza molto maggiore della media indicano raggruppamento (eventi a grappolo) e richiedono famiglie più flessibili, come la binomiale negativa, che non tratto qui.

Questo articolo lascia fuori anche la funzione generatrice dei momenti, l'ipergeometrica (estrazioni senza reinserimento) e la stima dei parametri da dati: sono argomenti di inferenza, dove la stima di p da una frequenza osservata riprende il calcolo di media e varianza fatto qui. Per il comportamento della binomiale quando n cresce con p fisso, il risultato giusto è il teorema centrale del limite, nell'articolo 05.

Domande frequenti

Che differenza c'è tra binomiale e Poisson?

La binomiale conta i successi in un numero fisso n di prove, con un limite superiore uguale a n. La Poisson non ha limite superiore e dipende da un solo parametro, la media; approssima bene la binomiale quando n è grande e p è piccolo.

Perché in una Poisson media e varianza coincidono?

Perché la Poisson è il limite di una binomiale con p che tende a zero: la varianza npq tende a np = λ, mentre la media resta np = λ. Se in dati reali la varianza supera di molto la media, il modello Poisson non è adeguato.

La geometrica parte da 0 o da 1?

Dipende dalla convenzione. Contando i tentativi, incluso il successo, parte da 1 e ha media 1/p; contando i soli fallimenti parte da 0 e ha media (1 − p)/p. NumPy usa la prima, come il programma sopra.

Che cosa sai fare adesso

Sai descrivere una variabile discreta con la sua PMF e calcolarne media e varianza dalla definizione, riconoscere quando un conteggio è binomiale e quando è ragionevole approssimarlo con una Poisson, e quantificare l'errore di quell'approssimazione invece di fidarti di una regola a occhio. Sai controllare una formula in due modi, con le frazioni esatte e con una simulazione a seme fisso e tolleranza dichiarata. L'articolo successivo sposta tutto sul continuo: densità, funzione di ripartizione, e le distribuzioni uniforme, esponenziale e normale.

Nella stessa serie trovi anche Analisi dati esplorativa (EDA) con Python. L'elenco completo è nella pagina della serie, Probabilità e Statistica.