Distribuzione Gamma: 3 esercizi svolti con Python su probabilità, latenza e affidabilità

Se un servizio online impiega in media mezzo secondo a rispondere, significa che quasi tutte le richieste vengono elaborate entro mezzo secondo?

E se un sistema informatico deve attraversare più fasi operative prima di completare un’attività, come possiamo stimare il tempo complessivo necessario?

La media è un indicatore utile, ma non racconta tutta la storia, soprattutto quando i tempi di attesa presentano una distribuzione asimmetrica e alcuni valori particolarmente elevati allungano la coda verso destra.

Nei sistemi informatici, nelle consegne e nei processi industriali, conoscere il tempo medio non basta per valutare il rischio di ritardi o stabilire una soglia operativa ragionevole.

La distribuzione Gamma permette di descrivere molte variabili casuali continue che assumono valori positivi e presentano un’asimmetria verso destra. In questo articolo la utilizzeremo per affrontare tre problemi progressivi: la probabilità associata ai tempi di consegna, il superamento di una soglia di latenza in un servizio di inferenza AI e la durata complessiva di una sequenza di fasi operative.

Per ciascun esercizio partiremo dalle formule, verificheremo i risultati con Python e concluderemo con un’interpretazione dei numeri, perché un calcolo statisticamente corretto diventa utile soltanto quando sappiamo che cosa ci permette di decidere e quali ipotesi stiamo accettando.

La distribuzione Gamma: definizione e significato dei parametri

Una variabile casuale continua [math]X[/math] segue una distribuzione Gamma con parametro di forma [math]\alpha[/math] e parametro di tasso [math]\beta[/math] quando scriviamo

[math]\displaystyle \begin{aligned}
X\sim\operatorname{Gamma}(\alpha,\beta)
\end{aligned}[/math]

con [math]\alpha>0[/math] e [math]\beta>0[/math]. Adottiamo qui la parametrizzazione basata sul tasso, chiamato rate nella letteratura statistica, anziché quella basata sulla scala.

La funzione di densità di probabilità è

[math]\displaystyle \begin{aligned}
f(x)=\frac{\beta^\alpha}{\Gamma(\alpha)}
x^{\alpha-1}e^{-\beta x},
\qquad x>0
\end{aligned}[/math]

dove [math]\Gamma(\alpha)[/math] è la funzione Gamma, che generalizza il fattoriale ai valori reali positivi. Quando [math]\alpha[/math] è un intero positivo, vale la relazione

[math]\displaystyle \begin{aligned}
\Gamma(\alpha)=(\alpha-1)!
\end{aligned}[/math]

Per esempio, se [math]\alpha=3[/math], allora [math]\Gamma(3)=2!=2[/math], mentre per [math]\alpha=4[/math] otteniamo [math]\Gamma(4)=3!=6[/math].

Il parametro di forma [math]\alpha[/math] controlla principalmente la geometria della distribuzione. Quando [math]\alpha<1[/math], la densità cresce senza limite avvicinandosi a zero; quando [math]\alpha=1[/math], la Gamma coincide con una distribuzione esponenziale; quando [math]\alpha>1[/math], la densità parte da zero, raggiunge un massimo e poi decresce. All’aumentare di [math]\alpha[/math], la distribuzione tende a diventare meno asimmetrica rispetto alla propria media, anche se non esiste una soglia universale oltre la quale l’approssimazione normale sia sempre adeguata.

Il parametro di tasso [math]\beta[/math] determina invece la scala dei valori osservati. A parità di [math]\alpha[/math], un tasso più elevato concentra la distribuzione su valori più piccoli, mentre un tasso più basso la allunga verso destra. Se preferiamo utilizzare il parametro di scala [math]\theta[/math], la relazione da ricordare è

[math]\displaystyle \begin{aligned}
\theta=\frac{1}{\beta}
\end{aligned}[/math]

Questa distinzione è importante quando si passa dalle formule a Python, perché SciPy utilizza il parametro di scala.

Media, varianza e moda

Per una variabile [math]X\sim\operatorname{Gamma}(\alpha,\beta)[/math], il valore atteso è

[math]\displaystyle \begin{aligned}
E[X]=\frac{\alpha}{\beta}
\end{aligned}[/math]

mentre la varianza è

[math]\displaystyle \begin{aligned}
\operatorname{Var}(X)=\frac{\alpha}{\beta^2}
\end{aligned}[/math]

La deviazione standard, che esprime la dispersione nella stessa unità di misura della variabile, vale

[math]\displaystyle \begin{aligned}
\sigma_X=\frac{\sqrt{\alpha}}{\beta}
\end{aligned}[/math]

Quando [math]\alpha>1[/math], la moda, cioè il punto in cui la densità raggiunge il massimo, è

[math]\displaystyle \begin{aligned}
\operatorname{Moda}(X)=\frac{\alpha-1}{\beta}
\end{aligned}[/math]

Per [math]\alpha=1[/math], la distribuzione esponenziale ha il massimo al bordo [math]x=0[/math]; per [math]0<\alpha<1[/math], la densità diverge avvicinandosi a zero e il massimo non è un punto interno del dominio.

Il confronto fra media e moda è particolarmente utile quando la distribuzione è asimmetrica. Per [math]\alpha>1[/math], la differenza è

[math]\displaystyle \begin{aligned}
E[X]-\operatorname{Moda}(X)=\frac{1}{\beta}
\end{aligned}[/math]

La media risente maggiormente dei valori elevati della coda destra, mentre la moda individua il punto di massima densità. Non bisogna quindi interpretare questi due indicatori come se descrivessero la stessa caratteristica dei dati.

Densità e probabilità non sono la stessa cosa

La funzione [math]f(x)[/math] descrive la densità della distribuzione, mentre la probabilità di osservare un valore all’interno di un intervallo si ottiene calcolando l’area sottesa alla curva.

La funzione di ripartizione, indicata con [math]F(x)[/math], è definita come

[math]\displaystyle \begin{aligned}
F(x)=P(X\leq x)
\end{aligned}[/math]

e permette di calcolare direttamente probabilità come

[math]\displaystyle \begin{aligned}
P(a<X<b)=F(b)-F(a)
\end{aligned}[/math]

Per una variabile casuale continua, la probabilità di osservare esattamente un valore prefissato è zero. Se la densità in corrispondenza di dieci minuti vale [math]0{,}0641[/math], per esempio, non significa che esista una probabilità del 6,41% di ottenere esattamente dieci minuti: il valore della densità deve essere integrato su un intervallo per ottenere una probabilità.

Perché la Gamma compare nei problemi reali

Una proprietà importante della distribuzione Gamma è che la somma di [math]n[/math] variabili esponenziali indipendenti, tutte con lo stesso tasso [math]\beta[/math], segue una distribuzione Gamma con forma [math]n[/math] e tasso [math]\beta[/math]:

[math]\displaystyle \begin{aligned}
X_i\sim\operatorname{Exp}(\beta)
\quad\Longrightarrow\quad
\sum_{i=1}^{n}X_i\sim\operatorname{Gamma}(n,\beta)
\end{aligned}[/math]

Questa proprietà descrive, per esempio, la durata complessiva di più fasi successive, purché le rispettive durate siano indipendenti e abbiano lo stesso tasso esponenziale.

La distribuzione Gamma ha anche un ruolo nell’inferenza bayesiana. Se il numero di eventi osservati segue un modello di Poisson e scegliamo una distribuzione a priori Gamma per il tasso degli eventi, la distribuzione a posteriori rimane Gamma. Con una parametrizzazione basata su forma e tasso, se la prior è [math]\operatorname{Gamma}(a,b)[/math], dopo aver osservato [math]N[/math] eventi durante un’esposizione totale [math]T[/math], la posterior è [math]\operatorname{Gamma}(a+N,b+T)[/math]. Questa proprietà consente di aggiornare le stime del tasso senza ricorrere necessariamente a integrazioni numeriche complesse.

Esercizio 1 — Tempi di consegna nel food delivery

Il problema

La startup QuickBite monitora i tempi di consegna in una zona urbana e, sulla base dei dati storici e di un modello statistico già stimato, assume che il tempo [math]X[/math], espresso in minuti, segua una distribuzione

[math]\displaystyle \begin{aligned}
X\sim\operatorname{Gamma}(3,0{,}25)
\end{aligned}[/math]

Il parametro di forma è quindi [math]\alpha=3[/math], mentre il parametro di tasso è [math]\beta=0{,}25[/math] minuti inversi.

Il responsabile delle operazioni vuole conoscere la densità in corrispondenza di dieci minuti, la probabilità di completare una consegna entro dieci minuti, la probabilità che il tempo sia compreso fra otto e quindici minuti e il tempo medio previsto dal modello.

Calcolo analitico della densità

Poiché [math]\alpha=3[/math], possiamo utilizzare [math]\Gamma(3)=2!=2[/math] nella formula della densità:

[math]\displaystyle \begin{aligned}
f(x)=\frac{0{,}25^3}{2}x^2e^{-0{,}25x}
\end{aligned}[/math]

Sostituendo [math]x=10[/math], otteniamo

[math]\displaystyle \begin{aligned}
f(10)
&=\frac{0{,}25^3}{2}\cdot10^2e^{-2{,}5}\\
&=0{,}0078125\cdot100\cdot e^{-2{,}5}\\
&\approx0{,}0641
\end{aligned}[/math]

Il risultato è una densità espressa in [math]\text{min}^{-1}[/math], non una probabilità. Per sapere quanto sia probabile osservare un tempo compreso, per esempio, fra nove e undici minuti, dovremmo calcolare [math]F(11)-F(9)[/math].

Calcolo delle probabilità

Per una Gamma con forma intera pari a tre, la funzione di ripartizione può essere espressa come

[math]\displaystyle \begin{aligned}
F(x)=1-e^{-\beta x}
\left(1+\beta x+\frac{(\beta x)^2}{2}\right)
\end{aligned}[/math]

La probabilità di completare una consegna in meno di dieci minuti è quindi

[math]\displaystyle \begin{aligned}
P(X<10)
&=F(10)\\
&=1-e^{-2{,}5}
\left(1+2{,}5+\frac{2{,}5^2}{2}\right)\\
&\approx0{,}4562
\end{aligned}[/math]

Secondo il modello, circa il 45,62% delle consegne viene completato prima dei dieci minuti.

Per calcolare la probabilità che il tempo sia compreso fra otto e quindici minuti, utilizziamo la differenza fra le due probabilità cumulative:

[math]\displaystyle \begin{aligned}
P(8<X<15)=F(15)-F(8)
\end{aligned}[/math]

Il calcolo numerico restituisce

[math]\displaystyle \begin{aligned}
P(8<X<15)\approx0{,}3997
\end{aligned}[/math]

Quasi il 40% delle consegne ricade quindi nell’intervallo considerato, sempre sotto l’ipotesi che la distribuzione Gamma scelta descriva adeguatamente i tempi osservati.

Calcolo del tempo medio

Il valore atteso è immediato da ottenere:

[math]\displaystyle \begin{aligned}
E[X]=\frac{\alpha}{\beta}
=\frac{3}{0{,}25}=12\text{ minuti}
\end{aligned}[/math]

La media teorica è dunque di dodici minuti, ma il risultato non significa che la maggior parte delle consegne duri esattamente dodici minuti.

La moda della distribuzione è invece

[math]\displaystyle \begin{aligned}
\operatorname{Moda}(X)
=\frac{3-1}{0{,}25}=8\text{ minuti}
\end{aligned}[/math]

La differenza fra otto e dodici minuti dipende dall’asimmetria della distribuzione: i tempi più lunghi contribuiscono ad aumentare la media, mentre il massimo della densità si trova in corrispondenza di un tempo inferiore.

Verifica con Python

Il codice seguente riproduce i calcoli analitici e genera due grafici, uno per la densità e uno per la funzione di ripartizione.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# Parametri della distribuzione Gamma
alpha = 3
beta = 0.25

# SciPy utilizza scale = 1 / beta
scale = 1 / beta

# Densità nel punto x = 10
densita_10 = gamma.pdf(10, a=alpha, scale=scale)

# Probabilità cumulative
prob_meno_10 = gamma.cdf(10, a=alpha, scale=scale)

# Probabilità che X sia compreso tra 8 e 15
prob_8_15 = (
    gamma.cdf(15, a=alpha, scale=scale)
    - gamma.cdf(8, a=alpha, scale=scale)
)

# Media, varianza e moda
media = alpha / beta
varianza = alpha / (beta ** 2)
moda = (alpha - 1) / beta

print(f"Densità in X = 10: {densita_10:.4f}")
print(f"P(X < 10): {prob_meno_10:.4f}")
print(f"P(8 < X < 15): {prob_8_15:.4f}")
print(f"Media: {media:.2f} minuti")
print(f"Varianza: {varianza:.2f} minuti²")
print(f"Moda: {moda:.2f} minuti")

# Grafico della densità
x = np.linspace(0, 40, 1000)
pdf = gamma.pdf(x, a=alpha, scale=scale)

plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="Densità Gamma", linewidth=2)

# Evidenzia l'intervallo tra 8 e 15 minuti
mask = (x >= 8) & (x <= 15)
plt.fill_between(
    x[mask],
    pdf[mask],
    alpha=0.35,
    label="Intervallo tra 8 e 15 minuti"
)

plt.axvline(
    media,
    linestyle="--",
    label=f"Media = {media:.0f} minuti"
)

plt.title("Distribuzione Gamma dei tempi di consegna")
plt.xlabel("Tempo di consegna (minuti)")
plt.ylabel("Densità di probabilità")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

I valori attesi in uscita sono:

f(10) = 0.0641 min^-1
P(X < 10) = 0.4562
P(8 < X < 15) = 0.3997
Media = 12.00 min
Varianza = 48.00 min^2
Moda = 8.00 min

Che cosa ci dicono i risultati?

Se il responsabile di QuickBite deve stabilire una promessa commerciale di consegna entro dieci minuti, il modello suggerisce che tale obiettivo sarebbe raggiunto soltanto da poco meno della metà degli ordini, assumendo che le condizioni operative restino compatibili con quelle dei dati storici.

Il valore medio di dodici minuti non è quindi sufficiente per definire una promessa di servizio affidabile. Per prendere una decisione più concreta, conviene calcolare il percentile 90 o 95, verificare se il modello riproduce adeguatamente la coda dei tempi osservati e confrontare i risultati fra zone, fasce orarie e giorni della settimana.

Riflessione

Se aumentiamo [math]\alpha[/math] mantenendo [math]\beta[/math] costante, che cosa succede alla distribuzione?

La media aumenta, perché [math]E[X]=\alpha/\beta[/math], e aumenta anche la varianza, perché [math]\operatorname{Var}(X)=\alpha/\beta^2[/math]. Tuttavia, il coefficiente di variazione

[math]\displaystyle \begin{aligned}
CV=\frac{\sigma_X}{E[X]}=\frac{1}{\sqrt{\alpha}}
\end{aligned}[/math]

diminuisce, indicando una dispersione relativa più contenuta. La distribuzione diventa progressivamente meno asimmetrica rispetto alla propria media, anche se la sua forma continua a dipendere dal parametro di forma.

Esercizio 2 — Latenza di un microservizio AI

Il problema

Un servizio di inferenza per modelli linguistici registra la latenza delle richieste API. Per costruire un esempio didattico, assumiamo che la latenza [math]X[/math], espressa in secondi, sia descritta da una distribuzione Gamma con parametri

[math]\displaystyle \begin{aligned}
X\sim\operatorname{Gamma}(4,8)
\end{aligned}[/math]

Il team di Reliability Engineering vuole valutare una soglia di 500 millisecondi, equivalente a [math]0{,}5[/math] secondi, e capire se il sistema rispetta un requisito temporale.

Le richieste consistono nel calcolare media, varianza, deviazione standard e moda, stimare la probabilità di superare la soglia e verificare i risultati con una simulazione Monte Carlo di 100.000 osservazioni.

Calcolo dei momenti

Con [math]\alpha=4[/math] e [math]\beta=8[/math], la media vale

[math]\displaystyle \begin{aligned}
E[X]=\frac{4}{8}=0{,}5\text{ s}
\end{aligned}[/math]

La varianza è

[math]\displaystyle \begin{aligned}
\operatorname{Var}(X)=\frac{4}{8^2}
=\frac{4}{64}=0{,}0625\text{ s}^2
\end{aligned}[/math]

La deviazione standard si ottiene estraendo la radice quadrata della varianza:

[math]\displaystyle \begin{aligned}
\sigma_X=\sqrt{0{,}0625}=0{,}25\text{ s}
\end{aligned}[/math]

La moda è

[math]\displaystyle \begin{aligned}
\operatorname{Moda}(X)=\frac{4-1}{8}
=0{,}375\text{ s}
\end{aligned}[/math]

La distribuzione ha quindi una media di 500 millisecondi, una deviazione standard di 250 millisecondi e una moda di 375 millisecondi.

Quanto spesso la latenza supera 500 millisecondi?

La probabilità di superare la soglia è

[math]\displaystyle \begin{aligned}
P(X>0{,}5)=1-F(0{,}5)
\end{aligned}[/math]

Poiché la forma è un intero positivo, possiamo utilizzare la formula della coda della distribuzione Gamma:

[math]\displaystyle \begin{aligned}
P(X>x)=e^{-\beta x}
\sum_{k=0}^{\alpha-1}\frac{(\beta x)^k}{k!}
\end{aligned}[/math]

Nel nostro caso, [math]\beta x=8\cdot0{,}5=4[/math], perciò

[math]\displaystyle \begin{aligned}
P(X>0{,}5)
&=e^{-4}\left(
1+4+\frac{4^2}{2!}+\frac{4^3}{3!}
\right)\\
&\approx0{,}4335
\end{aligned}[/math]

La probabilità teorica di superamento è quindi circa il 43,35%.

Questo risultato merita attenzione: la latenza media coincide con la soglia, ma il modello prevede che una quota consistente di richieste la superi. La media non è infatti una soglia che divide automaticamente le osservazioni in due gruppi uguali, soprattutto quando la distribuzione è asimmetrica.

Calcolo del percentile 95

Per descrivere le prestazioni della coda, calcoliamo anche il percentile 95, cioè il valore [math]x_{0,95}[/math] tale che

[math]\displaystyle \begin{aligned}
P(X\leq x_{0,95})=0{,}95
\end{aligned}[/math]

Il percentile si ricava invertendo numericamente la funzione di ripartizione. In Python possiamo utilizzare gamma.ppf(0.95, ...), ottenendo un valore di circa

[math]\displaystyle \begin{aligned}
x_{0,95}\approx0{,}9692\text{ s}
\end{aligned}[/math]

Secondo il modello, il 95% delle richieste ha una latenza non superiore a circa 969 millisecondi, mentre il restante 5% supera questo valore.

Un percentile non sostituisce automaticamente uno SLO, perché quest’ultimo deve precisare anche la finestra temporale di osservazione e la percentuale di richieste che devono rispettare il requisito. Il percentile è però un indicatore molto utile per descrivere il comportamento della coda e verificare quanto il sistema si discosti dalle aspettative.

Verifica con Python e simulazione Monte Carlo

La simulazione Monte Carlo genera 100.000 valori casuali dalla distribuzione ipotizzata, consentendo di confrontare la media campionaria e la quota di superamenti con le corrispondenti quantità teoriche.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# Parametri della distribuzione Gamma
alpha = 4
beta = 8
scale = 1 / beta

# Soglia di latenza in secondi
soglia = 0.5

# Probabilità teorica di superare la soglia
prob_teorica = gamma.sf(
    soglia,
    a=alpha,
    scale=scale
)

# 95° percentile della latenza
p95 = gamma.ppf(
    0.95,
    a=alpha,
    scale=scale
)

# Simulazione Monte Carlo
rng = np.random.default_rng(42)
n_simulazioni = 200_000

campione = rng.gamma(
    shape=alpha,
    scale=scale,
    size=n_simulazioni
)

prob_simulata = np.mean(campione > soglia)

print(f"P(X > {soglia}): {prob_teorica:.4f}")
print(f"Stima Monte Carlo: {prob_simulata:.4f}")
print(f"95° percentile (p95): {p95:.4f} secondi")
print(f"Latenza media teorica: {alpha / beta:.4f} secondi")

# Grafico della densità
x = np.linspace(0, 1.5, 1000)
pdf = gamma.pdf(x, a=alpha, scale=scale)

plt.figure(figsize=(10, 6))
plt.plot(x, pdf, linewidth=2, label="Densità Gamma")

# Evidenzia la coda oltre la soglia
mask = x >= soglia
plt.fill_between(
    x[mask],
    pdf[mask],
    alpha=0.35,
    label=f"Latenza oltre {soglia} secondi"
)

plt.axvline(
    soglia,
    linestyle="--",
    label=f"Soglia = {soglia} s"
)

plt.axvline(
    p95,
    linestyle=":",
    label=f"p95 = {p95:.3f} s"
)

plt.title("Distribuzione Gamma dei tempi di risposta")
plt.xlabel("Latenza (secondi)")
plt.ylabel("Densità di probabilità")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

I valori teorici sono:

Media teorica = 0.5000 s
Varianza teorica = 0.0625 s^2
Deviazione standard = 0.2500 s
Moda = 0.3750 s
P(X > 0.5) = 0.4335
Percentile 95 = 0.9692 s

La media simulata e la quota campionaria di superamenti saranno vicine ai valori teorici, ma non necessariamente identiche. La differenza dipende dalla variabilità del campione e tende generalmente a ridursi all’aumentare del numero di simulazioni, senza diventare necessariamente nulla.

Che cosa deve fare il team SRE?

Se il requisito è che almeno il 95% delle richieste venga completato entro 500 millisecondi, il modello attuale indica che il servizio è lontano dall’obiettivo: la quota teorica di richieste entro la soglia è soltanto il 56,65%.

Per individuare le cause, il team dovrebbe esaminare la latenza per tipo di richiesta, dimensione del prompt, lunghezza della risposta, carico del sistema e disponibilità delle risorse di calcolo. Prima di utilizzare la Gamma per previsioni operative, sarebbe inoltre necessario verificare la distribuzione sui dati reali, soprattutto nella coda destra, dove si concentrano i ritardi più problematici.

Riflessione

Perché la moda è inferiore alla media?

Per una Gamma con [math]\alpha>1[/math], la moda è [math](\alpha-1)/\beta[/math], mentre la media è [math]\alpha/\beta[/math]. La differenza è [math]1/\beta[/math], pari a [math]0{,}125[/math] secondi nel nostro esempio. La coda destra contiene valori elevati che contribuiscono ad aumentare la media, mentre il punto di massima densità rimane più vicino all’origine.

Esercizio 3 — Durata di una sequenza di fasi operative

Il problema

Un centro informatico deve completare un processo composto da quattro fasi consecutive prima di poter dichiarare conclusa un’operazione critica. Ogni fase può iniziare soltanto quando la precedente è terminata e la sua durata viene modellata come una variabile esponenziale indipendente con tasso

[math]\displaystyle \begin{aligned}
\lambda=\frac{1}{1500}\text{ h}^{-1}
\end{aligned}[/math]

Indichiamo con [math]X_i[/math] la durata della fase [math]i[/math], con [math]i=1,2,3,4[/math], e con [math]T[/math] il tempo totale necessario per completare l’intera sequenza:

[math]\displaystyle \begin{aligned}
T=X_1+X_2+X_3+X_4
\end{aligned}[/math]

Vogliamo calcolare la probabilità che il processo duri più di 3.000 ore, il tempo medio complessivo e il percentile 95.

È importante precisare che questo scenario descrive quattro fasi successive, non quattro server che lavorano contemporaneamente in ridondanza parallela. La distinzione cambia il modello probabilistico e quindi i risultati.

Perché la somma segue una distribuzione Gamma?

Le quattro durate sono indipendenti e seguono la stessa distribuzione esponenziale, con tasso [math]\lambda=1/1500[/math]. In queste condizioni, la loro somma segue una distribuzione Gamma con parametro di forma pari al numero delle fasi e parametro di tasso uguale a quello delle singole durate:

[math]\displaystyle \begin{aligned}
T\sim\operatorname{Gamma}
\left(4,\frac{1}{1500}\right)
\end{aligned}[/math]

Il parametro di scala corrispondente è

[math]\displaystyle \begin{aligned}
\theta=\frac{1}{\lambda}=1500\text{ ore}
\end{aligned}[/math]

Questa proprietà è valida per la somma di durate esponenziali indipendenti con lo stesso tasso; se i tassi fossero differenti, la somma non seguirebbe in generale una Gamma standard.

Probabilità di superare 3.000 ore

Per calcolare la probabilità che il processo non sia ancora terminato dopo 3.000 ore, utilizziamo la funzione di sopravvivenza:

[math]\displaystyle \begin{aligned}
P(T>3000)=1-F(3000)
\end{aligned}[/math]

Nel nostro caso,

[math]\displaystyle \begin{aligned}
\lambda t=\frac{3000}{1500}=2
\end{aligned}[/math]

Poiché il parametro di forma è quattro, possiamo applicare la formula della coda:

[math]\displaystyle \begin{aligned}
P(T>3000)
&=e^{-2}\left(
1+2+\frac{2^2}{2!}+\frac{2^3}{3!}
\right)\\
&=e^{-2}\left(1+2+2+\frac{4}{3}\right)\\
&\approx0{,}8571
\end{aligned}[/math]

La probabilità teorica di completare il processo dopo più di 3.000 ore è quindi circa l’85,71%.

Il risultato non è la probabilità che un singolo componente rimanga operativo per 3.000 ore, ma la probabilità che la somma delle quattro durate superi tale soglia.

Tempo medio e percentile 95

La media della somma è

[math]\displaystyle \begin{aligned}
E[T]=\frac{\alpha}{\beta}
=\frac{4}{1/1500}=6000\text{ ore}
\end{aligned}[/math]

Il percentile 95 è il tempo entro cui, secondo il modello, il 95% dei processi dovrebbe risultare completato. Per determinarlo, dobbiamo risolvere numericamente

[math]\displaystyle \begin{aligned}
F(t_{0,95})=0{,}95
\end{aligned}[/math]

Il calcolo con SciPy restituisce un valore di circa

[math]\displaystyle \begin{aligned}
t_{0,95}\approx11630\text{ ore}
\end{aligned}[/math]

Questo significa che il tempo medio è pari a 6.000 ore, mentre il 95° percentile si colloca intorno a 11.630 ore. La differenza dipende dalla dispersione e dall’asimmetria della distribuzione: una quota ridotta di processi può richiedere molto più tempo della media.

Verifica con Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# Parametri delle quattro fasi esponenziali
n_fasi = 4
lambda_rate = 1 / 1500  # Tasso per ora
scale = 1 / lambda_rate

# Orizzonte temporale
tempo = 3000

# Distribuzione del tempo totale: Gamma
prob_teorica = gamma.sf(
    tempo,
    a=n_fasi,
    scale=scale
)

# Media e varianza del tempo totale
media = n_fasi / lambda_rate
varianza = n_fasi / (lambda_rate ** 2)

# Simulazione Monte Carlo della somma delle quattro fasi
rng = np.random.default_rng(42)
n_simulazioni = 200_000

durate = rng.exponential(
    scale=scale,
    size=(n_simulazioni, n_fasi)
)

tempi_totali = durate.sum(axis=1)
prob_simulata = np.mean(tempi_totali > tempo)

print(f"P(T > {tempo} ore): {prob_teorica:.4f}")
print(f"Stima Monte Carlo: {prob_simulata:.4f}")
print(f"Tempo medio totale: {media:.2f} ore")
print(f"Varianza del tempo totale: {varianza:.2f} ore²")

# Grafico della densità
x = np.linspace(0, 10000, 1000)
pdf = gamma.pdf(
    x,
    a=n_fasi,
    scale=scale
)

plt.figure(figsize=(10, 6))
plt.plot(x, pdf, linewidth=2, label="Densità Gamma")

# Evidenzia i tempi superiori a 3000 ore
mask = x >= tempo
plt.fill_between(
    x[mask],
    pdf[mask],
    alpha=0.35,
    label=f"Tempo totale oltre {tempo} ore"
)

plt.axvline(
    tempo,
    linestyle="--",
    label=f"Soglia = {tempo} ore"
)

plt.axvline(
    media,
    linestyle=":",
    label=f"Media = {media:.0f} ore"
)

plt.title("Distribuzione Gamma del tempo totale delle quattro fasi")
plt.xlabel("Tempo totale (ore)")
plt.ylabel("Densità di probabilità")
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

L’output atteso è:

P(T > 3000 h) = 0.8571
Media = 6000 h
Varianza = 9000000 h^2
Percentile 95 = 11630 h

E se i quattro server funzionassero in parallelo?

Consideriamo ora un caso diverso, più vicino a un sistema con ridondanza parallela: quattro server indipendenti iniziano a funzionare contemporaneamente e il cluster rimane operativo finché almeno uno dei quattro server è ancora attivo.

In questo caso il tempo fino al guasto totale non è la somma delle durate, ma il massimo dei quattro tempi di vita:

[math]\displaystyle \begin{aligned}
T_{\max}=\max(X_1,X_2,X_3,X_4)
\end{aligned}[/math]

Se ciascun server ha una durata esponenziale con tasso [math]\lambda[/math], la probabilità che tutti e quattro siano già guasti al tempo [math]t[/math] è

[math]\displaystyle \begin{aligned}
P(T_{\max}\leq t)
=\left(1-e^{-\lambda t}\right)^4
\end{aligned}[/math]

Di conseguenza, la probabilità che il cluster rimanga operativo oltre [math]t[/math] è

[math]\displaystyle \begin{aligned}
P(T_{\max}>t)
=1-\left(1-e^{-\lambda t}\right)^4
\end{aligned}[/math]

Per [math]t=3000[/math] ore, con [math]\lambda=1/1500[/math], otteniamo

[math]\displaystyle \begin{aligned}
P(T_{\max}>3000)
=1-(1-e^{-2})^4
\approx0{,}4411
\end{aligned}[/math]

Il valore è molto diverso dall’85,71% calcolato per la somma di quattro fasi. Non si tratta di una contraddizione: stiamo rispondendo a due domande differenti, perché in un caso sommiamo le durate di attività successive e nell’altro consideriamo il tempo fino al guasto dell’ultimo server rimasto attivo.

Anche la media cambia. Per quattro server indipendenti con lo stesso tasso esponenziale, il tempo medio fino al guasto di tutti è

[math]\displaystyle \begin{aligned}
E[T_{\max}]
=\frac{1+\frac12+\frac13+\frac14}{\lambda}
=3125\text{ ore}
\end{aligned}[/math]

Questa distinzione è essenziale quando si modellano architetture ad alta disponibilità: prima di scegliere una distribuzione, bisogna definire con precisione quale evento determina la fine del servizio.

Riflessione

Che cosa succede se le quattro fasi consecutive hanno tassi di completamento differenti?

Se le durate rimangono indipendenti ma seguono distribuzioni esponenziali con tassi diversi, la somma non segue più, in generale, una Gamma standard. Nel caso di tassi distinti, la distribuzione risultante è detta ipoesponenziale e può essere descritta attraverso la funzione generatrice dei momenti:

[math]\displaystyle \begin{aligned}
M_T(s)=\prod_{i=1}^{n}\frac{\lambda_i}{\lambda_i-s}
\end{aligned}[/math]

La funzione è definita per i valori di [math]s[/math] per cui il prodotto converge. Per calcolare probabilità e quantili in questo scenario possiamo utilizzare formule specifiche oppure metodi numerici, mantenendo esplicite le ipotesi di indipendenza e di tasso costante.

Come scegliere e verificare un modello Gamma

I tre esercizi mostrano che la distribuzione Gamma può essere utile per descrivere tempi positivi e asimmetrici, ma la scelta del modello non dovrebbe dipendere soltanto dalla comodità delle formule.

Prima di utilizzarla su dati reali, è opportuno controllare che le osservazioni siano coerenti con le ipotesi del modello, confrontare la distribuzione empirica con quella teorica e verificare soprattutto il comportamento della coda destra, che spesso contiene gli eventi più rilevanti dal punto di vista operativo.

È altrettanto importante distinguere fra parametri stimati dai dati e valori assegnati a scopo didattico. Negli esercizi di questo articolo i parametri sono fissati per rendere trasparenti i calcoli; in un progetto reale dovrebbero invece essere stimati e sottoposti a verifiche di adeguatezza, tenendo conto anche di eventuali dipendenze, cambiamenti nel tempo e osservazioni anomale.

La distribuzione Gamma non è quindi una soluzione automatica per qualsiasi fenomeno caratterizzato da valori positivi. È uno strumento utile quando la struttura del problema e i dati disponibili ne sostengono l’impiego, mentre la lettura dei risultati deve sempre tenere conto delle ipotesi con cui sono stati ottenuti.

In sintesi

La distribuzione Gamma consente di descrivere tempi positivi e asimmetrici attraverso pochi parametri, collegando formule analitiche, calcolo delle probabilità e simulazione in Python.

Nel caso delle consegne, abbiamo distinto la densità dalla probabilità e osservato come la media possa differire dalla moda. Nell’esercizio sulla latenza AI, abbiamo visto che conoscere il tempo medio non basta per valutare il rispetto di una soglia e che i quantili forniscono informazioni più dirette sul comportamento della coda. Infine, il problema delle fasi operative ha mostrato quanto sia importante definire correttamente l’evento studiato: la somma di durate successive e il tempo fino al guasto di un sistema ridondante sono grandezze diverse, anche quando coinvolgono componenti con le stesse distribuzioni individuali.

La lezione più importante riguarda il rapporto fra matematica e decisioni: un modello probabilistico è utile quando consente di quantificare l’incertezza senza nascondere le ipotesi da cui dipendono i risultati.

---

Dalla teoria alla stima statistica: continua con la seconda parte

Conoscere la distribuzione Gamma e saper calcolare probabilità non basta quando dobbiamo analizzare dati reali. Occorre stimare i parametri a partire dalle osservazioni, valutare l’incertezza delle stime e capire se un unico modello descrive adeguatamente fenomeni che possono seguire regimi differenti.

Nella seconda parte passiamo dalla teoria all’inferenza statistica con tre esercizi svolti in Python, dedicati alla stima di massima verosimiglianza (MLE), all’approccio bayesiano Gamma-Poisson e alle mixture di distribuzioni Gamma, con esempi applicati a latenze dei servizi, sinistri assicurativi e precipitazioni.

👉 Continua con la seconda parte: Distribuzione Gamma in Python, MLE, Bayes e Mixture con 3 esercizi svolti

Pubblicità