Esercizio - Prior sampling sulla rete burglary (laboratorio 6)

Soluzione

Valore esatto. Prodotto delle CPT lungo la rete (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 →): $$P(j,m,a,\neg b,\neg e)=P(j\mid a)P(m\mid a)P(a\mid\neg b,\neg e)P(\neg b)P(\neg e)=0{,}9\cdot0{,}7\cdot0{,}001\cdot0{,}999\cdot0{,}998=0{,}000628,$$ che le slide arrotondano a $0{,}00063$.

Algoritmo (Inferenza approssimata nelle reti bayesiane - campionamentoL'inferenza esatta (enumerazione, eliminazione di variabili) diventa proibitiva nelle reti grandi; i metodi Monte Carlo la approssimano estraendo $N$ campioni dalla rete e contando. Prior sampling: si campiona ogni variabile dopo i genitori, $P(\text{evento})\approx N_{PS}/N$. Rejection sampling: per $\mathbf P(X\mid\mathbf e)$ si scartano i campioni incoerenti con $\mathbf e$ (spreca molto se $P(\mathbf e)$ è piccola). Likelihood weighting: si fissano le variabili di evidenza e ogni campione pesa il prodotto delle $P(e_i\mid\text{parents})$; non spreca campioni. Gibbs sampling (MCMC): partendo da uno stato, si ricampiona una variabile libera alla volta da $P(Z\mid\text{Markov blanket})$, dove il blanket è genitori, figli e altri genitori dei figli. Tutte convergono al valore vero per $N\to\infty$ con errore dell'ordine di $1/\sqrt N$.Inferenza approssimata nelle reti bayesiane - campionamento →). Per ogni campione si percorrono le variabili in ordine topologico ($B,E,A,J,M$) e si estrae ciascuna dalla sua CPT con i valori dei genitori già estratti. La probabilità dell'evento si stima con la frequenza: $\hat P=N_{PS}(j,m,a,\neg b,\neg e)/N$.

La funzione base samplegen estrae un booleano: restituisce vero se random.random() (uniforme in $[0,1)$) è minore di $P(X=\text{vero}\mid\text{genitori})$. Per leggere la riga giusta della CPT si usano i valori dei genitori come indici dell'array: Pdist[t][genitori], dove il primo indice è il valore della variabile ($t=0$ vero) e gli altri sono i valori dei genitori. Esempio: P_A_BE[t][f][t] è $P(a\mid\neg b,e)=0{,}29$.

Dettaglio utile del codice. Gli eventi sono salvati come liste di 't'/'f' nell'ordine ['B','E','A','J','M']: l'evento cercato è ['f','f','t','t','t'] (non $b$, non $e$, $a$, $j$, $m$). Attenzione all'inversione: nel codice t = 0 e f = 1, quindi 'f' if val[x] else 't' converte l'indice $1$ in 'f'.

Codice

import numpy as np
import random as rnd

t, f = 0, 1                                      # indici di 'vero' e 'falso'

def samplegen(Pdist, Parents=[]):
    """Campiona vero/falso da una CPT booleana; Pdist[t][genitori] = P(X = vero | genitori)."""
    assert len(Parents) < len(Pdist.shape)
    if rnd.random() < Pdist[t][tuple(Parents)]:
        return t
    return f

def parents(X):
    return [val[i] for i in par[X]]             # valori correnti dei genitori di X

def priorSampling(n_events):
    event = []
    for n in range(n_events):
        for x in var:                           # ordine topologico: B, E, A, J, M
            val[x] = samplegen(prd[x], parents(x))
        event.append(['f' if val[x] else 't' for x in var])
    return event

# rete burglary
P_B = np.array([0.001, 0.999])
P_E = np.array([0.002, 0.998])
P_A_BE = np.array([[[0.95, 0.94], [0.29, 0.001]], [[0.05, 0.06], [0.71, 0.999]]])   # [a][b][e]
P_J_A = np.array([[0.9, 0.05], [0.1, 0.95]])                                        # [j][a]
P_M_A = np.array([[0.7, 0.01], [0.3, 0.99]])                                        # [m][a]

var = ['B', 'E', 'A', 'J', 'M']
prd = {'B': P_B, 'E': P_E, 'A': P_A_BE, 'J': P_J_A, 'M': P_M_A}
par = {'B': [], 'E': [], 'A': ['B', 'E'], 'J': ['A'], 'M': ['A']}
val = {'B': f, 'E': f, 'A': f, 'J': f, 'M': f}

rnd.seed(42)
for N in [10, 100, 1000, 10000, 100000]:
    data = priorSampling(N)
    stima = data.count(['f', 'f', 't', 't', 't']) / len(data)    # B=f, E=f, A=t, J=t, M=t
    print(N, stima)

Output (seme $42$): 10 0.0, 100 0.0, 1000 0.002, 10000 0.0008, 100000 0.00071. La soluzione del corso (senza seme) dà $0$, $0$, $0{,}001$, $0{,}0007$, $0{,}0006$ e, con $10^6$ campioni, $0{,}000662$.

Cosa si osserva

L'evento ha probabilità $p=0{,}000628$: in $N$ campioni se ne aspettano $Np$.

• Con $N=10$ e $N=100$ ci si aspetta meno di un evento: la stima è quasi sempre $0$, che è troppo lontano per dire se l'algoritmo funziona.
• Con $N=1000$ basta un evento in più o in meno per cambiare la stima di $0{,}001$: una stima di $0{,}002$ è un errore relativo del $218\%$.
• Con $N\ge10^4$ la stima si avvicina a $0{,}00063$. Lo scarto tipico scende come $1/\sqrt N$: per dimezzarlo servono $4$ volte i campioni.
• Verifica: ripetendo $100$ volte l'esperimento (ogni stima è $\text{Binomiale}(N,p)/N$) la media delle stime coincide con il valore esatto e lo scarto tipico osservato coincide con quello atteso: a $N=10^5$ media $0{,}000629$, scarto $7{,}5\cdot10^{-5}$ (atteso $7{,}9\cdot10^{-5}$); a $N=10^6$ media $0{,}000628$, scarto $2{,}7\cdot10^{-5}$ (atteso $2{,}5\cdot10^{-5}$). Dunque la stima è corretta, solo rumorosa. Un campionatore indipendente scritto con vettori numpy dà $0{,}000641$ con $10^6$ campioni.

Conclusione: il prior sampling approssima correttamente $P(j,m,a,\neg b,\neg e)$, ma per eventi rari servono molti campioni; per probabilità condizionate con evidenza rara conviene il likelihood weighting (Esercizio - Campionamento - rejection sampling e likelihood weighting sulla rete sprinkler (slide del corso)).