Apprendimento delle reti bayesiane - parametri e struttura

Apprendere una rete bayesiana (Reti bayesiane - semantica, indipendenze e inferenza esattaUna rete bayesiana è un grafo diretto aciclico (DAG) con una variabile per nodo e una tabella delle probabilità condizionate (CPT) $P(X_i\mid\text{Parents}(X_i))$ per nodo; rappresenta la congiunta come prodotto $P(x_1,\dots,x_n)=\prod_i P(x_i\mid\text{parents}(X_i))$, con molti meno numeri della tabella completa (rete dell'allarme: 10 contro 31). Tre strutture elementari: catena $X\to Y\to Z$ e fork $Y\leftarrow X\to Z$ (estremi indipendenti dato il centro), collider $X\to Z\leftarrow Y$ (indipendenti, ma dipendenti se si osserva $Z$ o un suo discendente). Modelli causali (genitori = cause) danno reti più piccole. Inferenza esatta per enumerazione: $\mathbf P(X\mid\mathbf e)=\alpha\sum_{\mathbf y}\mathbf P(X,\mathbf e,\mathbf y)$ calcolata come somma di prodotti, in profondità e da sinistra a destra; per l'allarme $\mathbf P(B\mid j,m)=\langle0{,}284;\ 0{,}716\rangle$. Molti prodotti si ripetono: l'eliminazione di variabili li salva come fattori e li riusa.Reti bayesiane - semantica, indipendenze e inferenza esatta →) comprende due compiti:

• apprendimento dei parametri (parameter learning): la struttura è nota, si stimano le CPT dai dati;
• apprendimento della struttura (structure learning): si ricava anche il grafo, con metodi globali (ricerca e punteggio) o locali (test di indipendenza condizionata).

Esempio di dati usato sotto: il golf delle slide, $14$ righe con attributi Outlook (sunny, overcast, rainy), Temperature (low, medium, high), Humidity (high, normal), Wind (yes, no) e classe Play ($P$ = si gioca, $N$ = no). La rete è un classificatore bayesiano ingenuo (naive Bayes): $\text{Play}$ è genitore di tutti gli attributi, che sono condizionalmente indipendenti dato $\text{Play}$ (Incertezza - probabilità, regola di Bayes e indipendenza condizionataUn agente con osservazioni parziali non può usare regole logiche rigide ("mal di denti implica carie" è falso): usa la probabilità come grado di credenza e, con la teoria dell'utilità, sceglie l'azione con utilità attesa massima. Strumenti: distribuzione congiunta completa (tabella di $2^n$ valori), marginalizzazione $P(Y\mid e)=\alpha\sum_h P(Y,e,h)$ con costante di normalizzazione $\alpha=1/P(e)$, regola del prodotto $P(a,b)=P(a\mid b)P(b)$, regola di Bayes $P(\text{causa}\mid\text{effetto})=\frac{P(\text{effetto}\mid\text{causa})P(\text{causa})}{P(\text{effetto})}$ e indipendenza condizionata $P(a,b\mid c)=P(a\mid c)P(b\mid c)$, che permette di scomporre la congiunta in tabelle più piccole. Un test con sensibilità 73% e falsi positivi 12% su una malattia con prior 1/700 dà una probabilità a posteriori di appena 0,9%.Incertezza - probabilità, regola di Bayes e indipendenza condizionata →).

(Wind non è osservata in nessuna riga: è una variabile nascosta; due valori di Temperature sono mancanti.)

Apprendimento dei parametri con dati completi

Se i dati sono sufficienti e completi e la topologia è nota, la CPT di ogni nodo si stima dalle frequenze. È lo stimatore di massima verosimiglianza (maximum likelihood, ML): fra tutte le CPT possibili, la frequenza dei dati è quella che rende i dati più probabili.

Formula (stima ML). Per un nodo $B$ con genitori $A,C$: $$P(B=b_i\mid A=a_j,C=c_k)\approx\frac{N_{i,j,k}}{N_{j,k}},$$ dove $N_{i,j,k}$ conta le righe con $B=b_i,A=a_j,C=c_k$ e $N_{j,k}$ quelle con $A=a_j,C=c_k$ (qualunque $B$).

Esempio (golf). $P(\text{Play}=P)=9/14=0{,}643$, $P(N)=5/14=0{,}357$. Per Outlook dato $N$ ci sono $5$ righe di classe $N$ (righe $1,2,6,8,14$): $3$ sunny, $0$ overcast, $2$ rainy, quindi $P(\text{sunny}\mid N)=3/5=0{,}6$, $P(\text{overcast}\mid N)=0$, $P(\text{rainy}\mid N)=2/5=0{,}4$.

Il vantaggio della scomposizione è anche statistico: con $25$ campioni, la congiunta di tre booleane ha $8$ combinazioni (circa $3$ campioni ciascuna), mentre $P(A\mid D)$, $P(B\mid D)$, $P(D)$ hanno circa $6$, $6$, $12$ campioni per numero.

Smoothing di Laplace

Se un evento non compare nei dati, la sua stima è $0$, che significa "impossibile": nel golf $P(\text{overcast}\mid N)=0$ azzererebbe, nel prodotto del classificatore, ogni riga con overcast. Lo smoothing di Laplace parte da una distribuzione uniforme e la aggiorna con i dati.

Formula (Laplace). Variabile discreta con $k$ valori, $N$ campioni, in cui $x_i$ compare $m$ volte: $$P(x_i)=\frac{1+m}{k+N}.$$ (Per una CPT: $N$ = righe con quella configurazione dei genitori.)

Esempi. Moneta lanciata $10$ volte con sempre testa: $k=2$, $N=10$, $m_{\text{testa}}=10$, $m_{\text{croce}}=0$: $P(\text{testa})=11/12$, $P(\text{croce})=1/12\ne0$. Nel golf: $k=3$, $N=5$: $P(\text{sunny}\mid N)=4/8=0{,}5$, $P(\text{overcast}\mid N)=1/8=0{,}125$, $P(\text{rainy}\mid N)=3/8=0{,}375$ (somma $1$). Più dati ci sono, meno la correzione pesa: con $m=N$ la stima $(1+N)/(2+N)\to1$.

Dati incompleti

Due situazioni:

• valori mancanti: alcune celle di una o più variabili sono vuote (le due $\text{Temp}$ del golf);
• variabili nascoste: una variabile del modello non ha alcun dato (la $\text{Wind}$ del golf).

Valori mancanti

Alternative: eliminare la riga; usare un valore speciale "sconosciuto"; sostituire con il valore più frequente (moda); stimare il valore dalle altre variabili della riga, che è di solito la scelta migliore.

Procedura (completamento). (1) Stimare i parametri iniziali con le righe complete. (2) Per ogni riga con valori mancanti: istanziare le variabili note; con l'inferenza probabilistica calcolare la distribuzione a posteriori delle variabili mancanti; assegnare il valore più probabile. (3) Aggiungere la riga completata e ristimare i parametri.

Golf. Riga $1$: $\mathbf P(\text{Temp}\mid\text{sunny},\text{high},N)$. Poiché in un naive Bayes ogni attributo è indipendente dagli altri dato $\text{Play}$, vale $\mathbf P(\text{Temp}\mid\text{sunny},\text{high},N)=\mathbf P(\text{Temp}\mid N)$. Le righe complete di classe $N$ sono $4$ (temperature high, low, medium, medium): $\langle P(\text{low}),P(\text{medium}),P(\text{high})\rangle=\langle0{,}25;\ 0{,}5;\ 0{,}25\rangle$, quindi $\text{medium}$. Riga $9$: $\mathbf P(\text{Temp}\mid\text{sunny},\text{normal},P)=\mathbf P(\text{Temp}\mid P)$, dalle $8$ righe complete di classe $P$ (low $2$, medium $4$, high $2$): $\langle0{,}25;\ 0{,}5;\ 0{,}25\rangle$, di nuovo $\text{medium}$. Le slide completano entrambe con $\text{medium}$.

Variabili nascoste: l'algoritmo EM

Per stimare le CPT di una variabile nascosta si usa l'Expectation-Maximization (EM), che ripete due fasi:

• passo E (expectation): si stimano i valori mancanti con i parametri correnti, per inferenza;
• passo M (maximization): si aggiornano i parametri con i dati completati.

Algoritmo (EM). (1) Stimare con ML le CPT dei nodi completi (quelli per cui valore e genitori sono nei dati). (2) Inizializzare i parametri ignoti con valori casuali. (3) Con i parametri correnti stimare le variabili nascoste dalle altre, per inferenza. (4) Usare le stime per completare la tabella. (5) Ristimare i parametri delle nascoste. (6) Ripetere $3$-$5$ finché i parametri non cambiano in modo significativo.

Golf, variabile $\text{Wind}$. Si parte da $P(\text{Wind}\mid\text{Play})=\begin{pmatrix}0{,}5&0{,}5\\0{,}5&0{,}5\end{pmatrix}$. Con questa CPT si stima Wind in ogni riga e si sceglie il valore più probabile; nelle slide la tabella riempita ha Wind $=\text{no},\text{no},\text{no},\text{no},\text{yes},\text{yes},\text{yes},\text{no},\text{no},\text{no},\text{yes},\text{yes},\text{yes},\text{yes}$ (righe $1$-$14$). Ristimando:

• righe $N$ ($1,2,6,8,14$): Wind $=\text{no},\text{no},\text{yes},\text{no},\text{yes}$: $P(\text{no}\mid N)=3/5=0{,}60$, $P(\text{yes}\mid N)=0{,}40$;
• righe $P$ ($3,4,5,7,9,10,11,12,13$): $\text{no},\text{no},\text{yes},\text{yes},\text{no},\text{no},\text{yes},\text{yes},\text{yes}$: $P(\text{no}\mid P)=4/9=0{,}44$, $P(\text{yes}\mid P)=5/9=0{,}56$.

Nuova CPT $\begin{pmatrix}0{,}60&0{,}40\\0{,}44&0{,}56\end{pmatrix}$ (righe $N$, $P$; colonne no, yes), come nelle slide. Osservazione: con una CPT uniforme tutte le righe hanno la stessa probabilità $0{,}5$ per i due valori di Wind, quindi la scelta del primo riempimento è arbitraria; inoltre, essendo Wind una foglia senza dati, la verosimiglianza dei dati osservati non dipende affatto da $P(\text{Wind}\mid\text{Play})$ e la CPT non è identificabile. L'EM è utile quando la variabile nascosta ha come figli variabili osservate, che portano informazione.

Apprendimento della struttura

Obiettivo: ottenere il grafo dai dati. È difficile perché:

• il numero di strutture possibili cresce più che esponenzialmente nel numero di variabili: i DAG su $n=1,2,3,4,5,6,7$ nodi sono $1,\ 3,\ 25,\ 543,\ 29281,\ 3781503,\ 1138779265$ (già per $n=10$ superano $4\cdot10^{18}$);
• servono molti dati per misure statistiche affidabili.

Esistono tecniche per alberi, poli-alberi e DAG generici. Per il caso generale due classi:

• metodi globali: ricerca euristica nello spazio delle strutture a partire da una iniziale, con una variazione a ogni passo; la struttura migliore si sceglie con un punteggio (score, ad esempio BIC, MDL);
• metodi locali: si valutano le relazioni di (in)dipendenza tra sottoinsiemi di variabili, e si ricostruisce la rete passo passo. Il più noto è l'algoritmo PC, che si ritrova nella scoperta causale (Controfattuali e scoperta causaleUn controfattuale chiede «cosa sarebbe successo a questo individuo se X fosse stato x?»: $Y_x(u)$ è il valore di Y nel modello $M_x$ (equazione di X sostituita da X=x) con le esogene fissate al valore $u$ dell'individuo. Si calcola in tre passi: abduzione (aggiornare $P(U)$ con l'evidenza, nel caso deterministico ricavare $u$ dai dati), azione (sostituire l'equazione di X con la costante), predizione (calcolare Y in $M_x$). Differenza con gli interventi: $P(Y\mid do(x))\equiv P(Y_x)$ è a livello di popolazione, ma in generale $P(Y\mid do(x),z)\ne P(Y_x\mid z)$ perché nel secondo il condizionamento è prima dell'intervento. Nei modelli lineari gaussiani i coefficienti di regressione di ogni variabile sui suoi genitori sono i coefficienti strutturali, e l'effetto causale totale è la somma dei prodotti dei coefficienti lungo i cammini non backdoor. Scoperta causale (algoritmo PC): si parte dal grafo completo non orientato, si toglie ogni arco X−Y con $X\perp!!!\perp Y\mid Z$ a $|Z|$ crescente (test di Fisher), si orientano i collider (Z mai nell'insieme separatore), poi si propaga l'orientamento; catena e fork non sono distinguibili (stessa classe di equivalenza di Markov). Assunzioni: aciclicità, sufficienza, Markov, fedeltà. Argomento del corso da 9 CFU.Controfattuali e scoperta causale →).

Funzioni di punteggio e BIC

Il punteggio misura quanto bene una struttura $m$ spiega i dati $D$ e deve bilanciare precisione e complessità: la sola verosimiglianza premia sempre il grafo con più archi (più parametri adattano meglio i dati, fino all'overfitting, vedi Overfitting, ridge regression e cross-validationUna buona prestazione sul training non basta: serve stimare quella su dati nuovi. La cross-validation (K-fold: $k$ parti, ciascuna a turno come test, errore medio; Monte Carlo: $k$ divisioni casuali con quota di test $q$; leave-one-out se $k=n$) evita di dipendere da una sola divisione casuale. L'errore atteso si scompone in $\text{bias}^2+\text{varianza}+\sigma^2$: i modelli semplici fanno underfitting (bias alto), quelli complessi overfitting (varianza alta). La regolarizzazione aggiunge alla perdita una penalità: la ridge regression minimizza $|y-X\beta|^2+\lambda\sum_{j\ge1}\beta_j^2$ e ha soluzione $\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ (l'intercetta non si penalizza, le feature si standardizzano): riduce i coefficienti, rende l'inversa stabile con feature collineari, e $\lambda$ è un iperparametro scelto con la validazione (cross-validation annidata per non contaminare il test).Overfitting, ridge regression e cross-validation →).

Formula (BIC). $$\text{BIC}(m)=\log P(D\mid m)-\frac{|m|\,\log|D|}{2},$$ con $|D|$ numero di campioni e $|m|$ numero di parametri indipendenti del modello. Il secondo termine penalizza la complessità. Si sceglie il massimo.

Esempio svolto. Due booleane $A,B$, $N=20$ campioni: $(a,b)\times7$, $(a,\neg b)\times2$, $(\neg a,b)\times3$, $(\neg a,\neg b)\times8$. Confrontare $m_0$ (nessun arco) e $m_1$ ($A\to B$), con $\log$ naturale.

• $m_0$: $9$ campioni con $a$, $11$ con $\neg a$; $10$ con $b$, $10$ con $\neg b$. $\log P(D\mid m_0)=9\ln\frac9{20}+11\ln\frac{11}{20}+10\ln\frac12+10\ln\frac12=-7{,}187-6{,}576-13{,}863=-27{,}626$. Parametri $|m_0|=2$ ($P(a)$ e $P(b)$): penalità $2\ln20/2=2{,}996$. $\text{BIC}=-30{,}621$.
• $m_1$: stesso termine per $A$ ($-13{,}763$); per $B$ con $a$: $7\ln\frac79+2\ln\frac29=-4{,}767$; con $\neg a$: $3\ln\frac3{11}+8\ln\frac8{11}=-6{,}446$; $\log P(D\mid m_1)=-24{,}976$. Parametri $|m_1|=3$ ($P(a),P(b\mid a),P(b\mid\neg a)$): penalità $3\ln20/2=4{,}494$. $\text{BIC}=-29{,}469$.

$-29{,}469>-30{,}621$: i dati giustificano l'arco $A\to B$ (la verosimiglianza guadagna $2{,}65$, la penalità costa $1{,}50$). Con dati più deboli (per esempio $5,4,5,6$ campioni nelle quattro celle) la verosimiglianza guadagna solo $0{,}10$ e vince la rete senza arco ($\text{BIC}=-30{,}621$ contro $-32{,}018$).

Ricerca: hill climbing

Algoritmo (hill climbing). (1) Generare una struttura iniziale (per esempio un albero). (2) Calcolarne il punteggio. (3) Aggiungere, togliere o invertire un arco. (4) Calcolare il punteggio della nuova struttura. (5) Se migliora, tenere la modifica; altrimenti tornare alla precedente. (6) Ripetere $3$-$5$ finché non ci sono più miglioramenti.

È una ricerca locale (Ricerca informata - greedy best-first, A* ed euristicheLa ricerca informata usa una funzione euristica h(n) che stima il costo da n all'obiettivo. Greedy best-first espande il nodo con h minima (f = h): veloce ma non ottima né in generale completa. A∗ usa f = g + h e, con h ammissibile (mai sopra il costo vero) e, sul grafo con tabella dei raggiunti, consistente, è completa e ottima; il fattore di diramazione effettivo b∗ misura la qualità di h, e un'euristica che ne domina un'altra espande meno nodi. A∗ pesato (f = g + W h) rinuncia all'ottimalità per espandere meno nodi. Per i problemi in cui conta solo lo stato finale esiste la ricerca locale (hill climbing, simulated annealing, beam search, algoritmi genetici). In appendice, cenni di ricerca avversaria (minimax e potatura alfa-beta).Ricerca informata - greedy best-first, A* ed euristiche →): può fermarsi in un massimo locale. Fissare un ordinamento causale delle variabili (i genitori possono essere solo variabili precedenti) riduce il numero di strutture da valutare.

Codice

import math
from collections import Counter

# dati del golf (slide): Outlook, Temp (None = mancante), Humidity, Play; Wind non osservata
dati = [("sunny", None, "high", "N"), ("sunny", "high", "high", "N"),
        ("overcast", "high", "high", "P"), ("rainy", "medium", "high", "P"),
        ("rainy", "low", "normal", "P"), ("rainy", "low", "normal", "N"),
        ("overcast", "low", "normal", "P"), ("sunny", "medium", "high", "N"),
        ("sunny", None, "normal", "P"), ("rainy", "medium", "normal", "P"),
        ("sunny", "medium", "normal", "P"), ("overcast", "medium", "high", "P"),
        ("overcast", "high", "normal", "P"), ("rainy", "medium", "high", "N")]
valori = {0: ["sunny", "overcast", "rainy"], 1: ["low", "medium", "high"], 2: ["high", "normal"]}

def cpt(col, classe, laplace=False):
    """P(attributo | Play=classe) dalle righe complete; con laplace: (1+m)/(k+N)."""
    righe = [r for r in dati if r[3] == classe and r[col] is not None]
    n, k = len(righe), len(valori[col])
    cont = Counter(r[col] for r in righe)
    return {v: ((1 if laplace else 0) + cont[v]) / ((k if laplace else 0) + n)
            for v in valori[col]}

print("ML       P(Outlook|N):", cpt(0, "N"))
print("Laplace  P(Outlook|N):", cpt(0, "N", laplace=True))
print("ML       P(Temp|N)   :", cpt(1, "N"))

def bic(dataset, genitori):
    """BIC = log P(D|m) - |m| log|D| / 2, per variabili binarie A, B (genitori di B: [] o [0])."""
    N, ll, par = len(dataset), 0.0, 0
    for var, gen in ((0, []), (1, genitori)):
        conteggi = Counter((tuple(r[g] for g in gen), r[var]) for r in dataset)
        for g in {c[0] for c in conteggi}:
            tot = conteggi[(g, 0)] + conteggi[(g, 1)]
            for x in (0, 1):
                if conteggi[(g, x)]:
                    ll += conteggi[(g, x)] * math.log(conteggi[(g, x)] / tot)
        par += 2 ** len(gen)               # un parametro per ogni riga della CPT
    return ll - par * math.log(N) / 2

D = [(1, 1)] * 7 + [(1, 0)] * 2 + [(0, 1)] * 3 + [(0, 0)] * 8   # (A, B), N = 20
print(bic(D, []), bic(D, [0]))              # senza arco, con A -> B

Output: {'sunny': 0.6, 'overcast': 0.0, 'rainy': 0.4}, {'sunny': 0.5, 'overcast': 0.125, 'rainy': 0.375}, {'low': 0.25, 'medium': 0.5, 'high': 0.25}, -30.62145... -29.46920....

Domande d'esame

1. Con $N=10$ lanci e $m_{\text{testa}}=10$, stimare $P(\text{croce})$ con smoothing di Laplace. $(1+0)/(2+10)=1/12\approx0{,}083$.
2. A cosa serve il termine $|m|\log|D|/2$ nel BIC? Penalizza il numero di parametri: evita di scegliere sempre la rete con più archi.
3. Descrivere i due passi dell'algoritmo EM. E: con i parametri correnti si stimano (per inferenza) i valori delle variabili nascoste; M: si ristimano i parametri con i dati così completati; si ripete fino a convergenza.
4. Perché si impiega lo smoothing nella stima delle CPT? Per evitare probabilità zero per eventi non osservati, che renderebbero impossibili casi solo rari.