Esercizio - Pseudocodice del rejection sampling (esempio d'esame 1)

Teoria: 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 → (sezione sul rejection sampling). Esercizio collegato: Esercizio - Campionamento - rejection sampling e likelihood weighting sulla rete sprinkler (slide del corso). Domanda 8 del quiz da 9 CFU e domanda 5 di quello da 6 CFU: stesso pseudocodice, stesse dieci opzioni in ordine diverso.

Testo. Data una funzione Prior_Sample(BN) che restituisce un campione casuale dalla rete bayesiana BN, completare lo pseudocodice in modo che sia un algoritmo di rejection sampling.

function Rejection_Sampling(X, E, BN, N) returns P(X | E)
  X: variabile query    E: evidenza osservata    BN: rete bayesiana    N: numero di campioni
  C: vettore che conta ogni possibile valore di X
  P: vettore distribuzione di probabilità P(X | E)
  Blank 1
    Blank 2
    Blank 3
      Blank 4
    Blank 5
  Blank 6
  Blank 7
  return P

Ogni blank è un menu con le dieci opzioni: if P[k] > 1 then; C[j] <- C[j] * P[j]; for k = 1 to N do; endfor; C[j] <- C[j] + 1 (con $j$ valore attuale di $X$); X <- Prior_Sample(BN); if X is consistent with evidence E then; P <- Normalize(C); endif; P[j] <- P[j] + C[j]. Avvertenza: ogni for e if deve chiudersi con endfor e endif.

Soluzione

function Rejection_Sampling(X, E, BN, N) returns P(X | E)
  for k = 1 to N do                              # Blank 1
    X <- Prior_Sample(BN)                        # Blank 2
    if X is consistent with evidence E then      # Blank 3
      C[j] <- C[j] + 1   (j = valore di X)       # Blank 4
    endif                                        # Blank 5
  endfor                                         # Blank 6
  P <- Normalize(C)                              # Blank 7
  return P

Logica: si genera un campione completo dalla distribuzione a priori; se è coerente con l'evidenza lo si usa, altrimenti si scarta (non si tocca nessun contatore). Dei campioni conservati si contano i valori della variabile query; alla fine si normalizza (si divide per il totale dei campioni accettati) per ottenere $P(X\mid E)$. Gli ultimi tre blank sono forzati dalle parentesi: endif chiude l'if, endfor il for, e la normalizzazione va dopo il ciclo.

Perché le tre opzioni restanti sono sbagliate

• if P[k] > 1 then: una probabilità non supera mai $1$, quindi la condizione non è mai vera; inoltre $P$ non esiste ancora dentro il ciclo, perché si calcola solo alla fine.
• C[j] <- C[j] * P[j]: i contatori si incrementano di $1$, non si moltiplicano per probabilità (sarebbe il likelihood weighting, con pesi, non il rejection sampling).
• P[j] <- P[j] + C[j]: $P$ non si accumula dentro il ciclo; si ottiene una sola volta normalizzando $C$.

Controllo numerico

Rete $C\to S$, $C\to R$ con $P(C)=0{,}5$, $P(S\mid C)=0{,}1$, $P(S\mid\neg C)=0{,}5$, $P(R\mid C)=0{,}8$, $P(R\mid\neg C)=0{,}2$; si cerca $P(R\mid S=\text{vero})$. Esatto: $\dfrac{0{,}5\cdot0{,}1\cdot0{,}8+0{,}5\cdot0{,}5\cdot0{,}2}{0{,}5\cdot0{,}1+0{,}5\cdot0{,}5}=\dfrac{0{,}09}{0{,}30}=0{,}30$.

import random
rng = random.Random(0); cont = {True: 0, False: 0}
for _ in range(200000):                                        # for k = 1 to N
    c = rng.random() < .5; s = rng.random() < (.1 if c else .5); r = rng.random() < (.8 if c else .2)   # Prior_Sample
    if s:                                                      # coerente con l'evidenza S = vero
        cont[r] += 1                                           # C[j] <- C[j] + 1
tot = sum(cont.values()); print({k: v / tot for k, v in cont.items()})   # Normalize(C)

Output: {True: 0.3012, False: 0.6988}, vicino al valore esatto $0{,}30$ (si accettano circa 60 000 campioni su 200 000, cioè la frazione $P(S)=0{,}30$).