03 - Variabili aleatorie continue: Uniforme, Esponenziale, Normale

Aspetti un tram che passa ogni dieci minuti, ma non sai quando è passato l'ultimo. Qual è la probabilità di aspettare esattamente tre minuti? Zero. Qual è la probabilità di aspettare meno di tre minuti? Il 30%. Le due risposte sembrano contraddirsi finché non capisci come funziona una variabile che può assumere qualsiasi valore in un intervallo: non si assegna probabilità ai punti, la si assegna agli intervalli, calcolandola come un'area.

Questo articolo costruisce lo strumento (densità e funzione di ripartizione), lo applica alle tre distribuzioni più usate, uniforme, esponenziale e normale, e mostra come cambia la densità quando una variabile viene trasformata. Come nell'articolo sulle variabili discrete, ogni risultato è calcolato due volte, con la formula e con una simulazione a seme fisso, e la tolleranza della simulazione è dichiarata prima di leggerne l'esito. Ti serve saper calcolare un integrale elementare.

Densità e funzione di ripartizione

Una variabile aleatoria è continua se esiste una densità f(x) ≥ 0, con area totale 1, tale che la probabilità di cadere in un intervallo sia l'area sotto la densità: P(a < X ≤ b) = ∫_a^b f(x) dx.

La funzione di ripartizione (in inglese CDF) F(x) = P(X ≤ x) accumula quell'area da sinistra, quindi P(a < X ≤ b) = F(b) − F(a), e dove F è derivabile la densità è la sua derivata.

Due conseguenze vanno capite bene. Primo: la densità non è una probabilità, può superare 1 (una densità uniforme su un intervallo lungo 0,1 vale 10) ed è solo la probabilità per unità di lunghezza. Secondo: P(X = x) = 0 per ogni singolo valore x, perché l'area di un intervallo che si restringe a un punto tende a zero; questo non rende il valore impossibile, perché uno dei valori deve pur uscire.

Valore atteso e varianza si definiscono come per le variabili discrete, con l'integrale al posto della somma: E[X] = ∫ x f(x) dx e Var(X) = ∫ (x − E[X])² f(x) dx. Valgono le stesse regole: la media di una somma è la somma delle medie, la varianza si somma per variabili indipendenti.

Il tram è il caso più semplice, la distribuzione uniforme su [a, b]: densità costante 1/(b − a) in quell'intervallo e zero fuori, media (a + b)/2 e varianza (b − a)²/12. Il programma integra numericamente la densità del tram con scipy.integrate.quad per controllare area, media e varianza, calcola le probabilità con la CDF, mostra come la probabilità di un intervallo di ampiezza h attorno a 3 minuti si riduca con h, e simula un milione di arrivi.

from fractions import Fraction as F
import numpy as np
from scipy.integrate import quad

# Aspetti il tram: l'ultimo e' passato a un istante che non conosci, il prossimo arriva tra 0 e 10 minuti.
a, b = 0.0, 10.0
f = lambda x: 1 / (b - a) if a <= x <= b else 0.0      # densita'
Fc = lambda x: min(max((x - a) / (b - a), 0.0), 1.0)    # funzione di ripartizione

area, _ = quad(f, a, b)
media, _ = quad(lambda x: x * f(x), a, b)
var, _ = quad(lambda x: (x - media) ** 2 * f(x), a, b)
print(f"area della densita' = {area:.6f}  media = {media:.6f}  varianza = {var:.6f}  (formule: {(a+b)/2}, {(b-a)**2/12:.6f})")
print(f"P(attesa < 3 min) = {Fc(3):.4f}   P(2 < attesa <= 7) = {Fc(7) - Fc(2):.4f}")
print("P(attesa = 3 min esattamente) = 0, ma P(3 <= attesa <= 3 + h) vale h/10 =", [f"{h}: {Fc(3 + h) - Fc(3):.4f}" for h in (1, 0.1, 0.001)])

N = 10**6
rng = np.random.default_rng(21)
x = rng.uniform(a, b, size=N)
for nome, vera, oss in (("P(<3)", 0.3, np.mean(x < 3)), ("P(2<..<=7)", 0.5, np.mean((x > 2) & (x <= 7)))):
    tol = 3 * (vera * (1 - vera) / N) ** 0.5
    print(f"{nome}: Monte Carlo {oss:.5f} atteso {vera} +/- {tol:.5f}")
    assert abs(oss - vera) < tol
print(f"media campionaria {x.mean():.4f}, varianza campionaria {x.var():.4f}")

# Output:
# area della densita' = 1.000000  media = 5.000000  varianza = 8.333333  (formule: 5.0, 8.333333)
# P(attesa < 3 min) = 0.3000   P(2 < attesa <= 7) = 0.5000
# P(attesa = 3 min esattamente) = 0, ma P(3 <= attesa <= 3 + h) vale h/10 = ['1: 0.1000', '0.1: 0.0100', '0.001: 0.0001']
# P(<3): Monte Carlo 0.29924 atteso 0.3 +/- 0.00137
# P(2<..<=7): Monte Carlo 0.49996 atteso 0.5 +/- 0.00150
# media campionaria 5.0026, varianza campionaria 8.3285

L'area vale 1, la media 5 minuti e la varianza 8,33 minuti al quadrato, uguali alle formule. L'ultima serie di numeri della terza riga è la risposta al paradosso dell'inizio: l'intervallo di ampiezza h attorno a 3 minuti ha probabilità h/10, che tende a zero con h, quindi «esattamente 3 minuti» pesa zero, mentre «meno di 3» pesa 0,3.

Le simulazioni confermano entro tolleranza (scarti 0,00076 e 0,00004 contro tolleranze di 0,0014 e 0,0015). Un'avvertenza sul modello: l'uniforme descrive il tram solo se arriva puntuale e tu arrivi in un istante casuale e indipendente dall'orario; se arrivi per prendere proprio quel tram, la distribuzione non è più questa.

Esponenziale: il tempo di attesa senza memoria

La distribuzione esponenziale di parametro λ > 0 ha densità f(x) = λe^−λx per x ≥ 0, funzione di ripartizione F(x) = 1 − e^−λx, media 1/λ e varianza 1/λ². Modella un tempo di attesa.

La funzione S(t) = P(X > t) = e^−λt, la probabilità di sopravvivere oltre t, ha una proprietà notevole: P(X > s + t | X > s) = e^−λ(s+t)/e^−λs = e^−λt = P(X > t).

Chi è sopravvissuto fino a s ha davanti lo stesso futuro di chi parte da zero: è l'assenza di memoria, l'analogo continuo di quanto hai visto per la geometrica. È una proprietà forte e raramente vera nei sistemi reali: un componente che si usura non è senza memoria. L'esponenziale è il modello giusto solo per guasti che arrivano a caso, senza invecchiamento.

L'esempio è un componente con vita media di 1000 ore (λ = 0,001 per ora; ipotetico, scelto per rendere i numeri leggibili). La probabilità che duri oltre 1500 ore è e^−1,5, e quella di durare oltre 1500 ore sapendo che ha già superato le 1000 è e^−0,5, identica alla probabilità di durare 500 ore da nuovo.

Il programma campiona per inversione: se U è uniforme su (0, 1], allora X = −ln(U)/λ è esponenziale, perché inverte la CDF. L'istruzione 1 - rng.random(N) serve a evitare il logaritmo di zero, dato che random() restituisce valori in [0, 1). Il test di Kolmogorov-Smirnov misura la massima distanza fra la CDF empirica del campione e quella teorica.

import numpy as np
from math import exp, log
from scipy.stats import kstest, expon, uniform

# Vita di un componente: esponenziale con media 1000 ore, quindi lambda = 1/1000
lam = 1 / 1000
S = lambda t: exp(-lam * t)                      # P(X > t), funzione di sopravvivenza
print(f"P(X > 1500) = e^-1.5 = {S(1500):.5f}")
print(f"P(X > 1500 | X > 1000) = {S(1500) / S(1000):.5f}  e  P(X > 500) = {S(500):.5f}")

# Campionamento per inversione: se U e' uniforme su (0,1], X = -ln(U)/lambda e' esponenziale
N = 10**6
rng = np.random.default_rng(9)
u = 1 - rng.random(N)                            # rng.random() sta in [0,1): 1-u evita log(0)
x = -np.log(u) / lam
print(f"media {x.mean():.2f} (attesa 1000, errore standard {1000 / N**0.5:.2f}); varianza {x.var():.0f} (attesa 1000000)")

oss = np.mean(x > 1500)
tol = 3 * (S(1500) * (1 - S(1500)) / N) ** 0.5
print(f"frequenza X>1500: {oss:.5f} (atteso {S(1500):.5f} +/- {tol:.5f})")
assert abs(oss - S(1500)) < tol

dopo = x[x > 1000]                               # i componenti ancora vivi a 1000 ore
oss_c = np.mean(dopo > 1500)
tol_c = 3 * (S(500) * (1 - S(500)) / len(dopo)) ** 0.5
print(f"tra i {len(dopo)} componenti vivi a 1000 ore, quota viva a 1500: {oss_c:.5f} (atteso {S(500):.5f} +/- {tol_c:.5f})")
assert abs(oss_c - S(500)) < tol_c

# Il test di Kolmogorov-Smirnov confronta la funzione di ripartizione empirica con quella teorica
print(f"KS contro Exp(1/1000): D = {kstest(x, expon(scale=1000).cdf).statistic:.5f}")
# Trasformata integrale di probabilita': F(X) e' uniforme su (0,1) quando F e' la ripartizione di X
print(f"KS di F(X) contro l'uniforme: D = {kstest(1 - np.exp(-lam * x), uniform.cdf).statistic:.5f}")

# Output:
# P(X > 1500) = e^-1.5 = 0.22313
# P(X > 1500 | X > 1000) = 0.60653  e  P(X > 500) = 0.60653
# media 999.88 (attesa 1000, errore standard 1.00); varianza 1003337 (attesa 1000000)
# frequenza X>1500: 0.22330 (atteso 0.22313 +/- 0.00125)
# tra i 366897 componenti vivi a 1000 ore, quota viva a 1500: 0.60861 (atteso 0.60653 +/- 0.00242)
# KS contro Exp(1/1000): D = 0.00108
# KS di F(X) contro l'uniforme: D = 0.00108

I due valori teorici, 0,60653, coincidono come previsto, e le frequenze simulate sono 0,22330 contro 0,22313 e 0,60861 contro 0,60653. Quest'ultimo scarto è 2,6 errori standard, il più alto fra i quaranta semi che ho provato per controllare che non ci fosse una distorsione (media degli scarti −0,07 in unità di errore standard, deviazione 0,96): è dunque una oscillazione normale, e un esempio di perché la tolleranza è fissata a tre errori e non a uno.

La varianza campionaria, 1.003.337, differisce dello 0,33% da 10⁶, entro il suo errore (la varianza stimata di una variabile così asimmetrica converge lentamente, con errore standard di circa 2.800). Le due distanze di Kolmogorov-Smirnov, 0,00108, coincidono perché trasformare il campione con la sua CDF è una trasformazione crescente che non cambia l'ordine dei valori: se F è la CDF di X, F(X) è uniforme su (0, 1), ed è il fondamento del campionamento per inversione.

Normale: la distribuzione dell'errore di misura

La distribuzione normale (o gaussiana) di media μ e varianza σ² ha densità f(x) = (1/(σ√(2π))) exp(−(x − μ)²/(2σ²)): una campana simmetrica centrata in μ, con ampiezza regolata da σ. La normale standard è quella con μ = 0 e σ = 1; la sua CDF si indica Φ e non ha una forma elementare, quindi si calcola numericamente (con la funzione degli errori, erf, o con una libreria).

Valgono due proprietà che usi di continuo: se X è normale, anche aX + b lo è, con media aμ + b e varianza a²σ²; e la somma di normali indipendenti è normale, con medie e varianze che si sommano. Per questo basta conoscere la standard: ogni normale si riconduce a lei sottraendo la media e dividendo per la deviazione standard.

La probabilità di stare entro k deviazioni standard dalla media non dipende da μ né da σ: vale erf(k/√2). Il programma la calcola così, la confronta con SciPy e con un milione di estrazioni, e verifica che l'area della densità sia 1.

Poi applica tutto a un controllo di qualità ipotetico: pezzi con diametro normale di media 10,00 mm e deviazione standard 0,02 mm, tolleranza ±0,03 mm. Per simmetria la quota di pezzi fuori tolleranza è 2(1 − Φ(1,5)). Infine calcola l'intervallo centrale al 99% e verifica la regola sulla somma di due normali indipendenti.

import numpy as np
from math import erf, sqrt
from scipy.integrate import quad
from scipy.stats import norm

phi = lambda t: np.exp(-t * t / 2) / sqrt(2 * np.pi)       # densita' della normale standard
area, _ = quad(phi, -np.inf, np.inf)
print(f"integrale della densita' standard = {area:.10f}")

# Probabilita' entro k deviazioni standard: P(|Z|<k) = erf(k/sqrt(2)), senza tabelle
N = 10**6
rng = np.random.default_rng(13)
z = rng.standard_normal(N)
for k in (1, 2, 3):
    esatta = erf(k / sqrt(2))
    oss = np.mean(np.abs(z) < k)
    tol = 3 * (esatta * (1 - esatta) / N) ** 0.5
    print(f"|Z|<{k}: erf {esatta:.6f}  scipy {norm.cdf(k) - norm.cdf(-k):.6f}  Monte Carlo {oss:.6f} (+/- {tol:.6f})")
    assert abs(oss - esatta) < tol

# Controllo qualita' (caso ipotetico): diametro N(10.00, 0.02^2) mm, tolleranza 10.00 +/- 0.03 mm
mu, sigma, tol_mm = 10.00, 0.02, 0.03
fuori = 2 * (1 - norm.cdf(tol_mm / sigma))                  # simmetria: due code uguali
x = rng.normal(mu, sigma, size=N)
oss = np.mean(np.abs(x - mu) > tol_mm)
tol = 3 * (fuori * (1 - fuori) / N) ** 0.5
print(f"pezzi fuori tolleranza: esatto {fuori:.5f}, Monte Carlo {oss:.5f} (+/- {tol:.5f})")
assert abs(oss - fuori) < tol
print(f"intervallo centrale al 99%: media +/- {norm.ppf(0.995):.4f} sigma = [{mu - norm.ppf(0.995)*sigma:.4f}, {mu + norm.ppf(0.995)*sigma:.4f}] mm")

# Somma di normali indipendenti: N(mu1+mu2, s1^2+s2^2)
y = rng.normal(5, 0.015, size=N)
s = x + y
print(f"somma: media {s.mean():.4f} (15), deviazione {s.std():.5f} (attesa {sqrt(0.02**2 + 0.015**2):.5f})")

# Output:
# integrale della densita' standard = 1.0000000000
# |Z|<1: erf 0.682689  scipy 0.682689  Monte Carlo 0.682402 (+/- 0.001396)
# |Z|<2: erf 0.954500  scipy 0.954500  Monte Carlo 0.954419 (+/- 0.000625)
# |Z|<3: erf 0.997300  scipy 0.997300  Monte Carlo 0.997275 (+/- 0.000156)
# pezzi fuori tolleranza: esatto 0.13361, Monte Carlo 0.13355 (+/- 0.00102)
# intervallo centrale al 99%: media +/- 2.5758 sigma = [9.9485, 10.0515] mm
# somma: media 15.0000 (15), deviazione 0.02500 (attesa 0.02500)

Le probabilità entro uno, due e tre sigma sono 0,682689, 0,954500 e 0,997300, e la simulazione le conferma entro tolleranza. Nel caso dei pezzi, il 13,36% cade fuori tolleranza: il progettista legge in questo numero che una tolleranza di 1,5 sigma è troppo stretta, e che per scartare meno dell'1% occorre una tolleranza di almeno 2,58 sigma (l'intervallo al 99% va da 9,9485 a 10,0515 mm).

Sommando un'altra variabile normale indipendente di media 5 e deviazione 0,015, la deviazione standard del totale è √(0,02² + 0,015²) = 0,025, non 0,035: le deviazioni standard non si sommano, si sommano le varianze. La normale è un buon modello quando l'errore è somma di molti piccoli contributi indipendenti, come spiega l'articolo sul teorema centrale del limite; dove ci sono eventi estremi molto più frequenti di così (code pesanti), la sua stima delle probabilità rare è troppo ottimista.

Trasformazioni: come cambia la densità

Se Y = g(X) e g è monotona e derivabile, la densità di Y è f_Y(y) = f_X(g⁻¹(y)) · |dg⁻¹/dy|: il fattore con la derivata è lo «Jacobiano» in una variabile, e compensa il fatto che g allunga o accorcia l'asse. Se g non è monotona si sommano i contributi dei rami.

L'esempio classico è Y = X² con X normale standard: ogni valore positivo y proviene da due valori di X, +√y e −√y, ciascuno con peso φ(√y)/(2√y), da cui f_Y(y) = e^−y/2/√(2πy).

Questa è la distribuzione chi-quadro con un grado di libertà. Il programma controlla la formula in tre modi indipendenti, che devono dare lo stesso numero: l'integrale numerico della densità, la CDF ricavata direttamente come 2Φ(√y) − 1 (perché Y ≤ y equivale a |X| ≤ √y), e la CDF chi-quadro di SciPy; e li confronta con la simulazione.

import numpy as np
from math import sqrt, pi, exp
from scipy.integrate import quad
from scipy.stats import norm, chi2

# Y = X^2 con X normale standard. X^2 non e' monotona: si somma sui due rami X = +sqrt(y), X = -sqrt(y).
phi = lambda t: exp(-t * t / 2) / sqrt(2 * pi)
f_Y = lambda y: 2 * phi(sqrt(y)) / (2 * sqrt(y))        # f_X(sqrt y) |d sqrt(y)/dy| per ciascun ramo
print("f_Y(y) = exp(-y/2) / sqrt(2 pi y)?", all(abs(f_Y(y) - exp(-y / 2) / sqrt(2 * pi * y)) < 1e-15 for y in (0.1, 1, 4)))

N = 10**6
rng = np.random.default_rng(17)
y = rng.standard_normal(N) ** 2
print(f"{'intervallo':>12} {'integrale di f_Y':>17} {'2Phi(sqrt)-1':>14} {'chi2(1)':>9} {'Monte Carlo':>12}")
for a, b in ((0.0, 1.0), (1.0, 4.0), (4.0, 9.0)):
    integr, _ = quad(f_Y, a, b)
    dalla_cdf = (2 * norm.cdf(sqrt(b)) - 1) - (2 * norm.cdf(sqrt(a)) - 1)
    chi = chi2.cdf(b, 1) - chi2.cdf(a, 1)
    oss = np.mean((y > a) & (y <= b))
    tol = 3 * (dalla_cdf * (1 - dalla_cdf) / N) ** 0.5
    assert abs(integr - dalla_cdf) < 1e-9 and abs(chi - dalla_cdf) < 1e-12 and abs(oss - dalla_cdf) < tol
    print(f"({a:.0f}, {b:.0f}] {integr:17.6f} {dalla_cdf:14.6f} {chi:9.6f} {oss:12.6f}")
print(f"E[Y] = {y.mean():.4f} (esatto 1), Var(Y) = {y.var():.3f} (esatto 2)")

# Output:
# f_Y(y) = exp(-y/2) / sqrt(2 pi y)? True
#   intervallo  integrale di f_Y   2Phi(sqrt)-1   chi2(1)  Monte Carlo
# (0, 1]          0.682689       0.682689  0.682689     0.682203
# (1, 4]          0.271810       0.271810  0.271810     0.272241
# (4, 9]          0.042800       0.042800  0.042800     0.042816
# E[Y] = 1.0011 (esatto 1), Var(Y) = 2.000 (esatto 2)

I tre calcoli concordano alla sesta cifra decimale in ogni intervallo (0,682689 in (0, 1], 0,271810 in (1, 4], 0,042800 in (4, 9]) e le asserzioni verificano che concordino anche a 10⁻⁹; la simulazione è entro tolleranza. Il primo valore coincide con la probabilità di stare entro una deviazione standard della normale, e non per caso: X² ≤ 1 equivale a |X| ≤ 1.

La media di Y è 1, la varianza 2, valori che si ricavano dai momenti della normale (E[X²] = 1, E[X⁴] = 3, quindi Var(X²) = 3 − 1 = 2). L'errore tipico è dimenticare il ramo negativo e ottenere una densità che integra a 1/2: il controllo dell'area, come nel primo programma, lo rivela subito.

Che cosa resta fuori

Qui mancano la funzione generatrice dei momenti e quella caratteristica, il processo di Poisson (di cui l'esponenziale descrive i tempi fra arrivi successivi), le famiglie derivate dalla normale come la t di Student, il campionamento della normale con il metodo di Box-Muller e le trasformazioni in più dimensioni con il determinante dello Jacobiano. Per le trasformazioni di vettori vai all'articolo 04, dove compare la normale multivariata. Per tornare al caso discreto e ai suoi legami con questi modelli, l'articolo 01 contiene gli assiomi da cui anche questi integrali discendono.

Domande frequenti

Perché la probabilità di un singolo valore è zero?

Perché la probabilità è l'area sotto la densità, e l'area di un intervallo di ampiezza zero è zero. Questo non significa che il valore non possa uscire: significa che, fra infiniti valori possibili, nessuno ha una probabilità positiva da solo.

Una densità può essere maggiore di 1?

Sì. La densità è una probabilità per unità di lunghezza: su un intervallo corto può superare 1, purché l'area totale resti 1. Una uniforme su un intervallo lungo 0,1 ha densità 10.

Che cosa vuol dire che l'esponenziale non ha memoria?

Che la probabilità di attendere ancora un tempo t non dipende da quanto hai già atteso. Vale per guasti puramente casuali, non per componenti che si usurano.

Che cosa sai fare adesso

Sai passare da densità a CDF e viceversa, calcolare la probabilità di un intervallo come un'area e capire perché un singolo punto pesa zero. Sai usare uniforme, esponenziale e normale, campionare per inversione, ridurre una normale qualsiasi alla standard e trovare la densità di una funzione di una variabile, controllando ogni risultato con un secondo metodo. Nel prossimo articolo si passa da una variabile a più variabili insieme: la distribuzione congiunta, la covarianza e la normale multivariata.

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