42 La verosimiglianza

In questo capitolo, imparerai a:

Comprendere il concetto di verosimiglianza: Scoprirai il ruolo fondamentale che la verosimiglianza svolge nella stima dei parametri statistici.
Generare grafici della funzione di verosimiglianza:
- Implementare grafici per la funzione di verosimiglianza nel caso binomiale.
Interpretare i grafici della funzione di verosimiglianza: Sviluppare le competenze necessarie per analizzare e trarre conclusioni dai grafici generati.
Capire il principio di stima di massima verosimiglianza (MLE): Approfondiremo il metodo di stima di massima verosimiglianza.

Prerequisiti

Leggere il capitolo Estimation (Schervish & DeGroot, 2014).
Leggere il capitolo Bayes’ rule (Johnson et al., 2022).

Preparazione del Notebook

here::here("code", "_common.R") |> 
  source()

Introduzione

I ricercatori utilizzano diversi modelli matematici per descrivere e prevedere il comportamento dei dati osservati. Questi modelli si distinguono tra loro per la struttura funzionale, ovvero il modo in cui collegano le variabili osservate con parametri teorici. La scelta del modello migliore avviene confrontando le previsioni teoriche generate dal modello con i dati effettivamente osservati. Il modello che produce previsioni più vicine ai dati reali viene considerato il più adeguato per descrivere il fenomeno studiato.

In questo processo di confronto, la funzione di verosimiglianza gioca un ruolo fondamentale. Essa quantifica la probabilità che i dati osservati siano stati generati da un particolare modello con determinati valori dei suoi parametri. In altre parole, la verosimiglianza misura quanto i dati siano compatibili con il modello ipotizzato.

42.1 Il Principio della Verosimiglianza

La verosimiglianza misura quanto sia plausibile ciascun valore dei parametri alla luce dei dati osservati. In altre parole, indica quanto ogni possibile valore dei parametri sia compatibile con i dati raccolti.

Definizione 42.1 Consideriamo un vettore aleatorio \(Y\), la cui distribuzione è descritta da una funzione di densità (nel caso continuo) oppure da una funzione di massa di probabilità (nel caso discreto), indicata con \(f(y \mid \theta)\), dove \(\theta\) è un vettore di parametri appartenente allo spazio parametrico \(\Theta\).

Una volta osservato un valore specifico \(y\) del vettore \(Y\), definiamo la funzione di verosimiglianza come una funzione che associa a ciascun valore possibile dei parametri \(\theta\) la plausibilità di aver osservato proprio quei dati:

\[ L(\theta; y) = f(y \mid \theta). \]

In questa espressione, i dati osservati \(y\) sono considerati fissi, mentre la variabile di interesse è \(\theta\). La funzione di verosimiglianza esprime quindi quanto ogni possibile valore di \(\theta\) sia compatibile con i dati osservati.

42.1.1 Relazione tra Verosimiglianza e Funzione di Probabilità

Sia la funzione di probabilità (o densità) sia la funzione di verosimiglianza sono costruite sulla stessa espressione matematica: \(f(y \mid \theta)\). Tuttavia, il significato che attribuiamo a questa espressione cambia radicalmente a seconda del contesto inferenziale in cui ci troviamo.

La differenza tra funzione di probabilità e funzione di verosimiglianza riguarda il ruolo epistemologico assegnato ai dati e ai parametri:

Funzione di densità (o massa) di probabilità:
In questo caso, assumiamo che i parametri \(\theta\) siano noti e consideriamo i dati \(y\) come variabili aleatorie. La funzione \(f(y \mid \theta)\) rappresenta quindi il meccanismo generativo dei dati: ci dice quanto è probabile (o densa) l’osservazione di \(y\), se \(\theta\) è fissato.
Funzione di verosimiglianza:
Qui la prospettiva si inverte: i dati \(y\) sono fissi perché già osservati, mentre i parametri \(\theta\) sono incogniti e rappresentano l’oggetto dell’inferenza. La funzione \(L(\theta; y) = f(y \mid \theta)\) misura quanto ciascun valore possibile di \(\theta\) sia compatibile con i dati osservati.

Formalmente, la relazione è:

\[ L(\theta; y) = f(y \mid \theta) \]

ma l’interpretazione è diversa:

\(f(y \mid \theta)\) → probabilità di osservare \(y\), dato \(\theta\) (parametri fissi, dati variabili);
\(L(\theta; y)\) → plausibilità di \(\theta\), dati gli \(y\) osservati (dati fissi, parametri variabili).

In sintesi:

La funzione di probabilità risponde alla domanda:
“Se i parametri fossero questi, quanto è probabile osservare questi dati?”
La funzione di verosimiglianza risponde alla domanda:
“Dati questi dati, quali valori dei parametri sono più plausibili?”

Questa distinzione è fondamentale per l’inferenza statistica: mentre la funzione di probabilità descrive il processo generativo dei dati, la verosimiglianza gioca un ruolo centrale nell’aggiornamento delle credenze sui parametri.

In particolare, nella prospettiva bayesiana, la verosimiglianza fornisce l’informazione derivante dai dati osservati, che viene combinata con le credenze precedenti (prior) per ottenere la distribuzione a posteriori dei parametri.

In sintesi, la verosimiglianza ci dice quanto i dati supportano i diversi valori dei parametri. Nella visione bayesiana, essa è lo strumento attraverso cui i dati modificano le nostre credenze.

42.2 La Log-Verosimiglianza

Per ragioni sia analitiche che computazionali, si lavora spesso con la log-verosimiglianza, ovvero il logaritmo della funzione di verosimiglianza:

\[ \ell(\theta; y) = \log L(\theta; y) = \log f(y \mid \theta). \]

La log-verosimiglianza presenta numerosi vantaggi:

Numerici: i prodotti di probabilità molto piccole possono causare problemi di underflow numerico. Il logaritmo trasforma questi prodotti in somme, che sono numericamente più stabili.
Analitici: derivare la log-verosimiglianza rispetto ai parametri è spesso più semplice, facilitando la stima mediante metodi come il massimo della verosimiglianza (MLE).
Additività: se i dati sono indipendenti, la log-verosimiglianza totale è la somma delle log-verosimiglianze individuali:

\[ \ell(\theta; y_1, \dots, y_n) = \sum_{i=1}^{n} \log f(y_i \mid \theta). \]

Questa proprietà è fondamentale nei modelli statistici in cui si assumono osservazioni indipendenti e identicamente distribuite (i.i.d.).

Infine, è importante ricordare che massimizzare la log-verosimiglianza è equivalente a massimizzare la verosimiglianza, poiché il logaritmo è una funzione monotona crescente. Per questo motivo, molte tecniche di stima e modellazione in statistica e machine learning sono formulate in termini di log-verosimiglianza.

42.3 Sequenza di Lanci di una Moneta

Per comprendere il concetto di verosimiglianza, iniziamo con un esempio semplice e intuitivo. Immagina di voler stimare la probabilità che una moneta cada su testa, e chiamiamo questa probabilità \(p_H\).

Il nostro obiettivo è capire quali valori di \(p_H\) rendono più plausibile ciò che abbiamo osservato. Per farlo, vedremo come calcolare la probabilità di osservare determinate sequenze di lanci, assumendo che ogni lancio sia indipendente dagli altri.

42.3.1 Perché moltiplichiamo le probabilità?

Se lanciamo una moneta più volte, ogni risultato è indipendente dai precedenti: il fatto che esca “testa” o “croce” in un lancio non influenza il risultato del lancio successivo. Questo significa che la probabilità di osservare una specifica sequenza si ottiene moltiplicando le probabilità di ciascun singolo lancio.

Ad esempio, supponiamo che la probabilità di testa sia \(p_H\), e quindi la probabilità di croce sarà \(1 - p_H\). Se osserviamo la sequenza HTHT, allora la probabilità di ottenere questa sequenza è:

\[ P(\text{HTHT} \mid p_H) = p_H \cdot (1 - p_H) \cdot p_H \cdot (1 - p_H) = p_H^2 (1 - p_H)^2. \]

Se invece osserviamo una sequenza come THTH o HHTT, la probabilità è la stessa: \(p_H^2 (1 - p_H)^2.\) Questo perché non ci interessa l’ordine dei lanci, ma solo quante teste e quante croci abbiamo ottenuto.

42.3.2 Generalizzazione

In generale, se lanciamo una moneta \(n\) volte e osserviamo \(y\) teste (successi) e \(n - y\) croci (insuccessi), la probabilità di ottenere una qualsiasi sequenza con questa configurazione è:

\[ P(Y = y \mid p_H) = p_H^y (1 - p_H)^{n - y}. \tag{42.1}\]

Questa espressione è la funzione di verosimiglianza, perché ci dice quanto un certo valore di \(p_H\) è compatibile con i dati osservati.

42.3.3 Connessione con la Distribuzione Binomiale

L’Equazione 42.1 è anche il nucleo della distribuzione binomiale. La formula completa per la probabilità binomiale è:

\[ P(Y = y \mid p_H) = \binom{n}{y} p_H^y (1 - p_H)^{n - y}. \]

L’unica differenza è il coefficiente binomiale \(\binom{n}{y}\), che conta quante sequenze diverse portano allo stesso numero di successi e insuccessi.

Quando calcoliamo la verosimiglianza, possiamo ignorare questo coefficiente: non dipende da \(p_H\), quindi non influisce sul valore che massimizza la funzione. Per questo motivo, usiamo solo il “nucleo” della funzione binomiale:

\[ L(p_H \mid y) = p_H^y (1 - p_H)^{n - y}. \tag{42.2}\]

42.3.4 Esempio 1: Due lanci

Immaginiamo di lanciare una moneta due volte e di osservare una testa e una croce. Il nostro obiettivo è stimare \(p_H\), cioè la probabilità che la moneta cada su testa. Per iniziare, valutiamo quanto sono plausibili due diversi valori di \(p_H\) alla luce dei dati osservati, calcolando la funzione di verosimiglianza.

Se \(p_H = 0.5\) (una moneta equa), allora:

\[ L(0.5) = 0.5^1 \cdot (1 - 0.5)^1 = 0.25 . \]
Se \(p_H = 0.4\), allora:

\[ L(0.4) = 0.4^1 \cdot 0.6^1 = 0.24 . \]

In entrambi i casi, la funzione di verosimiglianza ci dice quanto sia plausibile quel valore di \(p_H\), dati i risultati che abbiamo osservato (una testa e una croce). Come possiamo notare, il valore \(p_H = 0.5\) è più compatibile con i dati osservati rispetto a \(p_H = 0.4\).

42.3.5 Calcolo della verosimiglianza su una griglia di valori

Ora ripetiamo questo calcolo per molti valori diversi di \(p_H\), compresi tra 0 e 1. Questo ci permette di visualizzare la funzione di verosimiglianza e di vedere per quali valori di \(p_H\) essa è più alta.

# Parametri osservati
n <- 2             # Numero totale di lanci
y <- 1             # Numero di teste osservate

# Griglia di valori per p_H da 0 a 1
p_H <- seq(0, 1, length.out = 100)

# Calcolo della funzione di verosimiglianza
likelihood <- p_H^y * (1 - p_H)^(n - y)

42.3.6 Rappresentazione grafica

Il grafico seguente mostra la funzione di verosimiglianza per i 100 valori di \(p_H\). Ogni punto della curva indica quanto è compatibile quel valore di \(p_H\) con i dati che abbiamo osservato.

ggplot(data.frame(p_H, likelihood), aes(x = p_H, y = likelihood)) +
  geom_line(color = okabe_ito_palette[2], linewidth = 1.2) +
  labs(
    title = "Funzione di Verosimiglianza per 2 Lanci\n(1 Testa, 1 Croce)",
    x = expression(p[H]),
    y = "Verosimiglianza"
  )

42.3.7 Cosa ci dice il grafico?

La funzione di verosimiglianza raggiunge il suo massimo per \(p_H = 0.5\), che corrisponde alla proporzione osservata (1 testa su 2 lanci).
Questo valore di \(p_H\) è quindi il più plausibile secondo i dati: rappresenta la stima di massima verosimiglianza.
I valori estremi (vicini a 0 o a 1) hanno verosimiglianza molto bassa: non spiegano bene il fatto che abbiamo ottenuto una testa e una croce.

Ecco una versione migliorata e più didattica del tuo Esempio 2: Tre lanci, pensata per studenti che stanno imparando il concetto di verosimiglianza passo dopo passo. Mantiene l’approccio del primo esempio e rafforza la continuità logica, l’intuizione e l’interpretazione dei risultati, con un linguaggio semplice ma preciso.

42.3.8 Esempio 2: Tre lanci

Proseguiamo con un secondo esperimento: lanciamo una moneta tre volte, e otteniamo una testa e due croci. Anche in questo caso, vogliamo stimare \(p_H\), la probabilità che la moneta cada su testa, e valutare quali valori di \(p_H\) sono più compatibili con i dati osservati.

Iniziamo calcolando la verosimiglianza per due valori specifici:

Se \(p_H = 0.5\) (moneta equa): \[ L(0.5) = 0.5^1 \cdot (1 - 0.5)^2 = 0.5 \cdot 0.25 = 0.125 \]
Se \(p_H = 0.4\): \[ L(0.4) = 0.4^1 \cdot 0.6^2 = 0.4 \cdot 0.36 = 0.144 \]

Come vediamo, in questo caso \(p_H = 0.4\) ha una verosimiglianza maggiore rispetto a \(p_H = 0.5\): questo suggerisce che il valore 0.4 spiega meglio i dati (1 testa su 3 lanci) rispetto alla moneta equa.

42.3.9 Calcolo su una griglia di valori

Calcoliamo ora la verosimiglianza per 100 valori compresi tra 0 e 1 per visualizzare la funzione su tutto l’intervallo di probabilità.

# Parametri osservati
n <- 3             # Numero totale di lanci
y <- 1             # Numero di teste osservate

# Sequenza di valori possibili per p_H
p_H <- seq(0, 1, length.out = 100)

# Calcolo della funzione di verosimiglianza
likelihood <- p_H^y * (1 - p_H)^(n - y)

42.3.10 Rappresentazione grafica

Il grafico seguente mostra la funzione di verosimiglianza: ci dice quanto ogni valore di \(p_H\) è compatibile con l’osservazione di 1 testa e 2 croci.

# Mostra la curva della verosimiglianza
ggplot(data.frame(p_H, likelihood), aes(x = p_H, y = likelihood)) +
  geom_line(color = okabe_ito_palette[2], linewidth = 1.2) +
  labs(
    title = "Funzione di Verosimiglianza per 3 Lanci\n(1 Testa, 2 Croci)",
    x = expression(p[H]),
    y = "Verosimiglianza"
  )

42.3.11 Cosa osserviamo?

Il massimo della funzione di verosimiglianza non è più in \(p_H = 0.5\), ma più vicino a 0.33, cioè alla proporzione osservata (1 testa su 3 lanci).
Questo riflette il fatto che il valore di \(p_H\) che meglio spiega i dati è quello che riproduce la frequenza osservata.
La curva è più stretta rispetto al caso con 2 lanci: abbiamo più informazioni, quindi possiamo stimare \(p_H\) con maggiore precisione.
Anche qui, i valori estremi di \(p_H\) (vicino a 0 o 1) hanno verosimiglianze basse, perché non giustificano bene l’osservazione di una testa e due croci.

42.3.12 Interpretazione dei Risultati

La funzione di verosimiglianza raggiunge il massimo per valori di \(p_H\) vicini alla proporzione di teste osservata.
Quando il numero di lanci aumenta, la curva diventa più “stretta”: i dati ci permettono di stimare \(p_H\) in modo più preciso.
Valori estremi di \(p_H\) (vicini a 0 o 1) hanno verosimiglianze basse: non spiegano bene i dati osservati.

In sintesi, abbiamo visto come costruire la funzione di verosimiglianza a partire dai dati osservati, senza usare direttamente la distribuzione binomiale completa.
La tua sezione è già molto chiara, ma possiamo migliorarla ulteriormente per renderla ancora più didattica e accessibile a studenti principianti, mantenendo un linguaggio preciso e coerente con le sezioni precedenti.
L’obiettivo è aiutare gli studenti a vedere come la distribuzione binomiale si collega alla verosimiglianza e alla stima del parametro.

42.4 Verosimiglianza binomiale

Torniamo al nostro esperimento con la moneta, ma questa volta lo rendiamo un po’ più realistico: lanciamo la moneta \(n = 30\) volte e otteniamo \(y = 23\) teste.

Per modellare il numero totale di teste osservate, utilizziamo la distribuzione binomiale, che ci dice qual è la probabilità di ottenere esattamente \(y\) successi su \(n\) prove, quando la probabilità di successo in ogni singolo lancio è \(\theta\):

\[ P(Y = y \mid \theta) = \binom{n}{y} \theta^y (1 - \theta)^{n - y}. \]

In questo contesto:

\(Y\) è una variabile casuale che rappresenta il numero di teste ottenute,
\(y = 23\) è il valore osservato,
\(\theta\) (o \(p_H\)) è la probabilità incognita di ottenere testa in un singolo lancio.

42.4.1 Dalla distribuzione alla verosimiglianza

Una volta osservati i dati (\(y = 23\)), possiamo considerarli fissati e analizzare quanto ciascun valore possibile del parametro \(\theta\) (la probabilità di ottenere testa) sia compatibile con questi dati. Per farlo, possiamo riutilizzare direttamente la formula della distribuzione binomiale, trattandola come funzione di \(\theta\) anziché come funzione di \(y\):

\[ L(\theta \mid y = 23) = \binom{30}{23} \theta^{23} (1 - \theta)^7. \]

Questa è la funzione di verosimiglianza, che ci dice quanto ogni valore di \(\theta\) sia plausibile alla luce dei dati osservati. A differenza degli esempi precedenti (in cui abbiamo ignorato le costanti moltiplicative), qui non abbiamo bisogno di semplificare la formula: la funzione dbinom() in R calcola automaticamente l’intera espressione, costante inclusa.

42.4.2 Visualizzazione della funzione di verosimiglianza

Il codice seguente costruisce il grafico della funzione di verosimiglianza, calcolando per ogni valore di \(\theta\) la probabilità di osservare 23 successi su 30 prove:

# Parametri osservati
n <- 30       # Numero totale di lanci
y <- 23       # Numero di teste osservate

# Griglia di valori possibili per p_H
p_H <- seq(0, 1, length.out = 1000)

# Calcolo della funzione di verosimiglianza (usando la binomiale completa)
likelihood <- dbinom(y, size = n, prob = p_H)

# Creazione del data frame
data <- data.frame(p_H, likelihood)

# Grafico della funzione di verosimiglianza
ggplot(data, aes(x = p_H, y = likelihood)) +
  geom_line(color = okabe_ito_palette[2], linewidth = 1.2) +
  labs(
    title = "Funzione di Verosimiglianza per 30 Lanci di Moneta (23 Teste)",
    x = expression(p[H]),
    y = "Verosimiglianza"
  )

42.4.3 Interpretazione del risultato

Osservando il grafico, vediamo che:

la funzione di verosimiglianza raggiunge il suo massimo in corrispondenza di \(p_H \approx 0.77\);
questo significa che il valore di \(\theta\) (cioè \(p_H\)) che rende più plausibile l’osservazione di 23 teste su 30 è circa 0.77;
in altre parole, la stima di massima verosimiglianza (MLE) per la probabilità di ottenere testa è: \[ \hat{p}_H = \frac{23}{30} \approx 0.767. \]

Questo risultato è del tutto intuitivo: la stima migliore è la proporzione osservata di teste.
La verosimiglianza ci mostra quali altri valori di \(\theta\) sono meno compatibili con i dati.

In sintesi, abbiamo visto come utilizzare direttamente la distribuzione binomiale per costruire la funzione di verosimiglianza. Questa funzione è uno strumento fondamentale per confrontare diversi valori possibili del parametro \(\theta\) (cioè la probabilità di successo) alla luce dei dati osservati.

Nel caso della moneta lanciata 30 volte con 23 teste, abbiamo trattato il numero di successi osservati (\(y = 23\)) come un dato fisso, e abbiamo valutato per quali valori di \(\theta\) questa osservazione sarebbe più plausibile.

In R, per calcolare la funzione di verosimiglianza non serve riscrivere manualmente la formula della binomiale. Possiamo usare la funzione dbinom(), che calcola le probabilità (o masse) della distribuzione binomiale per un dato numero di successi \(y\), un numero totale di prove \(n\), e una certa probabilità di successo \(\theta\).

Ad esempio, nel codice:

likelihood <- dbinom(y, size = n, prob = p_H)

y è il numero di successi osservati (23),
n è il numero totale di prove (30),
p_H è un vettore di valori di \(\theta\) tra 0 e 1.

Il risultato likelihood è un vettore che contiene, per ciascun valore di \(\theta\), la verosimiglianza associata a quel valore: cioè, quanto quel valore di \(\theta\) è compatibile con i dati che abbiamo osservato.

Riassumendo:

dbinom() fornisce la funzione di probabilità della binomiale,
fissando \(y\) e variando \(\theta\), usiamo dbinom() per costruire la funzione di verosimiglianza,
possiamo poi tracciare un grafico di questa funzione per visualizzare quali valori di \(\theta\) sono più plausibili.

Questa strategia mostra come la verosimiglianza possa essere costruita direttamente a partire da una distribuzione conosciuta (in questo caso, la binomiale) e implementata facilmente in R con strumenti standard.

42.5 La Stima di Massima Verosimiglianza

Nel momento in cui osserviamo dei dati e vogliamo stimare un parametro incognito (ad esempio, la probabilità \(\theta\) che una moneta dia “testa”), un metodo classico è la stima di massima verosimiglianza (Maximum Likelihood Estimation, o MLE).

Anche se nella prospettiva bayesiana ci interessa la distribuzione completa dei valori plausibili del parametro (e non un singolo valore stimato), è utile conoscere il concetto di MLE, che rappresenta il valore di \(\theta\) che rende i dati osservati più compatibili con il modello. Più avanti vedremo come la MLE corrisponde, nel caso di una prior uniforme, al valore massimo della distribuzione a posteriori.

42.5.1 L’idea di Fondo

La logica è semplice. Una volta osservati i dati, ci chiediamo: quali valori del parametro sono più compatibili con ciò che abbiamo visto? Il valore che risulta più plausibile è la stima di massima verosimiglianza.

Immagina di “provare” tanti valori di \(\theta\), e per ciascuno chiederti: “se fosse questo il vero valore di \(\theta\), quanto sarebbe plausibile osservare questi dati?” Il valore che meglio spiega i dati è quello scelto come stima.

42.5.2 Un Esempio Concreto

Hai lanciato una moneta 30 volte e hai osservato 23 teste. La tua domanda è: quanto è sbilanciata la moneta?

Utilizziamo la funzione di verosimiglianza, che ci dice quanto ogni valore possibile di \(\theta\) è compatibile con i dati osservati. Più alta è la verosimiglianza, più plausibile è il valore di \(\theta\).

42.5.3 Una metafora visiva

Immagina di osservare una curva che rappresenta la verosimiglianza al variare del parametro \(\theta\).
La curva si alza, raggiunge un punto massimo, e poi scende: proprio come una collina.

Quel punto più alto della curva rappresenta il valore di \(\theta\) che è più compatibile con i dati osservati: è il valore per cui i dati risultano più verosimili, secondo il modello.

Quel valore è chiamato stima di massima verosimiglianza (maximum likelihood estimate): è il punto in cima alla collina della verosimiglianza.

42.5.4 Come si Trova il Massimo?

Dal punto di vista matematico, il massimo di una funzione si trova dove la sua pendenza (la derivata) è zero: cioè, dove smette di salire e inizia a scendere.

Per la verosimiglianza binomiale (trasformata in log-verosimiglianza), questo calcolo si può fare in modo esatto. Se hai osservato \(y\) successi su \(n\) prove, la log-verosimiglianza è:

\[ \ell(\theta) = y \log \theta + (n - y) \log(1 - \theta), \]

e la derivata si annulla quando:

\[ \hat{\theta} = \frac{y}{n}. \]

In altre parole, la MLE è semplicemente la proporzione osservata di successi.

42.5.5 Un punto Importante per l’Approccio Bayesiano

Nel contesto bayesiano, non ci limitiamo a cercare il valore più plausibile del parametro, ma costruiamo una distribuzione completa che rappresenta l’incertezza su \(\theta\). Tuttavia, il punto in cui la distribuzione a posteriori raggiunge il suo massimo viene chiamato MAP (Maximum A Posteriori), e se assumiamo una distribuzione a priori uniforme, MLE e MAP coincidono. Quindi, la MLE può essere vista come un caso speciale del ragionamento bayesiano, utile per introdurre il concetto di compatibilità tra parametri e dati.

42.6 Calcolare la MLE in R

42.6.1 Metodo 1: Valutazione su Griglia

# Parametri
n <- 30
y <- 23
theta <- seq(0, 1, length.out = 10000)

# Calcolo della verosimiglianza binomiale
likelihood <- dbinom(y, size = n, prob = theta)

# Trova il massimo
max_index <- which.max(likelihood)
optimal_theta <- theta[max_index]

# Risultato
optimal_theta
#> [1] 0.7667

Spiegazione:

dbinom() calcola la verosimiglianza per ogni valore di \(\theta\);
which.max() individua il massimo;
theta[max_index] restituisce la stima MLE.

42.6.2 Metodo 2: Ottimizzazione Numerica

In alternativa, possiamo trovare la MLE senza usare griglie, con un approccio più efficiente.

# Funzione log-verosimiglianza negativa
neg_log_likelihood <- function(theta) {
  - (y * log(theta) + (n - y) * log(1 - theta))
}

# Ottimizzazione numerica
result <- optim(
  par = 0.5,
  fn = neg_log_likelihood,
  method = "Brent",
  lower = 1e-6,
  upper = 1 - 1e-6
)

optimal_theta_numerical <- result$par
optimal_theta_numerical
#> [1] 0.7667

Nota: abbiamo calcolato la log-verosimiglianza negativa, perché optim() cerca minimi per default.

42.6.3 Confronto tra le Soluzioni

c(
  "Griglia" = optimal_theta, 
  "Ottimizzazione" = optimal_theta_numerical, 
  "Analitica" = y / n
)
#>        Griglia Ottimizzazione      Analitica 
#>         0.7667         0.7667         0.7667

Tutti i metodi restituiscono lo stesso risultato:
\[ \hat{\theta} = \frac{23}{30} \approx 0.767. \]

42.7 Verosimiglianza Congiunta

Abbiamo visto che, nel caso di una sequenza di \(n\) lanci di una moneta, la funzione di verosimiglianza si basa sulla distribuzione binomiale. In questo caso, trattiamo un esperimento Bernoulliano ripetuto \(n\) volte, e la nostra osservazione è il numero totale di successi (teste). Il numero complessivo di successi segue una distribuzione binomiale, e la funzione di verosimiglianza assume la forma:

\[ \mathcal{L}(\theta) = P(Y = y \mid \theta) = \binom{n}{y} \theta^y (1 - \theta)^{n - y}. \]

Qui la verosimiglianza è espressa direttamente in termini del numero totale di successi e insuccessi, senza dover scrivere il contributo di ogni singolo lancio.

Tuttavia, possiamo affrontare la questione da una prospettiva diversa: invece di considerare il numero totale di successi, possiamo pensare alla verosimiglianza come il prodotto delle probabilità di ogni singolo lancio. Questo ci porta a una generalizzazione importante: la verosimiglianza congiunta di più osservazioni indipendenti.

42.7.1 Dal Caso Binomiale alla Verosimiglianza Congiunta

Nel caso dei lanci della moneta, le singole osservazioni sono prove Bernoulliane indipendenti, ovvero ogni singolo lancio è un’osservazione indipendente che segue una distribuzione Bernoulli con parametro \(\theta\):

\[ P(Y_i = 1 \mid \theta) = \theta, \quad P(Y_i = 0 \mid \theta) = 1 - \theta. \]

Se trattiamo ogni prova individualmente, la funzione di verosimiglianza per una singola osservazione è:

\[ \mathcal{L}(\theta \mid y_i) = \theta^{y_i} (1 - \theta)^{1 - y_i}. \]

Ora, per un campione di \(n\) osservazioni indipendenti, la verosimiglianza congiunta è il prodotto delle verosimiglianze delle singole osservazioni:

\[ \mathcal{L}(\theta \mid y_1, y_2, \dots, y_n) = \prod_{i=1}^{n} \theta^{y_i} (1 - \theta)^{1 - y_i}. \]

Riconosciamo che questa espressione è identica alla funzione di verosimiglianza della distribuzione binomiale, perché il numero totale di successi è:

\[ y = \sum_{i=1}^{n} y_i. \]

Quindi, riscrivendo la verosimiglianza congiunta, otteniamo:

\[ \mathcal{L}(\theta) = \theta^{\sum y_i} (1 - \theta)^{n - \sum y_i} = \theta^y (1 - \theta)^{n - y}. \]

Questa è proprio la verosimiglianza della distribuzione binomiale! Questo mostra che il caso binomiale può essere visto come una forma compatta di verosimiglianza congiunta per prove Bernoulliane indipendenti.

42.8 Perché è Importante la Verosimiglianza Congiunta?

L’idea della verosimiglianza congiunta è fondamentale perché ci permette di estendere i concetti di verosimiglianza dal caso di una singola osservazione al caso di molte osservazioni indipendenti. Questo è utile in molti contesti statistici:

stimare parametri basandosi su un intero campione invece che su una singola osservazione;
definire modelli statistici più complessi, in cui le osservazioni sono indipendenti ma possono avere diverse distribuzioni.

In sintesi, l’esempio binomiale mostra che la verosimiglianza congiunta di prove Bernoulliane indipendenti si riconduce alla verosimiglianza binomiale. Tuttavia, la vera potenza della verosimiglianza congiunta si manifesta in distribuzioni continue come la normale, dove la produttoria delle densità di probabilità per singole osservazioni è chiaramente distinta dalla funzione di verosimiglianza per il campione intero.

La chiave per comprendere il concetto è rendersi conto che la verosimiglianza di un’intera sequenza di prove indipendenti è il prodotto delle singole verosimiglianze.

42.9 Esempio: Osservazioni Raggruppate

Per illustrare il concetto di verosimiglianza congiunta nel caso della distribuzione binomiale, consideriamo quattro gruppi distinti di osservazioni binomiali indipendenti:

gruppo 1: 23 successi su 30 prove,
gruppo 2: 20 successi su 28 prove,
gruppo 3: 29 successi su 40 prove,
gruppo 4: 29 successi su 36 prove.

Poiché ogni gruppo segue una distribuzione binomiale indipendente con lo stesso parametro \(\theta\), la log-verosimiglianza congiunta si ottiene sommando le log-verosimiglianze di ciascun gruppo:

\[ \log \mathcal{L}(\theta) = \sum_{i=1}^{4} \left[ y_i \log(\theta) + (n_i - y_i) \log(1 - \theta) \right], \]

dove:

\(n_i\) è il numero totale di prove nel gruppo \(i\),
\(y_i\) è il numero di successi nel gruppo \(i\).

Sostituendo i valori specifici:

\[ \begin{aligned} \log \mathcal{L}(\theta) &= [23\log(\theta) + (30-23)\log(1 - \theta)] \\ &\quad + [20\log(\theta) + (28-20)\log(1 - \theta)] \\ &\quad + [29\log(\theta) + (40-29)\log(1 - \theta)] \\ &\quad + [29\log(\theta) + (36-29)\log(1 - \theta)]. \end{aligned} \]

Questa formula ci permette di calcolare quanto è plausibile il valore del parametro \(\theta\), tenendo conto simultaneamente di tutte le osservazioni nei quattro gruppi.

42.10 Implementazione in R

Supponiamo di avere i dati di quattro gruppi indipendenti, e vogliamo trovare la stima del parametro \(\theta\) che massimizza la log-verosimiglianza congiunta. Procederemo passo dopo passo.

1. Definire una funzione per la log-verosimiglianza congiunta.

Questa funzione riceve un valore di \(\theta\) e una lista di gruppi. Ogni gruppo contiene il numero totale di prove e il numero di successi. La funzione calcola la somma delle log-verosimiglianze per ciascun gruppo.

# Funzione che calcola la log-verosimiglianza congiunta
log_verosimiglianza_congiunta <- function(theta, dati) {
  # Limitiamo theta per evitare log(0)
  if (theta <= 0) theta <- 1e-10
  if (theta >= 1) theta <- 1 - 1e-10

  # Inizializziamo il totale
  somma_loglik <- 0

  # Per ogni gruppo, calcoliamo il contributo alla log-verosimiglianza
  for (i in 1:length(dati)) {
    gruppo <- dati[[i]]   # Estrae il gruppo i-esimo
    n <- gruppo[1]        # Numero totale di prove
    y <- gruppo[2]        # Numero di successi

    # Aggiunge il contributo del gruppo alla somma totale
    somma_loglik <- somma_loglik + y * log(theta) + (n - y) * log(1 - theta)
  }

  # Restituiamo il valore negativo (perché optim() cerca minimi)
  return(-somma_loglik)
}

2. Inserire i dati.

Qui definiamo i dati per ciascun gruppo come coppie (numero di prove, numero di successi):

# Dati osservati per ciascun gruppo
dati_gruppi <- list(
  c(30, 23),
  c(28, 20),
  c(40, 29),
  c(36, 29)
)

3. Trovare la stima di θ che massimizza la verosimiglianza.

Usiamo optim() per trovare numericamente il valore di \(\theta\) che minimizza il valore negativo della log-verosimiglianza, ovvero che massimizza la log-verosimiglianza.

# Ricerca del valore ottimale di theta
result <- optim(
  par = 0.5,                            # Valore iniziale
  fn = log_verosimiglianza_congiunta,  # Funzione da minimizzare
  dati = dati_gruppi,                  # Passiamo i dati
  method = "L-BFGS-B",                 # Metodo con vincoli
  lower = 0, upper = 1                 # Vincoli: theta tra 0 e 1
)

# Stima ottimale trovata
result$par
#> [1] 0.7537

4. Visualizzare la log-verosimiglianza.

Costruiamo ora un grafico che mostri come varia la log-verosimiglianza al variare di \(\theta\).

# Valori di theta da esplorare
theta_values <- seq(0.01, 0.99, length.out = 100)

# Vettore per salvare i valori di log-verosimiglianza
log_likelihood_values <- numeric(length(theta_values))

# Calcoliamo il valore della funzione per ogni valore di theta
for (i in 1:length(theta_values)) {
  t <- theta_values[i]
  log_likelihood_values[i] <- log_verosimiglianza_congiunta(t, dati_gruppi)
}

Ora costruiamo il grafico:

ggplot(
  data.frame(theta = theta_values, log_likelihood = log_likelihood_values),
  aes(x = theta, y = log_likelihood)
) +
  geom_line(color = "blue", linewidth = 1.2) +
  labs(
    title = "Funzione di Log-Verosimiglianza Congiunta",
    x = expression(theta),
    y = "Log-verosimiglianza negativa"
  )

Con questo esempio abbiamo visto che:

la log-verosimiglianza congiunta si ottiene sommando le log-verosimiglianze dei singoli gruppi;
è possibile trovare la stima ottimale di \(\theta\) numericamente, senza formule complicate;
il grafico della log-verosimiglianza ci aiuta a visualizzare il punto in cui il modello è più compatibile con i dati.

Questo approccio — basato sull’uso della verosimiglianza e della sua somma tra gruppi indipendenti — sarà anche la base per il passaggio all’inferenza bayesiana, dove aggiungeremo una distribuzione a priori per ottenere una distribuzione a posteriori di \(\theta\).

42.11 La Verosimiglianza Marginale

La verosimiglianza marginale è un concetto fondamentale nell’inferenza bayesiana, utilizzato per valutare quanto un modello sia compatibile con i dati osservati, tenendo conto dell’incertezza sui parametri.

A differenza della verosimiglianza standard, che misura la plausibilità dei dati per un valore fisso del parametro, la verosimiglianza marginale considera tutti i possibili valori del parametro, pesandoli in base alla loro probabilità a priori. Questo approccio permette di integrare l’incertezza nella valutazione del modello.

42.11.1 Caso con Parametri Discreti

Per comprendere meglio il concetto, consideriamo un esperimento in cui eseguiamo 10 tentativi e otteniamo 7 successi. Supponiamo che la probabilità di successo \(\theta\) possa assumere solo tre valori discreti:

\[ \theta \in \{0.1, 0.5, 0.9\}. \]

Per calcolare la verosimiglianza marginale, dobbiamo:

Assegnare una probabilità a priori a ciascun valore di \(\theta\), ad esempio:
- Distribuzione uniforme: \[ p(\theta = 0.1) = p(\theta = 0.5) = p(\theta = 0.9) = \frac{1}{3}. \]
- Distribuzione non uniforme (ad esempio, dando più peso a \(\theta = 0.5\)): \[ p(\theta = 0.1) = \frac{1}{4}, \quad p(\theta = 0.5) = \frac{1}{2}, \quad p(\theta = 0.9) = \frac{1}{4}. \]
Calcolare la probabilità di osservare 7 successi su 10 prove per ogni valore di \(\theta\):

\[ p(k=7 \mid \theta) = \binom{10}{7} \theta^7 (1 - \theta)^3. \]
Moltiplicare ciascuna di queste probabilità per la corrispondente probabilità a priori e sommare i risultati:

\[ p(k=7 \mid n=10) = \sum_{i} p(k=7 \mid \theta_i) p(\theta_i). \]

Sostituendo i valori per la distribuzione uniforme:

\[ p(k=7 \mid n=10) = \binom{10}{7} 0.1^7 (0.9)^3 \cdot \frac{1}{3} + \binom{10}{7} 0.5^7 (0.5)^3 \cdot \frac{1}{3} + \binom{10}{7} 0.9^7 (0.1)^3 \cdot \frac{1}{3}. \]

Questa somma rappresenta la verosimiglianza marginale, ossia la probabilità di ottenere 7 successi su 10, considerando tutte le possibili incertezze su \(\theta\).

42.11.2 Caso con Parametri Continui

Nella maggior parte delle situazioni, il parametro \(\theta\) non assume solo pochi valori discreti, ma può variare continuamente in un intervallo (ad esempio, tra 0 e 1). In questo caso, invece di sommare, dobbiamo integrare:

\[ p(k=7 \mid n=10) = \int_{0}^{1} \binom{10}{7} \theta^7 (1 - \theta)^3 p(\theta) \, d\theta. \]

Qui:

\(p(\theta)\) è la distribuzione a priori di \(\theta\).
L’integrale rappresenta una media ponderata della probabilità di ottenere i dati, considerando tutti i valori di \(\theta\).

Ad esempio, se \(\theta \sim \text{Beta}(2,2)\), la verosimiglianza marginale diventa:

\[ p(k=7 \mid n=10) = \int_{0}^{1} \binom{10}{7} \theta^7 (1 - \theta)^3 \frac{\theta (1-\theta)}{B(2,2)} \, d\theta. \]

Questo tipo di calcolo viene spesso risolto numericamente.

42.11.3 Calcolo Numerico della Verosimiglianza Marginale in R

Se vogliamo calcolare la verosimiglianza marginale numericamente, possiamo usare l’integrazione numerica in R.

Caso con Parametri Discreti.

# Definiamo i valori possibili di theta e le probabilità a priori
theta_vals <- c(0.1, 0.5, 0.9)
prior_probs <- c(1/3, 1/3, 1/3)  # Distribuzione uniforme

# Calcoliamo la verosimiglianza per ciascun valore di theta
likelihoods <- dbinom(7, size = 10, prob = theta_vals)

# Calcoliamo la verosimiglianza marginale sommando i contributi ponderati
marginal_likelihood <- sum(likelihoods * prior_probs)
print(marginal_likelihood)
#> [1] 0.0582

Caso con Parametri Continui.

# Definiamo la funzione di verosimiglianza pesata dalla prior
integrand <- function(theta) {
  dbinom(7, size = 10, prob = theta) * dbeta(theta, shape1 = 2, shape2 = 2)
}

# Eseguiamo l'integrazione numerica
marginal_likelihood <- integrate(integrand, lower = 0, upper = 1)$value
print(marginal_likelihood)
#> [1] 0.1119

42.11.4 Interpretazione della Verosimiglianza Marginale

La verosimiglianza marginale rappresenta la probabilità complessiva dei dati, tenendo conto di tutte le possibili incertezze sul parametro \(\theta\).

Se la verosimiglianza marginale è alta, significa che il modello nel suo insieme è compatibile con i dati osservati.
Se la verosimiglianza marginale è bassa, significa che, indipendentemente dal valore di \(\theta\), il modello non spiega bene i dati.

A differenza della verosimiglianza classica, che valuta quanto siano probabili i dati per un singolo valore di \(\theta\), la verosimiglianza marginale considera tutte le possibili ipotesi sul parametro.

42.11.5 Ruolo nella Statistica Bayesiana

La verosimiglianza marginale svolge un ruolo cruciale nell’inferenza bayesiana perché appare nella formula di Bayes:

\[ p(\theta \mid D) = \frac{p(D \mid \theta) p(\theta)}{p(D)}, \]

dove \(p(D)\) è la verosimiglianza marginale. Questa quantità:

serve da fattore di normalizzazione per la distribuzione a posteriori \(p(\theta \mid D)\);
permette di confrontare modelli diversi, attraverso il fattore di Bayes:

\[ BF = \frac{p(D \mid M_1)}{p(D \mid M_2)}, \]

dove \(M_1\) e \(M_2\) sono due modelli diversi.

In conclusione, la verosimiglianza marginale è un concetto chiave nell’inferenza bayesiana, che ci permette di valutare quanto un modello sia coerente con i dati, tenendo conto di tutte le possibili incertezze sui parametri:

per parametri discreti, si calcola come una somma ponderata;
per parametri continui, si calcola con un integrale;
è essenziale per il calcolo della distribuzione a posteriori e per il confronto tra modelli.

42.12 Riflessioni Conclusive

La funzione di verosimiglianza rappresenta un ponte fondamentale tra i dati osservati e i parametri di un modello statistico, offrendo una misura della plausibilità dei dati rispetto a diversi valori possibili dei parametri. La sua costruzione richiede l’integrazione di tre elementi chiave: il modello statistico ipotizzato come generatore dei dati, lo spazio dei parametri associato al modello e le osservazioni empiriche disponibili.

Nell’ambito dell’inferenza statistica, la funzione di verosimiglianza svolge un ruolo centrale. Essa consente di valutare quanto bene diversi valori dei parametri siano in grado di spiegare i dati osservati, diventando così uno strumento essenziale per la stima dei parametri e la selezione del modello. La sua corretta applicazione è determinante per garantire analisi dati rigorose e interpretazioni affidabili dei risultati.

In sintesi, padroneggiare il concetto di verosimiglianza e saperlo applicare in modo appropriato sono competenze indispensabili per chi si occupa di ricerca empirica e analisi di dati complessi. La verosimiglianza non è solo uno strumento tecnico, ma un pilastro metodologico che supporta la comprensione e l’interpretazione dei fenomeni osservati attraverso modelli statistici.

Esercizi

Problemi

Spiega ciascuno dei concetti seguenti con una frase:

probabilità.
funzione di massa di probabilità.
funzione di densità di probabilità.
distribuzione di probabilità.
distribuzione di probabilità discreta.
distribuzione di probabilità continua.
funzione di distribuzione cumulativa (cdf).
verosimiglianza

Problemi

All’esame ti verrà chiesto di:

Calcolare la funzione di verosimiglianza binomiale e riportare il valore della funzione in corrispondenza di specifici valori \(\theta\).
Calcolare la stima di massima verosimiglianza.
Rispondere a domande che implicano una adeguata comprensione del concetto di funzione di verosimiglianza.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.4
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
#> 
#> locale:
#> [1] C/UTF-8/C/C/C/C
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0 see_0.11.0      
#>  [5] gridExtra_2.3    patchwork_1.3.0  bayesplot_1.11.1 psych_2.5.3     
#>  [9] scales_1.3.0     markdown_2.0     knitr_1.50       lubridate_1.9.4 
#> [13] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.4     
#> [17] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1   
#> [21] tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.3    stringi_1.8.4     lattice_0.22-6    hms_1.1.3        
#>  [5] digest_0.6.37     magrittr_2.0.3    evaluate_1.0.3    grid_4.4.2       
#>  [9] timechange_0.3.0  fastmap_1.2.0     rprojroot_2.0.4   jsonlite_1.9.1   
#> [13] mnormt_2.1.1      cli_3.6.4         rlang_1.1.5       munsell_0.5.1    
#> [17] withr_3.0.2       tools_4.4.2       parallel_4.4.2    tzdb_0.5.0       
#> [21] colorspace_2.1-1  pacman_0.5.1      vctrs_0.6.5       R6_2.6.1         
#> [25] lifecycle_1.0.4   htmlwidgets_1.6.4 pkgconfig_2.0.3   pillar_1.10.1    
#> [29] gtable_0.3.6      glue_1.8.0        xfun_0.51         tidyselect_1.2.1 
#> [33] rstudioapi_0.17.1 farver_2.1.2      htmltools_0.5.8.1 nlme_3.1-167     
#> [37] labeling_0.4.3    rmarkdown_2.29    compiler_4.4.2

Bibliografia

Johnson, A. A., Ott, M., & Dogucu, M. (2022). Bayes Rules! An Introduction to Bayesian Modeling with R. CRC Press.

Schervish, M. J., & DeGroot, M. H. (2014). Probability and statistics (Vol. 563). Pearson Education London, UK: