Distribuzione Gamma: 6 Esercizi con Soluzioni e Codice per Calcolare Probabilità, Media e Inferenza Bayesiana

Quante volte, aspettando in fila o usando un nuovo dispositivo, ti sei chiesto: “Quanto ancora ci vorrà?” oppure “Quanto durerà questo prodotto prima di rompersi?”. Queste non sono solo domande quotidiane, ma problemi centrali in molti campi, dall’ingegneria alla finanza. La statistica ci offre uno strumento incredibilmente versatile per dare una risposta: la distribuzione Gamma.

Lungi dall’essere un concetto astratto, la distribuzione Gamma modella il tempo d’attesa non per un singolo evento, ma per una serie di eventi. Immaginala come la “cugina” più esperta della distribuzione esponenziale. Se l’esponenziale calcola l’attesa per il prossimo cliente, la Gamma calcola l’attesa per il quinto, o il ventesimo cliente.

In questa guida, esploreremo la sua potenza attraverso 6 esercizi pratici di difficoltà crescente. Partiremo da un semplice smoothie per arrivare a complessi problemi di controllo qualità e inferenza bayesiana, mostrando come questo modello passi dai banchi dell’università alle decisioni strategiche di un’azienda.

Esercizio 1 – Facile

Il tempo di attesa del tuo smoothie

Al bar dell’istituto il tempo di attesa X (in minuti) per un smoothie segue una distribuzione Gamma di parametri shape [math]k=2[/math] e rate [math]\theta=0.5\ \text{min}^{-1}[/math].

1. Qual è la probabilità di aspettare meno di 3 minuti?
2. Calcola il tempo medio di attesa e la deviazione standard.

Risoluzione

1) Richiamo teorico
Una v.a. [math]X \sim \text{Gamma}(k,\theta)[/math] ha densità

[math]\displaystyle f(x)=\frac{\theta^k x^{k-1}e^{-\theta x}}{\Gamma(k)},\qquad x\ge 0[/math]

dove [math]\Gamma(k)[/math] è la funzione gamma di Eulero.
Per [math]k[/math] intero si ha [math]\Gamma(k)=(k-1)![/math].

2) Probabilità [math]P(X<3)[/math]
Poiché [math]k=2[/math] è intero, la Gamma è la somma di 2 esponenziali indipendenti di parametro [math]\theta=0.5[/math].
La funzione di ripartizione ha formula chiusa:

[math]\displaystyle F(x)=1-e^{-\theta x}\sum_{n=0}^{k-1}\frac{(\theta x)^n}{n!}[/math]

Sostituendo:

[math]\displaystyle \begin{aligned}
P(X<3)&=F(3)=1-e^{-0.5\cdot 3}\left(1+\frac{0.5\cdot 3}{1!}\right)\\
&=1-e^{-1.5}(1+1.5)=1-2.5e^{-1.5}\\
&\approx 1-2.5\cdot 0.22313=0.442\ (44.2\%).
\end{aligned}[/math]

3) Media e deviazione standard
Formule generali:

[math]\displaystyle \mathbb E[X]=\frac{k}{\theta},\qquad \text{SD}[X]=\frac{\sqrt{k}}{\theta}[/math]
[math]\displaystyle \Rightarrow\ \mathbb E[X]=\frac{2}{0.5}=4\ \text{min},\qquad \text{SD}[X]=\frac{\sqrt{2}}{0.5}=2.83\ \text{min}.[/math]

Mini-quiz
✍️ Cosa cambierebbe nella formula di [math]F(x)[/math] se [math]k=1.5[/math] (non intero)?
Spunto: la somma finita diventa…

Esercizio 2 – Medio-facile

La garanzia dello smartphone

La durata in anni T della batteria di uno smartphone è modellata da una Gamma di media=3 anni e varianza=1.5 anni. Il produttore vuole offrire una garanzia g tale che solo il 5 % dei dispositivi venga sostituito. Trova g.

Risoluzione

1) Ricavo parametri k,θ
Sistema

[math]\displaystyle \mathbb E[T]=\frac k\theta=3,\qquad \text{Var}(T)=\frac k{\theta^2}=1.5[/math]

Divido membro a membro:

[math]\displaystyle \frac{k/\theta}{k/\theta^2}=\theta=\frac{3}{1.5}=2\ \text{anni}^{-1}\ \Rightarrow\ k=6.[/math]

2) Percentile 5 °
Devo risolvere

[math]\displaystyle P(T\le g)=0.05,\qquad T\sim\text{Gamma}(k=6,\theta=2).[/math]

# Python (SciPy)
from scipy.stats import gamma
g = gamma.ppf(0.05, a=6, scale=1/2)   # scale = 1/θ
print(g)        # ≈ 1.3033230482592888 anni

Pertanto, la garanzia g che il produttore dovrebbe offrire affinché solo il 5% dei dispositivi venga sostituito è di circa 1.30 anni (corrispondenti a circa 15 mesi e mezzo).

Mini-quiz
✍️ Come cambierebbe g se la varianza fosse minore (maggiore precisione di produzione)?
Spunto: disegna le due densità…

Esercizio 3 – Medio

Il contatore Geiger

Una sorgente radioattiva emette in media 60 particelle al secondo. Il tempo che intercorre tra l’inizio dell’osservazione e il rilevamento della 20° particella è modellato da una Gamma.

• Scrivi i parametri della distribuzione.
• Calcola la probabilità che il contatore impieghi più di [math]0.4\ \text{s}[/math].
• Determina l’intervallo di tempo [math][t_1,t_2][/math] entro cui cade l’evento con probabilità [math]0.90[/math] (simmetrico rispetto alla mediana).

Risoluzione

1) Parametri

I tempi tra arrivi successivi sono [math]\text{Exp}(\lambda=60\ \text{s}^{-1})[/math]. La somma di [math]n=20[/math] esponenziali indipendenti è [math]X\sim\text{Gamma}(k=20,\theta=60).[/math]

2) [math]P(X>0.4)[/math]

Uso la coda superiore: [math]\displaystyle P(X>0.4)=1-F(0.4;20,60).[/math]

# R
1-pgamma(0.4, shape=20, rate=60)  # ≈ 0.021

3) Intervallo [math]90%[/math] simmetrico

Voglio [math]t_1, t_2[/math] tali che [math]\displaystyle P(t_1\le X\le t_2)=0.90,\qquad P(X<t_1)=P(X>t_2)=0.05.[/math]

t1 <- qgamma(0.05, 20, 60)
t2 <- qgamma(0.95, 20, 60)

t1 ≈ 0.231 s, t2 ≈ 0.424 s

💡 Osservazioni

Per [math]k[/math] “grande” la Gamma si avvicina a una Normale (CLT): confronta i percentili con [math]\mu\pm1.645\sigma[/math]. La media teorica è [math]k/\theta=20/60\approx 0.333\ \text{s}[/math], perfettamente centrata nell’intervallo trovato.

Mini-quiz

✍️ Cosa succede alla larghezza dell’intervallo [math]90%[/math] se raddoppio il tempo di osservazione (quindi [math]k=40[/math])?

Spunto: calcola il coefficiente di variazione…

Esercizio 4 – Medio-difficile

Il test di resistenza dei tessuti

Un laboratorio tessile sottopone un campione a cicli di trazione fino a rottura. Il numero di cicli [math]N[/math] segue una Gamma di media=[math]400[/math] cicli e varianza=[math]8000[/math] cicli.

Per superare la certificazione, il campione deve resistere almeno [math]350[/math] cicli.

• Calcola la probabilità di superare il test.
• Se testo 10 campioni indipendenti, qual è la probabilità che almeno 8 superino la certificazione?

Risoluzione

1) Parametri

[math]\displaystyle \frac k\theta=400,\quad \frac k{\theta^2}=8000\ \Rightarrow\ \theta=0.05,\ k=20.[/math]

Probabilità di superare 350 cicli: [math]\displaystyle P(N\ge 350)=1-F(350;20,0.05)=1-\text{pgamma}(350,20,0.05)\approx 0.762.[/math]

2) Modello binomiale

Sia [math]Y=\text{# successi in 10 prove}[/math], [math]Y\sim\text{Bin}(n=10,p=0.762).[/math]

[math]\displaystyle P(Y\ge 8)=\sum_{y=8}^{10}\binom{10}{y}p^y(1-p)^{10-y}=0.47\ (\text{calcolo diretto}).[/math]

💡 Osservazioni

Abbiamo “innestato” una binomiale sopra una Gamma: esempio classico di modello gerarchico. Se [math]p[/math] fosse stato [math]<0.5[/math], la probabilità di “almeno 8 su 10” sarebbe crollata: prova a cambiare soglia a 450 cicli.

Mini-quiz

✍️ Come cambierebbe la risposta al punto 2 se i campioni fossero correlati positivamente (stessa macchina, stesso lotto)?

Spunto: la varianza di [math]Y[/math] aumenterebbe o diminuirebbe?

Esercizio 5 – Difficile

Bayes e la durata delle lampadine

Il tempo di vita (ore) di un nuovo tipo di lampadina è X~Gamma(k,θ).
Prima di raccogliere dati, l’ingegnore esprime l’incertezza su θ con un prior θ~Gamma(α=4, β=1000) (parametrizzazione shape–rate).
Osserva poi n=8 lampadine che durano in media 1200 ore e decide di usare k=3 fisso (tecnologia nota).

1. Determina il posterior di θ.
2. Calcola il valor medio a posteriori di θ e confrontalo con la stima di massima verosimiglianza (MLE).

Risoluzione

1) Modello coniugato
Gamma è coniugata al rate di un’altra Gamma.
Likelihood:

[math]\displaystyle L(\theta)\propto \theta^{nk}\exp(-\theta\sum x_i)=\theta^{24}\exp(-\theta\cdot 9600).[/math]

Prior:

[math]\displaystyle \pi(\theta)\propto \theta^{\alpha-1}e^{-\beta\theta}=\theta^{3}e^{-1000\theta}.[/math]

Posterior:

[math]\displaystyle \pi(\theta|\text{data})\propto \theta^{24+3}e^{-(9600+1000)\theta}= \theta^{27}e^{-10600\theta}.[/math]

La densità proporzionale a θ²⁷ e^−10600θ è

[math]\displaystyle \theta|\text{data}\sim\text{Gamma}(\alpha’=28,\beta’=10600).[/math]

2) Media a posteriori

[math]\displaystyle \mathbb E[\theta|\text{data}]=\frac{28}{10600}\approx 0.00264\ \text{ore}^{-1}.[/math]

MLE:

[math]\displaystyle \hat\theta_{\text{MLE}}=\frac{nk}{\sum x_i}=\frac{24}{9600}=0.00250\ \text{ore}^{-1}.[/math]

Il prior ha “spostato” leggermente la stima verso l’alto (shrinkage bayesiano).

Perché la coniugazione Gamma-Gamma  semplifica così tanto i calcoli nell’inferenza bayesiana

Il Problema di Partenza: Il Teorema di Bayes e la sua Complessità

Per capire il “trucco”, dobbiamo prima capire il problema che risolve. Il cuore dell’inferenza bayesiana è il Teorema di Bayes, che ci dice come aggiornare le nostre credenze su un parametro [math]\theta[/math] (ad esempio, il tasso di guasto di una lampadina) dopo aver osservato dei dati.

La sua forma è:

[math]
P(\theta|dati) = \frac{P(dati|\theta) \cdot P(\theta)}{P(dati)}
[/math]

Dove:

• [math]P(\theta|dati)[/math] è la distribuzione a posteriori (la nostra credenza aggiornata su [math]\theta[/math]).
• [math]P(dati|\theta)[/math] è la verosimiglianza (likelihood), che ci dice quanto sono probabili i dati che abbiamo osservato, dato un certo valore di [math]\theta[/math].
• [math]P(\theta)[/math] è la distribuzione a priori (la nostra credenza iniziale su [math]\theta[/math], prima di vedere i dati).
• [math]P(dati)[/math] è l’evidenza marginale, ovvero la probabilità dei dati.

Il vero problema computazionale è il denominatore, [math]P(dati)[/math]. Per calcolarlo, dovremmo risolvere un integrale spesso molto complesso su tutti i possibili valori di [math]\theta[/math]:

[math]
P(dati) = \int P(dati|\theta) \cdot P(\theta)d\theta
[/math]

Risolvere questo integrale è spesso difficile o addirittura impossibile analiticamente. Questo è il collo di bottiglia che ha reso l’approccio bayesiano poco pratico per decenni.

La Soluzione Elegante: La Coniugazione

Qui entra in gioco la coniugazione.

Una distribuzione a priori si dice coniugata a una certa verosimiglianza se la distribuzione a posteriori risultante appartiene alla stessa famiglia della distribuzione a priori.

È come una reazione chimica perfetta: mescoli un ingrediente di tipo A (il prior) con un reagente di tipo B (la verosimiglianza) e ottieni un prodotto che è ancora di tipo A, solo con proprietà leggermente modificate.

Il Caso Specifico: Gamma-Gamma

Vediamo come funziona nel nostro esempio (Esercizio 5).

La Verosimiglianza (Likelihood)

I nostri dati (durata delle lampadine) provengono da una distribuzione [math]\text{Gamma}(k, \theta)[/math], dove [math]k[/math] è noto e [math]\theta[/math] è il parametro incognito. La funzione di verosimiglianza per [math]\theta[/math], dopo aver osservato [math]n[/math] lampadine, ha una forma matematica proporzionale a:

[math]
L(\theta) \propto \theta^{nk} \exp(-\theta \sum x_i)
[/math]

Nota la struttura: è una potenza di [math]\theta[/math] moltiplicata per un’esponenziale di [math]-\theta[/math].

La Scelta del Prior “Furbo”

Dobbiamo scegliere un prior per [math]\theta[/math]. Osservando la forma della verosimiglianza, ci chiediamo: “Esiste una famiglia di distribuzioni di probabilità che ha proprio quella forma, [math]\theta[/math] elevato a qualcosa per [math]e[/math] elevato a [math]-\theta[/math] per qualcos’altro?”. La risposta è sì: la distribuzione Gamma!

Scegliamo un prior [math]\text{Gamma}(\alpha, \beta)[/math]. La sua densità di probabilità è proporzionale a:

[math]
P(\theta) \propto \theta^{\alpha-1} e^{-\beta\theta}
[/math]

Il “Miracolo” Matematico:

Ora calcoliamo il posteriore. Poiché il denominatore [math]P(dati)[/math] è solo una costante di normalizzazione, possiamo ignorarlo e lavorare con le proporzionalità:

[math]
\displaystyle
\begin{aligned}
P(\theta|dati) &\propto L(\theta) \cdot P(\theta) \\
&\propto [\theta^{nk} \exp(-\theta \sum x_i)] \cdot [\theta^{\alpha-1} e^{-\beta\theta}]
\end{aligned}
[/math]

Usando le proprietà delle potenze e degli esponenziali, raggruppiamo i termini:

[math]
P(\theta|dati) \propto \theta^{(nk+\alpha)-1} \exp(-(\sum x_i + \beta)\theta)
[/math]

Osserva attentamente questa formula finale. Ha esattamente la stessa struttura matematica di una distribuzione Gamma! Abbiamo trovato la nostra distribuzione a posteriori senza risolvere alcun integrale.

Il posteriore è una distribuzione [math]\text{Gamma}(\alpha’, \beta’)[/math], dove i nuovi parametri sono ottenuti con una semplice somma:

Perché è un “Trucco Elegante”?

• Evita l’Integrazione Complessa: Abbiamo completamente aggirato il calcolo dell’integrale al denominatore. L’eleganza sta nel fatto che, scegliendo il prior giusto, sappiamo già a priori quale sarà la “forma” della soluzione.
• Trasforma un Problema di Calcolo in un Semplice Aggiornamento:  L’inferenza bayesiana diventa una sorta di ricetta. Parti con dei parametri “iniziali” ([math]\alpha[/math], [math]\beta[/math]). I dati ti forniscono due numeri riassuntivi ([math]nk[/math] e [math]\sum x_i[/math]). Per ottenere i parametri “aggiornati” ([math]\alpha’[/math], [math]\beta’[/math]), devi solo fare due addizioni. È un processo computazionalmente banale.
• È Intuitivo: I parametri hanno un’interpretazione chiara. [math]\alpha[/math] può essere visto come il numero di “osservazioni fittizie” della nostra conoscenza pregressa, e [math]\beta[/math] come la loro “somma fittizia”. L’aggiornamento bayesiano non fa altro che aggiungere le nuove osservazioni ([math]nk[/math]) e le nuove somme ([math]\sum x_i[/math]) a quelle che già “conoscevamo”.

In sintesi, la coniugazione Gamma-Gamma è un “trucco” perché sfrutta la somiglianza strutturale tra la verosimiglianza e una famiglia di distribuzioni (la Gamma) per trasformare un problema di integrazione analiticamente intrattabile in un semplice esercizio di algebra. Questo non solo fa risparmiare tempo e fatica, ma rende il processo di aggiornamento delle credenze trasparente e interpretabile.

Mini-quiz
✍️ Cosa succederebbe al posterior se k fosse incerto e volessimo stimare anche k?
Spunto: la coniugazione decade…

Esercizio 6 – Difficile (+)

Controllo qualità in tempo reale

Una macchina produce viti in serie. Il tempo di rottura T (ore) è Gamma(k=5,θ=0.2).
Il capo-reparto vuole un allarme che scatti quando la probabilità condizionata che la vite duri ancora meno di 10 ore, dato che ha già lavorato s ore, supera il 10 %.
Determina la soglia di durata s* oltre la quale l’allarme deve entrare in azione.

Risoluzione

1) Proprietà di “memoryless” generalizzata
Per una Gamma la funzione di hazard è

[math]\displaystyle h(s)=\frac{f(s)}{1-F(s)}=\frac{\theta^k s^{k-1}e^{-\theta s}}{\Gamma(k)\ P(T>s)}.[/math]

Ci serve la probabilità condizionata

[math]\displaystyle P(T<s+10\mid T>s)=\frac{F(s+10)-F(s)}{1-F(s)}\ge 0.10.[/math]

2) Soluzione numerica 

import numpy as np
from scipy.stats import gamma
from scipy.optimize import brentq

def probabilita_condizionata(s, k=5, scala=5):
    """
    Calcola la probabilità condizionata P(T < s + 10 | T > s) per una distribuzione Gamma.
    
    Parametri:
    - s: Tempo di utilizzo in ore.
    - k: Parametro di forma della distribuzione Gamma.
    - scala: Parametro di scala della distribuzione Gamma (1/theta).
    
    Ritorna:
    - Probabilità condizionata.
    """
    cdf_s = gamma.cdf(s, a=k, scale=scala)
    cdf_s_plus_10 = gamma.cdf(s + 10, a=k, scale=scala)
    return (cdf_s_plus_10 - cdf_s) / (1 - cdf_s)

def funzione_obiettivo(s, k=5, scala=5, soglia=0.10):
    """
    Funzione obiettivo per trovare s* dove P(T < s + 10 | T > s) = soglia.
    
    Parametri:
    - s: Tempo in ore.
    - k: Parametro di forma della distribuzione Gamma.
    - scala: Parametro di scala della distribuzione Gamma.
    - soglia: Soglia della probabilità condizionata (default 0.10).
    
    Ritorna:
    - Valore della funzione obiettivo (zero quando s = s*).
    """
    return probabilita_condizionata(s, k, scala) - soglia

# Parametri della distribuzione
k = 5  # parametro di forma
theta = 0.2  # parametro di tasso
scala = 1 / theta  # parametro di scala = 1/theta

# Trova s* usando brentq
try:
    s_star = brentq(funzione_obiettivo, 0, 50, args=(k, scala))
    print(f"La soglia s* è: {s_star:.4f} ore")
except ValueError as e:
    print(f"Errore nella ricerca della radice: {e}")

# Verifica il risultato
prob_s_star = probabilita_condizionata(s_star, k, scala)
print(f"Probabilità condizionata a s* = {s_star:.4f}: {prob_s_star:.4f}")

Output:

La soglia s* è: 2.1661 ore

Probabilità condizionata a s* = 2.1661: 0.1000

Pertanto l’allarme scatta dopo circa 2.17 ore di lavoro.