3.8 Cenni ai metodi numerici

I risultati introdotti finora, riguardanti le variabili aleatorie, in particolare la formula di Bayes, permettono almeno in linea teorica di studiare in modo rigoroso moltissimi problemi concreti, in particolare qualora si possano ridurre a calcolare la densità (discreta o continua) di una variabile aleatoria \(Y\) (i “parametri” di un modello del problema), sulla base di informazione \(I\) inizialmente disponibile a cui si aggiunge un’informazione associata all’osservazione di una variabile \(X\) (i “dati” del modello), ovviamente rispetto alla quale \(Y\) non sia indipendente. Abbiamo infatti già visto diversi esempi in cui l’applicazione diretta della formula di Bayes (ad esempio, nel caso continuo) \[ p(Y = y | I, X=x) \propto p(Y=y|I) p( X=x| I, Y=y)\] permette di giustificare risultati riguardanti semplici situazioni, come ad esempio nel modello delle estrazioni dall’urna, in cui i “parametri” sono quantità relative allo stato dell’urna e i “dati” provengono dalle osservazioni delle estrazioni.

Dovrebbe anche essere evidente dagli stessi esempi che i metodi che abbiamo introdotto, sono generalizzabili in linea teorica a modelli arbitrariamente complessi: si pensi ad una rete bayesiana con tantissimi nodi e connessioni, in cui “dati” osservati \(X\) corrispondono alla variabile congiunta associata ad una famiglia di nodi, e si chiede di determinare la densità dei “parametri” \(Y\), la congiunta associata ai nodi rimanenti. Al crescere della numerosità dei dati e dei parametri, ci si scontra tuttavia con il problema di calcolare o almeno approssimare in modo computazionalmente efficiente una densità congiunta in uno spazio \(\R^d\) di dimensione estremamente elevata (diciamo dell’ordine dei parametri). Anche limitandosi alla stima di massimo a posteriori (o quella di massima verosimiglianza) \(y_{\mle}\), che fornisce comunque una informazione utile su \(Y\), si incontra il problema di massimizzare una funzione di \(y\), quindi di nuovo in uno spazio \(\R^d\) di dimensione molto elevata. Anche la numerosità dei dati \(X=x\) diventa problematica, ossia se \(X \in \R^D\) con \(D\) estremamente grande: diventa praticamente impossibile anche solo “definire” analiticamente la funzione di verosimiglianza \(p(X=x| I, Y=y)\), e quindi operare su di essa (ad esempio trovare la stima di massima verosimiglianza).

Questo problema era evidente ancor di più prima dell’uso dei computer, e classicamente era affrontato introducendo particolari densità a priori, in base alle specifiche funzioni di verosimiglianza in un problema, in modo che i calcoli delle densità a posteriori fossero analiticamente trattabili. Tali scelte speciali di densità a priori sono dette coniugate¹⁰ e sono ancora utili come particolari esempi e per situazioni semplici, ma limitano molto l’applicabilità del metodo. Inoltre non risolvono il problema quando comunque la numerosità delle osservazioni \(X\) diventa intrattabile.

Esempio 3.26 In \(n\) esperimenti indipendenti, ciascuno con probabilità di successo \(Y =y\in (0,1)\), il numero \(X\) di successi ha una densità binomiale di parametri \((n,y)\), pertanto se si è interessati a stimare \(Y\) dalle osservazioni di \(X\), la verosimiglianza è \[ L(Y = y; X =k) = P(X=k | Y = y) = {n \choose k} y^k (1-y)^{n-k}.\] Se la densità di \(Y\) a priori è della famiglia Beta di parametri \(\alpha\), \(\beta>0\), ossia \[ p( Y = y) \propto y^{\alpha-1} (1-y)^{\beta-1},\] (per \(y \in (0,1)\) e zero altrimenti) allora la densità a posteriori avendo osservato \(X=k\) successi è ancora dello stesso tipo, con parametri \(\alpha+k\), \(\beta+n-k\), \[ p(Y= y| X=k) \propto y^{\alpha+k-1} (1-y)^{\beta+n-k-1}.\]

Negli anni sono state introdotte e perfezionate molteplici tecniche numeriche che cercano di superare queste difficoltà, e la ricerca in questo ambito è estremamente attuale e rilevante. Volendo quindi accennare ad alcuni degli approcci principali, essi si possono dividere in due gruppi, secondo l’obiettivo che si pongono: approssimare tutta la densità a posteriori di \(Y\), oppure determinarne solamente la stima di massima verosimiglianza \(y_{\mle}\) (o di massimo a posteriori).

Nel primo caso, il problema consiste nell’approssimare una densità (di solito continua) in uno spazio di dimensione alta (anche la variante discreta comunque è importante, qualora si voglia ridurre la numerosità dei possibili valori). Negli esempi che abbiamo considerato, ci siamo limitati ad approssimare la densità valutandola in una griglia di punti equispaziati. È evidente tuttavia che questa scelta non è ottimale: vi sono regioni dove la densità è bassa e quindi sono poco rilevanti per l’approssimazione, mentre dove la densità è alta dovremmo viceversa infittire la griglia. In generale, il problema di approssimare una densità continua con un’opportuna densità discreta è detto di quantizzazione, e vi sono diversi algoritmi per ottenere soluzioni in modo efficace. Uno molto popolare nell’ambito della statistica bayesiana è il cosiddetto metodo Monte Carlo (più precisamente i metodi Markov chain Monte Carlo MCMC) in cui si “simulano” al computer un gran numero di variabili aleatorie indipendenti tutte con stessa densità (quella da approssimare). I valori ottenuti da tali simulazioni sostituiscono quindi la griglia di punti equispaziati. Questo fatto è una conseguenza dell’intepretazione della probabilità come frequenza (la legge dei grandi numeri), che dal punto di vista matematico è un teorema come vedremo nel Capitolo 8.

Remark. In R vi sono diverse librerie dedicate ai metodi Monte Carlo con applicazioni alla statistica bayesiana per l’approsimazione delle densità a posteriori. Due sono JAGS e Stan (che poi sono dei veri e propri linguaggi a sé, quindi non per ragioni di tempo non li illustreremo).

Nel secondo caso, ossia la determinazione della stima di massima verosimiglianza \(y_{\mle}\), si ricade in un contesto più ampio che è quello dell’ottimizzazione, ossia la determinazione di punti (e valori) di massimo o minimo di funzioni. Vi sono tantissimi algoritmi generali, molti dei quali sfruttano il calcolo in più variabili (nel caso appunto in cui \(Y \in \R^d\)), come ad esempio quelli di ascesa gradiente (o discesa se l’obiettivo è un minimo invece di un massimo). Qui l’idea di base è di avvicinarsi a \(y_{\mle}\) compiendo dei piccoli “passi” in ogni punto verso la direzione in cui la funzione aumenta maggiormente (appunto data dal gradiente della funzione nel punto). Il problema principale qui è che la “passeggiata” si potrebbe bloccare in un massimo locale (e non globale). Vi sono diversi metodi generali, spesso “probabilistici” come il simulated annealing che cercano di risolvere questa difficoltà.

In R le funzioni nlm() e optim() permettono di usare diversi metodi per l’ottimizzazione di funzioni (ossia la determinazione di minimi o massimi). Vediamo un esempio di applicazione.

Esempio 3.27 Consideriamo le densità dell’Esempio 3.20. Usando il comando optim() determiniamo numericamente la stima di massmia verosimiglianza per la frazione di palline rosse.

# Per determinare numericamente la
# stima di massima verosmiglianza
# definiamo come funzione da
# minimizzare l'opposto del logaritmo
# della verosimiglianza:

log_likelihood <- function(x) {
  -log(dbinom(3, 10, x))
}


# determiniamone il punto di minimo la
# funzione optim(), specificando un
# valore iniziale 1/2 per il metodo
# iterativo e un intervallo di valori
# deltax,1 -deltax.

deltax <- 0.01

plot(seq(0, 1, deltax), log_likelihood(seq(0,
  1, deltax)), col = miei_colori[1], lwd = 3,
  type = "l", xlab = "-log(L)", ylab = "")

x_mle <- optim(1/2, log_likelihood, method = "L-BFGS-B",
  lower = deltax, upper = 1 - deltax)

# il parametro trovato dal metodo è

x_mle$par

## [1] 0.3000006

# da confrontare con quello teorico
# 3/10

Concludiamo notando che vi sono anche casi di metodi “misti”. Spesso la variabile \(Y\) è interpretata come una variabile congiunta \[Y = (Y_{\operatorname{par}}, Y_{\operatorname{hid}}),\] in cui la prima marginale sono propriamente i “parametri” che si vogliono stimare, mentre la seconda sono variabili “nascoste” (in inglese latent variables) che non interessano direttamente (ma neppure sono i dati osservati \(X\)). La formula per la densità della marginale diventa in questo caso \[ p(Y_{\operatorname{par}} = y | I, X=x)= \int p(Y_{\operatorname{par}} = y | I, X=x, Y_{\operatorname{hid}} =z) p(Y_{\operatorname{hid}} = z | I, X=x) dz,\] dove l’integrale si estende a tutti i possibili valori delle variabili nascoste, possibilmente di dimensione molto grande. Inoltre la formula sopra richiede di calcolare (o almeno stimare) la densità a posteriori delle variabili nascoste. Per superare queste difficoltà è possibile limitarsi alla determinazione della stima di massima verosimiglianza per i parametri (quindi di \(Y_{\operatorname{par}}\)) utilizzando algoritmi che “alternano” l’uso di metodi del primo caso con metodi del secondo caso (il più popolare tra questi è l’algoritmo EM).

3.8.1 Esercizi

Esercizio 3.13 Usare la funzione rbinom() per simulare una variabile avente densità binomiale. Osservare che la frequenza dei valori osservati tende alla probabilità per un grande numero di simulazioni.

Esercizio 3.14 Usare la funzione optim() per calcolare numericamente \[ \arg \min \cur{ (x-1)^2 + (y-2)^2 - xy : (x,y ) \in \R^2}.\]