Esercizio - Filtraggio nell'ombrello con pomegranate (laboratorio 6)

Soluzione

Le tre domande a mano (filtraggio, $\mathbf f_t=\alpha\,\mathbf O_t\mathbf T^\top\mathbf f_{t-1}$). Con $\mathbf f_0=\langle0{,}5;\ 0{,}5\rangle$ la predizione del giorno $1$ è $\langle0{,}5;\ 0{,}5\rangle$.

1. $u_1$ e $u_2$. Giorno $1$: $\langle0{,}9\cdot0{,}5;\,0{,}2\cdot0{,}5\rangle=\langle0{,}45;\,0{,}10\rangle\to\mathbf f_1=\langle0{,}818;\ 0{,}182\rangle$. Predizione del giorno $2$: $\langle0{,}7\cdot0{,}818+0{,}3\cdot0{,}182;\ \dots\rangle=\langle0{,}627;\ 0{,}373\rangle$; con $u_2$: $\langle0{,}9\cdot0{,}627;\,0{,}2\cdot0{,}373\rangle=\langle0{,}5645;\,0{,}0745\rangle$, somma $0{,}639$: $$P(r_2\mid u_1,u_2)=0{,}8834.$$
2. $\neg u_1$ e $u_2$. Giorno $1$: $\langle0{,}1\cdot0{,}5;\,0{,}8\cdot0{,}5\rangle=\langle0{,}05;\,0{,}40\rangle\to\mathbf f_1=\langle0{,}111;\ 0{,}889\rangle$ (senza ombrello è improbabile che piova). Predizione del giorno $2$: $\langle0{,}7\cdot0{,}111+0{,}3\cdot0{,}889;\ \dots\rangle=\langle0{,}344;\ 0{,}656\rangle$; con $u_2$: $\langle0{,}9\cdot0{,}344;\,0{,}2\cdot0{,}656\rangle=\langle0{,}310;\,0{,}131\rangle$: $$P(r_2\mid\neg u_1,u_2)=0{,}310/0{,}441=0{,}7028.$$
3. Solo $u_2$. La predizione del giorno $2$ con il prior uniforme resta $\langle0{,}5;\ 0{,}5\rangle$ (la matrice $\mathbf T$ è simmetrica, quindi il prior uniforme è stazionario), poi con $u_2$: $\langle0{,}45;\,0{,}10\rangle\to$ $$P(r_2\mid u_2)=0{,}8182.$$

Confronto: $0{,}883$ (due ombrelli) $>0{,}818$ (un solo ombrello, nessuna informazione sull'altro giorno) $>0{,}703$ (ombrello mancante ieri, presente oggi). Il giorno $1$ senza ombrello abbassa la credenza su $r_1$ e, per la persistenza ($0{,}7$), anche su $r_2$.

Codice del corso (pomegranate 0.15)

from pomegranate import *

Rain0 = DiscreteDistribution({'y': 0.5, 'n': 0.5})          # prior su Rain0

Rain1 = ConditionalProbabilityTable(                         # transizione: [Rain_N, Rain_N+1, prob]
        [['y', 'y', 0.7], ['y', 'n', 0.3],
         ['n', 'y', 0.3], ['n', 'n', 0.7]], [Rain0])
Rain2 = ConditionalProbabilityTable(
        [['y', 'y', 0.7], ['y', 'n', 0.3],
         ['n', 'y', 0.3], ['n', 'n', 0.7]], [Rain1])
Umbrella1 = ConditionalProbabilityTable(                     # sensore: [Umbrella, Rain, prob]
        [['y', 'y', 0.9], ['y', 'n', 0.1],
         ['n', 'y', 0.2], ['n', 'n', 0.8]], [Rain1])
Umbrella2 = ConditionalProbabilityTable(
        [['y', 'y', 0.9], ['y', 'n', 0.1],
         ['n', 'y', 0.2], ['n', 'n', 0.8]], [Rain2])

In ogni riga di una ConditionalProbabilityTable le prime colonne sono i valori dei genitori, nell'ordine della lista finale (qui [Rain1]), poi il valore del nodo, poi la probabilità. Quindi ['y', 'n', 0.1] per Umbrella1 è $P(\text{Umbrella}=\text{n}\mid\text{Rain}=\text{y})=0{,}1$ e ['n', 'y', 0.2] è $P(\text{Umbrella}=\text{y}\mid\text{Rain}=\text{n})=0{,}2$. Il commento del corso ([ 'Umbrella', 'Rain', <probability>]) scrive l'ordine inverso, ma i valori sono giusti: danno $P(u\mid r)=0{,}9$ e $P(u\mid\neg r)=0{,}2$, come conferma l'output $0{,}8834$ uguale al calcolo a mano. Lo stesso vale per le transizioni: ['y', 'n', 0.3] per Rain1 è $P(\text{Rain}_1=\text{n}\mid\text{Rain}_0=\text{y})$.

Costruzione della rete, nodi (Rain0, Rain1, Umbrella1, Rain2, Umbrella2) con gli archi s1->s2, s2->s3, s2->s4, s4->s5; model.bake() ne fissa la topologia. Poi si aggiunge $\text{Rain}_3$ (s6, figlio di s4) e si richiama bake(). Le domande sono query con evidenza: una lista con un elemento per nodo (None = non osservato), nell'ordine [Rain0, Rain1, Umbrella1, Rain2, Umbrella2, Rain3].

scenario = [[None, None, 'y', None, 'y', None]]     # Q2.1: ombrello i giorni 1 e 2
results = model.predict_proba(scenario)
print(results[0][3].items())                         # indice 3 = Rain2

Output del corso:

• Q2.1 ($u_1,u_2$): (('y', 0.8833570412517777), ('n', 0.11664295874822232));
• Q2.2 ($\neg u_1$, scenario [None, None, 'n', None, 'y', None]): (('y', 0.7027707808564231), ('n', 0.2972292191435768));
• Q2.3 (nessuna informazione su $U_1$, [None, None, None, None, 'y', None]): (('y', 0.8181818181818179), ('n', 0.1818181818181821)).

Verifica con numpy

pomegranate non è installato: i tre valori sono stati ricalcolati con codice proprio. Si scrive la congiunta della rete srotolata $P(r_0)\prod_{k=1}^3P(r_k\mid r_{k-1})\prod_kP(u_k\mid r_k)$ come prodotto di CPT, si sommano le $2^4$ assegnazioni di $\text{Rain}_0,\dots,\text{Rain}_3$ che sono coerenti con l'evidenza e si normalizza:

import itertools
import numpy as np

Tm = {(1, 1): .7, (1, 0): .3, (0, 1): .3, (0, 0): .7}      # P(R_k | R_k-1), 1 = pioggia
Pu = {1: .9, 0: .2}                                           # P(u | R)

def post(profondita, evid):                                   # evid: {'U1': 1, ...}
    out = [np.zeros(2) for _ in range(profondita + 1)]
    for rs in itertools.product((1, 0), repeat=profondita + 1):
        p = .5                                                # prior su Rain0
        for a, b in zip(rs[:-1], rs[1:]): p *= Tm[(a, b)]
        for k in range(1, profondita + 1):
            if 'U%d' % k in evid:                             # evidenza sull'ombrello k
                p *= Pu[rs[k]] if evid['U%d' % k] else 1 - Pu[rs[k]]
        for k in range(profondita + 1): out[k][0 if rs[k] else 1] += p
    return [o / o.sum() for o in out]

print(post(3, {'U1': 1, 'U2': 1})[2])      # Q2.1
print(post(3, {'U1': 0, 'U2': 1})[2])      # Q2.2
print(post(3, {'U2': 1})[2])               # Q2.3

Output: [0.8834 0.1166], [0.7028 0.2972], [0.8182 0.1818], identici ai tre output di pomegranate.