Esercizio - Campionamento - rejection sampling e likelihood weighting sulla rete sprinkler (slide del corso)

Soluzione

I valori esatti si ottengono con l'enumerazione (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 →).

1. Rejection sampling

Stima. Si normalizza il vettore dei conteggi $\langle8,19\rangle$ sui $27$ campioni accettati: $$\hat{\mathbf P}(R\mid s)=\frac{\langle8,\,19\rangle}{27}=\langle0{,}296;\ 0{,}704\rangle.$$ Esatto. $P(r,s)=\sum_cP(c)P(s\mid c)P(r\mid c)=0{,}5\cdot0{,}1\cdot0{,}8+0{,}5\cdot0{,}5\cdot0{,}1=0{,}04+0{,}025=0{,}065$ e $P(s)=0{,}5\cdot0{,}1+0{,}5\cdot0{,}5=0{,}3$, quindi $P(r\mid s)=0{,}065/0{,}3=0{,}2167$. La stima ($0{,}296$) è troppo alta: con soli $27$ campioni utili lo scarto tipico è $\sqrt{0{,}217\cdot0{,}783/27}=0{,}079$, e $0{,}296$ dista circa uno scarto tipico dal valore esatto. Inoltre la frazione di campioni accettati attesa è $P(s)=0{,}3$ ($27$ su $100$: in linea). Con $10^5$ campioni la stima è $0{,}2127$–$0{,}2149$ a seconda del seme, con circa $30\,000$ campioni accettati.

2. Likelihood weighting

Peso di un campione. Evidenza $c,g$; si campionano $S$ e $R$. Peso $w=P(c)\cdot P(g\mid s,r)$ (il fattore $P(c)=0{,}5$ per $C$ e $P(g\mid\text{parents})$ per $G$):

Stima dai $10$ campioni. Le slide mostrano $6$ campioni $(\neg s,r)$ con $w=0{,}45$, $2$ campioni $(s,r)$ con $w=0{,}475$ e $2$ campioni $(\neg s,\neg r)$ con $w=0{,}05$: $$W_{r}=6\cdot0{,}45+2\cdot0{,}475=2{,}7+0{,}95=3{,}65,\qquad W_{\neg r}=2\cdot0{,}05=0{,}10,$$ $$\hat{\mathbf P}(R\mid c,g)=\frac{\langle3{,}65,\ 0{,}10\rangle}{3{,}75}=\langle0{,}9733;\ 0{,}0267\rangle.$$ Esatto (enumerazione, $S$ nascosta): $$P(r,c,g)=P(c)P(r\mid c)\sum_sP(s\mid c)P(g\mid s,r)=0{,}5\cdot0{,}8\cdot(0{,}1\cdot0{,}95+0{,}9\cdot0{,}90)=0{,}4\cdot0{,}905=0{,}362,$$ $$P(\neg r,c,g)=0{,}5\cdot0{,}2\cdot(0{,}1\cdot0{,}90+0{,}9\cdot0{,}10)=0{,}1\cdot0{,}18=0{,}018,\qquad\mathbf P(R\mid c,g)=\langle0{,}362;\,0{,}018\rangle/0{,}380=\langle0{,}9526;\ 0{,}0474\rangle.$$ La stima da $10$ campioni ($0{,}973$) è vicina al valore esatto ($0{,}953$). Con $N=10^5$ il codice sotto dà $0{,}9533$.

3. Gibbs sampling

Variabili libere: $C$ e $R$ (evidenza $s,g$ fissata). Stato iniziale $[c,s,\neg r,g]$.

1. Ricampionare $C$ con blanket $\{S,R\}$ (genitori nessuno, figli $S$ e $R$): $P(C\mid s,\neg r)\propto P(C)P(s\mid C)P(\neg r\mid C)$: per $c$ $0{,}5\cdot0{,}1\cdot0{,}2=0{,}01$, per $\neg c$ $0{,}5\cdot0{,}5\cdot0{,}9=0{,}225$; normalizzando, $P(c\mid s,\neg r)=0{,}01/0{,}235=0{,}0426$. Si estrae da $\langle0{,}0426;\ 0{,}9574\rangle$: nell'esempio esce $\neg c$.
2. Ricampionare $R$ con blanket $\{C,S,G\}$: $P(R\mid\neg c,s,g)\propto P(R\mid\neg c)P(g\mid s,R)$: per $r$ $0{,}1\cdot0{,}95=0{,}095$, per $\neg r$ $0{,}9\cdot0{,}90=0{,}81$; $P(r\mid\neg c,s,g)=0{,}095/0{,}905=0{,}105$. Nell'esempio esce $r$: stato $[\neg c,s,r,g]$.

Esatto. $P(r,s,g)=\sum_cP(c)P(s\mid c)P(r\mid c)P(g\mid s,r)=0{,}95\cdot0{,}065=0{,}06175$ e $P(\neg r,s,g)=0{,}90\cdot(0{,}5\cdot0{,}1\cdot0{,}2+0{,}5\cdot0{,}5\cdot0{,}9)=0{,}90\cdot0{,}235=0{,}2115$: $P(r\mid s,g)=0{,}06175/0{,}27325=0{,}2260$. La stima delle slide ($20$ volte $r$ su $80$ passi: $0{,}25$) è vicina. Con $10^5$ passi dopo $1000$ di burn-in: $0{,}2234$.

4. Markov blanket di $S$

$S$ ha genitore $C$, figlio $G$ e $G$ ha un altro genitore, $R$. Dunque il blanket è $\{C,R,G\}$ = $\{\text{Cloudy},\text{Rain},\text{GrassWet}\}$: sono corrette tutte e tre le prime opzioni (e quindi "nessuna delle precedenti" è sbagliata). Per $S$ vale $P(s\mid C,R,G)\propto P(s\mid C)P(g\mid s,R)$.

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.9533..., 0.22344; valori esatti $0{,}2167$, $0{,}9526$, $0{,}2260$. Con altri semi cambiano le ultime cifre: ad esempio la rejection con $10^5$ campioni dà $0{,}2149$ (seme $3$) e il likelihood weighting con $10$ campioni $0{,}878$ (seme $3$), con $100$ campioni $0{,}929$, con $1000$ $0{,}948$.