Inferenza approssimata nelle reti bayesiane - campionamento

L'inferenza esatta (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 →) ha costo esponenziale nel caso peggiore. Con reti grandi si accetta una risposta approssimata: si estraggono $N$ campioni dalla distribuzione descritta dalla rete e si stima la probabilità cercata dalle frequenze. Per $N\to\infty$ la stima converge al valore vero (Legge dei grandi numeri e metodo Monte CarloSe X₁, X₂, ... sono i.i.d. con media μ, la media campionaria X̄ₙ = (X₁ + ... + Xₙ)/n converge a μ: in probabilità (legge debole, dimostrata con Chebyshev se la varianza è finita: P(|X̄ₙ − μ| > ε) ≤ σ²/(nε²)) e quasi certamente (legge forte). Metodo Monte Carlo: ∫ g = E[g(U)] si stima con la media di g(U₁), ..., g(Uₙ) per uniformi indipendenti.Legge dei grandi numeri e metodo Monte Carlo →): i metodi basati su campionamento casuale si chiamano algoritmi Monte Carlo. Due famiglie:

• campionamento diretto: prior sampling, rejection sampling, likelihood weighting;
• Markov chain Monte Carlo (MCMC): ogni campione deriva dal precedente con una piccola modifica casuale (Gibbs sampling).

Tutte richiedono di saper campionare da una distribuzione nota.

Campionare da una distribuzione discreta

Per una variabile booleana con $P(x)=p$: si estrae $u$ uniforme in $[0,1)$ (in Python random.random()) e si risponde "vero" se $u<p$, "falso" altrimenti; la probabilità di "vero" è esattamente $p$. Per più valori si ordinano le probabilità, si costruisce la funzione di distribuzione cumulata e si prende il primo valore la cui cumulata supera $u$ (inversione della cumulata). Con $20$ campioni da $\langle0{,}5;\ 0{,}2;\ 0{,}3\rangle$ non si ottengono esattamente quelle frequenze, ma ci si avvicina al crescere di $N$.

La rete di esempio: sprinkler

Il prato è bagnato (GrassWet, $G$) se piove (Rain, $R$) o se è acceso lo sprinkler (Sprinkler, $S$); entrambi dipendono dal cielo nuvoloso (Cloudy, $C$). Archi $C\to S$, $C\to R$, $S\to G$, $R\to G$.

Ordine topologico: $[C,S,R,G]$. Valori esatti che serviranno per il confronto: $P(s)=0{,}3$, $P(r)=0{,}45$, $P(c,r)=0{,}4$, $P(c,\neg s,r,g)=0{,}5\cdot0{,}9\cdot0{,}8\cdot0{,}9=0{,}324$, $P(r\mid s)=0{,}2167$, $P(r\mid c,g)=0{,}9526$, $P(r\mid s,g)=0{,}2260$.

Prior sampling

Algoritmo (prior-sample). Per ogni variabile nell'ordine topologico, campionare un valore da $\mathbf P(X_i\mid\text{parents}(X_i))$ con i valori già estratti dei genitori. Restituire l'evento completo.

Un campione. Da $P(C)$ esce $\text{true}$; da $P(S\mid c)$, $P(s\mid c)=0{,}1$ esce $\text{false}$; da $P(R\mid c)$ ($0{,}8$) esce $\text{true}$; da $P(G\mid\neg s,r)$ ($0{,}9$) esce $\text{true}$: evento $[c,\neg s,r,g]$. La probabilità che un campione sia proprio questo è quella della congiunta, $P(c)P(\neg s\mid c)P(r\mid c)P(g\mid\neg s,r)=0{,}324$: in $N$ campioni ci si aspetta il $32{,}4\%$ di eventi $[c,\neg s,r,g]$.

Stima. Con $N$ campioni, $\hat P(\text{evento})=N_{PS}(\text{evento})/N$, dove $N_{PS}$ conta i campioni in cui l'evento è vero.

Esempio (slide). Su $10$ campioni, $4$ hanno $C=\text{true}$ e $R=\text{true}$: $\hat P(c,r)=4/10=0{,}4$, uguale al valore esatto $0{,}4$. Ricampionando $1000$ volte e contando $[c,\neg s,r,g]$ si ottiene un numero vicino a $0{,}324$ (nel laboratorio $0{,}328$).

Precisione. Una frequenza su $N$ prove ha scarto tipico $\sqrt{p(1-p)/N}$, che per eventi rari è grande in valore relativo: per $p=0{,}000628$ (evento $[j,m,a,\neg b,\neg e]$ della rete dell'allarme) lo scarto relativo è $\sqrt{(1-p)/(pN)}\approx1{,}26$ per $N=1000$ ($126\%$) e $0{,}04$ per $N=10^6$. Ripetendo con $N=10,\,100,\,1000,\,10^4,\,10^5,\,10^6$ (seme $42$) si trovano $0,\ 0,\ 0{,}002,\ 0{,}001,\ 0{,}00072,\ 0{,}000644$ contro $0{,}000628$: servono molti campioni, come nell'esercizio Esercizio - Prior sampling sulla rete burglary (laboratorio 6).

Rejection sampling

Per una probabilità condizionata $\mathbf P(X\mid\mathbf e)=\alpha\,\mathbf P(X,\mathbf e)$ si conta solo nei campioni coerenti con l'evidenza.

Algoritmo (rejection sampling). Per $N$ volte: eseguire prior-sample; se il campione è coerente con $\mathbf e$, incrementare il contatore del valore di $X$ nel campione. Alla fine normalizzare il vettore dei contatori.

Esempio (slide). $P(\text{Rain}\mid s)$: $100$ campioni, $73$ con $S=\text{false}$ vengono scartati; dei $27$ con $s$, $8$ hanno $r$ e $19$ hanno $\neg r$. Stima $\langle8,19\rangle/27=\langle0{,}296;\ 0{,}704\rangle$ (valore esatto $P(r\mid s)=0{,}217$: con soli $27$ campioni utili l'errore è grande). Con $10^5$ campioni, circa il $30\%$ ($=P(s)$) è accettato e la stima è $0{,}215$.

Difetto. La frazione accettata è $P(\mathbf e)$: se l'evidenza è rara, quasi tutti i campioni si buttano. Nell'allarme, per $\mathbf P(B\mid j,m)$ solo $2109$ campioni su $10^6$ sono coerenti ($P(j,m)=0{,}0021$) e la stima di $P(b\mid j,m)$ è $0{,}275$ contro il valore esatto $0{,}284$.

Likelihood weighting

Idea: non generare campioni da scartare. Si fissano le variabili di evidenza ai valori osservati e si campionano solo le altre. Un campione che "contiene" l'evidenza non è però un campione della distribuzione condizionata: va pesato con la verosimiglianza dell'evidenza.

Algoritmo (likelihood weighting). Per ogni campione: $w\leftarrow1$; per ogni variabile $X_i$ nell'ordine topologico: se $X_i$ è di evidenza con valore $e_i$, $w\leftarrow w\cdot P(e_i\mid\text{parents}(X_i))$ (non si campiona); altrimenti si campiona da $\mathbf P(X_i\mid\text{parents}(X_i))$. Si somma $w$ nel contatore del valore di $X$. Alla fine si normalizza.

Esempio: $\mathbf P(\text{Rain}\mid c,g)$. Evidenza $C=\text{true}$, $G=\text{true}$; si campionano $S$ e $R$.

1. $w=1$; $C$ è evidenza: $w=1\cdot P(c)=0{,}5$.
2. $S$: campiono da $P(S\mid c)$, sia $\neg s$ (probabilità $0{,}9$).
3. $R$: campiono da $P(R\mid c)$, sia $r$ ($0{,}8$).
4. $G$ è evidenza: $w=0{,}5\cdot P(g\mid\neg s,r)=0{,}5\cdot0{,}9=0{,}45$.

Evento $[c,\neg s,r,g]$ con peso $0{,}45$. Se invece fosse uscito $s,r$, il peso sarebbe $0{,}5\cdot0{,}95=0{,}475$; con $\neg s,\neg r$, $0{,}5\cdot0{,}1=0{,}05$. Dieci campioni (tabella delle slide): sei con $[\neg s,r]$ (peso $0{,}45$), due con $[s,r]$ ($0{,}475$), due con $[\neg s,\neg r]$ ($0{,}05$). Somme dei pesi: $$W_{r}=6\cdot0{,}45+2\cdot0{,}475=3{,}65,\qquad W_{\neg r}=2\cdot0{,}05=0{,}1,\qquad \mathbf P(R\mid c,g)\approx\frac{\langle3{,}65;\ 0{,}1\rangle}{3{,}75}=\langle0{,}973;\ 0{,}027\rangle.$$ Esatto: $0{,}9526$ (enumerazione: $\mathbf P(R\mid c,g)=\alpha\langle0{,}8\cdot(0{,}1\cdot0{,}95+0{,}9\cdot0{,}9),\ 0{,}2\cdot(0{,}1\cdot0{,}9+0{,}9\cdot0{,}1)\rangle=\alpha\langle0{,}724;\,0{,}036\rangle=\langle0{,}9526;\ 0{,}0474\rangle$). Con $N=10^5$ si ottiene $0{,}9521$.

Vantaggi e limiti. Nessun campione scartato, quindi più efficiente del rejection. Ma con molte variabili di evidenza i pesi diventano piccoli e pochi campioni dominano la somma; inoltre l'evidenza posta "a valle" non influenza il campionamento delle variabili a monte, che sono estratte senza guardarla.

Markov blanket

Definizione (Markov blanket). Il Markov blanket di un nodo $X$ è l'insieme dei suoi genitori, dei suoi figli e degli altri genitori dei suoi figli. Dato il suo blanket, $X$ è condizionalmente indipendente da tutti gli altri nodi della rete.

Per lo sprinkler il blanket di $S$ è $\{C,R,G\}$ (genitore $C$, figlio $G$, altro genitore del figlio $R$); quello di $C$ è $\{S,R\}$; quello di $R$ è $\{C,S,G\}$. Dalle indipendenze: $$P(x\mid\text{mb}(X))\ \propto\ P(x\mid\text{parents}(X))\prod_{Y\in\text{children}(X)}P(y\mid\text{parents}(Y)).$$ Sono i soli fattori della congiunta in cui compare $X$.

Gibbs sampling

Algoritmo (Gibbs sampling). Si parte da uno stato arbitrario coerente con l'evidenza $\mathbf e$ (variabili di evidenza fissate). A ogni passo si sceglie a caso una variabile non di evidenza $Z_i$ e si estrae un nuovo valore da $\mathbf P(Z_i\mid\text{mb}(Z_i))$, tenendo fissate tutte le altre. Ogni stato visitato è un campione; la frequenza dei valori di $X$ negli stati stima $\mathbf P(X\mid\mathbf e)$.

Esempio: $\mathbf P(R\mid s,g)$. Libere: $C$ e $R$. Stato iniziale $[c,s,\neg r,g]$.

1. Scelgo $C$ (blanket $\{S,R\}$): $P(c\mid s,\neg r)\propto P(c)P(s\mid c)P(\neg r\mid c)=0{,}5\cdot0{,}1\cdot0{,}2=0{,}01$ e $P(\neg c\mid s,\neg r)\propto0{,}5\cdot0{,}5\cdot0{,}9=0{,}225$, quindi $P(c\mid s,\neg r)=0{,}0426$. Estraggo, sia $\neg c$: stato $[\neg c,s,\neg r,g]$.
2. Scelgo $R$ (blanket $\{C,S,G\}$): $P(r\mid\neg c,s,g)\propto P(r\mid\neg c)P(g\mid s,r)=0{,}1\cdot0{,}95=0{,}095$ e $P(\neg r\mid\neg c,s,g)\propto0{,}9\cdot P(g\mid s,\neg r)=0{,}9\cdot0{,}9=0{,}81$: $P(r\mid\neg c,s,g)=0{,}105$. Estraggo, sia $r$: stato $[\neg c,s,r,g]$.
3. Si continua scegliendo a caso $C$ o $R$.

Se in $80$ passi $R$ vale $\text{true}$ $20$ volte, la stima è $\langle0{,}25;\ 0{,}75\rangle$ (esatto $0{,}226$). Con $10^5$ passi (dopo $1000$ di scarto) si ottiene $0{,}228$.

Particolarità rispetto al campionamento diretto: i campioni sono correlati (ognuno differisce poco dal precedente); i primi stati dipendono dal punto di partenza e si scartano (burn-in); l'evidenza è sempre rispettata, e anche quella a valle influenza le variabili a monte. Per questo con pochi passi la stima è molto instabile (con $100$ passi e seme $3$ si è ottenuto $0{,}69$).

Confronto e codice

import random

# rete sprinkler: (genitori, CPT = P(X=vero | genitori))
bn = {
    'C': ([], {(): 0.5}),
    'S': (['C'], {(True,): 0.1, (False,): 0.5}),
    'R': (['C'], {(True,): 0.8, (False,): 0.1}),
    'G': (['S', 'R'], {(True, True): 0.95, (True, False): 0.90,
                       (False, True): 0.90, (False, False): 0.10}),
}
ordine = ['C', 'S', 'R', 'G']                  # ordine topologico

def p_vero(var, val):
    genitori, cpt = bn[var]
    return cpt[tuple(val[g] for g in genitori)]

def prior_sample(rng):
    val = {}
    for x in ordine:                           # dai genitori ai figli
        val[x] = rng.random() < p_vero(x, val)
    return val

def rejection(X, ev, N, rng):
    c = {True: 0, False: 0}
    for _ in range(N):
        s = prior_sample(rng)
        if all(s[k] == v for k, v in ev.items()):   # scarta se contraddice l'evidenza
            c[s[X]] += 1
    tot = c[True] + c[False]
    return c[True] / tot, tot

def likelihood_weighting(X, ev, N, rng):
    W = {True: 0.0, False: 0.0}
    for _ in range(N):
        val, w = dict(ev), 1.0
        for x in ordine:
            p = p_vero(x, val)
            if x in ev:
                w *= p if ev[x] else 1 - p     # evidenza: si pesa, non si campiona
            else:
                val[x] = rng.random() < p      # altrimenti si campiona
        W[val[X]] += w
    return W[True] / (W[True] + W[False])

def gibbs(X, ev, N, rng, burn=1000):
    val = {x: rng.random() < 0.5 for x in ordine}
    val.update(ev)
    libere = [x for x in ordine if x not in ev]
    figli = {x: [y for y in ordine if x in bn[y][0]] for x in ordine}
    def p_cond(x):                             # P(x=vero | Markov blanket)
        pesi = []
        for v in (True, False):
            val[x] = v
            p = p_vero(x, val) if v else 1 - p_vero(x, val)
            for y in figli[x]:
                py = p_vero(y, val)
                p *= py if val[y] else 1 - py
            pesi.append(p)
        return pesi[0] / (pesi[0] + pesi[1])
    cont = 0
    for i in range(N + burn):
        x = rng.choice(libere)                 # una variabile libera alla volta
        val[x] = rng.random() < p_cond(x)
        if i >= burn:
            cont += val[X]
    return cont / N

rng = random.Random(1)
print(rejection('R', {'S': True}, 100000, rng))
print(likelihood_weighting('R', {'C': True, 'G': True}, 100000, rng))
print(gibbs('R', {'S': True, 'G': True}, 100000, rng))

Output (seme $1$): (0.2127..., 30204), 0.9532..., 0.22344, da confrontare con i valori esatti $0{,}2167$, $0{,}9526$, $0{,}2260$. Esempi svolti in Esercizio - Campionamento - rejection sampling e likelihood weighting sulla rete sprinkler (slide del corso).

Domande d'esame

1. Perché il rejection sampling può essere inefficiente? La frazione di campioni accettati è $P(\mathbf e)$: con evidenza rara quasi tutti i campioni vengono scartati.
2. Calcolare il peso di un campione di likelihood weighting per $\mathbf P(R\mid c,g)$ con $s=\text{false}$, $r=\text{true}$. $w=P(c)\cdot P(g\mid\neg s,r)=0{,}5\cdot0{,}9=0{,}45$.
3. Qual è il Markov blanket di $S$ nello sprinkler? $\{C,R,G\}$: genitore $C$, figlio $G$, altro genitore del figlio $R$.
4. Su $100$ campioni di prior sampling, $73$ hanno $S=\text{false}$; dei restanti, $8$ hanno $R=\text{true}$. Stima di $P(r\mid s)$ con rejection sampling? $8/27\approx0{,}296$.