Indépendances conditionnelles et réseaux bayésiens

Partie obligatoire

Dans ce TME, l'objectif est d'apprendre des réseaux bayésiens à partir de bases de données. Hormis la base asia, un exemple jouet relativement petit qui vous permettra de mettre au point les différents algorithmes du TME, et car, les autres bases correspondront à des distributions de probabilité de tailles raisonnables :

Apprendre un réseau bayésien consiste à apprendre sa structure graphique ainsi que les paramètres de ses distributions de probabilité conditionnelles. Pour réaliser la deuxième tâche, il suffit d'estimer les paramètres de chaque distribution conditionnelle par maximum de vraisemblance, comme vous l'avez fait dans le TME 3. Ici, nous nous focaliserons donc plutôt sur l'apprentissage de structure. Celle-ci reflétant des indépendances conditionnelles entre variables aléatoires, vous devrez exploiter des tests d'indépendance du {$\chi^2$} afin d'obtenir des structures graphiques les moins denses possibles (en termes de nombres d'arcs). Ainsi, alarm représente une distribution jointe de plus de 10¹⁶ événements élémentaires mais, quand cette distribution est décomposée grâce au graphe ci-dessous (les noeuds représentant les variables aléatoires), elle peut être décrite (sans perte d'informations) à l'aide de seulement 752 paramètres. Comme nous l'avons vu en cours, cette représentation permet également d'effectuer très rapidement des calculs probabilistes.

1. Lecture des données

Dans le code ci-dessous, la fonction read_csv : string -> (string np.array, int np.2D-array, dico{string -> int} np.array) vous permettra de lire les données des bases sur lesquelles vous allez travailler, et de les organiser sous une forme adéquate. Par exemple, une base de données est un fichier de la forme:

 X_0,X_1,X_2,X_3
 haut,gauche,petit,bas
 bas,droite,grand,gauche
 bas,gauche,moyen,bas

Dans cette base, nous avons 4 variables aléatoires nommées X_0, X_1, X_2, et X_3, et 3 enregistrements qui représentent des instanciations (observées) de ces 4 variables. Ainsi, X_0 a pour valeurs haut, bas et bas, X_1 a pour valeurs gauche, droite, gauche, etc.

La fonction read_csv prend en argument le nom d'un fichier CSV contenant une base de données et renvoie un triplet composé de :

• 1 tableau numpy de strings contenant les noms des variables aléatoires. Par exemple, pour la base ci-dessus, ce tableau correspond à:

 n.array (['X_0', 'X_1', 'X_2', 'X_3'])

• un tableau numpy 2D contenant les données du fichier CSV encodées sous forme numérique (les valeurs des variables aléatoires sont transformées en nombres entiers): chaque ligne de ce tableau représente les intanciations d'une variable aléatoire et chaque colonne représente un enregistrement de la base de données, c'est-à-dire une instanciation/observation de toutes les variables aléatoires. Pour la base ci-dessus, nous obtiendrions le tableau ci-dessous (la signification des nombres est indiquée dans le dictionnaire précisé plus bas):

 np.array ( [ [0, 1, 1],    # instanciations de la variable X_0
              [0, 1, 0],    # instanciations de la variable X_1
              [0, 1, 2],    # instanciations de la variable X_2
              [0, 1, 0]] )  # instanciations de la variable X_3

• un tableau numpy de dictionnaires faisant la correspondance, pour chaque variable aléatoire, entre l'encodage numérique du tableau 2D ci-dessus et les données du fichier CSV (le 1er dictionnaire correspond à la variable de la 1ère colonne du CSV, le 2ème dictionnaire à celle de la 2ème colonne, etc.). Ainsi, le dictionnaire est égal à :

 np.array( [ {'haut': 0, 'bas': 1},                  # encodage variable X_0
             {'gauche': 0, 'droite': 1},             # encodage variable X_1
             {'petit': 0, 'grand': 1, 'moyen': 2 },  # encodage variable X_2
             {'bas': 0, 'gauche': 1} ] )             # encodage variable X_3

Téléchargez le fichier 2015_tme5_asia.csv et lisez-le à l'aide de la fonction read_csv ci-dessous: la dernière instruction, names, data, dico = read_csv ( "2015_tme5_asia.csv" ), vous permettra de récupérer, séparément, les trois champs du triplet renvoyé par la fonction read_csv.

2. Statistique du {$\chi^2$} conditionnel

Soit deux variables aléatoires {$X$} et {$Y$}. Appelons {$N_{xy}$}, {$N_x$} et {$N_y$}, respectivement, le nombre d'occurrences du couple {$(X=x,Y=y)$} et des singletons {$X=x$} et {$Y=y$} dans la base de données. Alors, comme indiqué dans le cours 5, la statistique du {$\chi^2$} de {$X$} et {$Y$} est égale à :

{$$\chi^2_{X,Y} = \sum_x\sum_y\frac{\left(N_{xy} - \frac{N_x \times N_y}{N}\right)^2}{\frac{N_x \times N_y}{N}}$$}

où {$N$} représente le nombre de lignes de la base de données. Cette formule permet de tester l'indépendance entre les deux variables {$X$} et {$Y$}. On peut aisément généraliser celle-ci pour tester des indépendances conditionnellement à un ensemble de variables {$\mathbf{Z}$}:

{$$\chi^2_{X,Y|\mathbf{Z}} = \sum_x\sum_y\sum_{\mathbf{z}}\frac{\left(N_{xy\mathbf{z}} - \frac{N_{x\mathbf{z}} \times N_{y\mathbf{z}}}{N_{\mathbf{z}}}\right)^2}{\frac{N_{x\mathbf{z}} \times N_{y\mathbf{z}}}{N_{\mathbf{z}}}}$$}

où {$N_{xy\mathbf{z}}$}, {$N_{x\mathbf{z}}$}, {$N_{y\mathbf{z}}$} et {$N_{\mathbf{z}}$} représentent, respectivement, le nombre d'occurrences du triplet {$(X=x,Y=y,\mathbf{Z} = \mathbf{z})$}, des couples {$(X=x,\mathbf{Z} = \mathbf{z})$} et {$(Y=y,\mathbf{Z} = \mathbf{z})$}, et du singleton {$\mathbf{Z} = \mathbf{z}$}. Ainsi, si {$\mathbf{Z}$} est un ensemble de 3 variables aléatoires {$(A,B,C)$}, les valeurs {$\mathbf{z}$} seront des triplets {$(a,b,c)$}.

Afin de vous aider à calculer ces {$\chi^2$}, vous pourrez utiliser la fonction create_contingency_table : int np.2D-array x dico{string -> int} np.array x int x int x int list -> (int, np.2D-array) np.array ci-dessous. Celle-ci prend en argument le tableau 2D numpy data et le tableau de dictionnaires dico retournés à la fin de la question 1, ainsi que l'index x d'une variable aléatoire (0 = 1ère variable aléatoire (celle de la 1ère ligne de data), 1 = 2ème variable, etc.), l'index y d'une autre variable et une liste z d'index d'autres variables aléatoires. Elle renvoie un tableau de couples ({$N_{\mathbf{z}}, T_{X,Y,\mathbf{z}})$}, pour tous les {$\mathbf{z} \in\mathbf{Z}$}, où:

• {$N_{\mathbf{z}}$} représente le nombre d'occurences de {$Z=z$} dans la base de données. Par exemple, si la base de données est la suivante :

 X_0,X_1,X_2,X_3
 haut,gauche,petit,bas
 bas,droite,grand,gauche
 bas,gauche,moyen,bas

 data = np.array ( [ [0, 1, 1],    # instanciations de la variable X_0
                     [0, 1, 0],    # instanciations de la variable X_1
                     [0, 1, 2],    # instanciations de la variable X_2
                     [0, 1, 0]] )  # instanciations de la variable X_3

 resultat = array([ (2, array([[ 1.,  0.,  0.],        # Z = 0 => N_{Z=0} = 2
                               [ 0.,  0.,  1.]])),
                    (1, array([[ 0.,  0.,  0.],        # Z = 1 => N_{Z=1} = 1
                               [ 0.,  1.,  0.]])) ])

• {$T_{X,Y,\mathbf{z}}$} est un tableau 2D contenant le nombre d'occurrences {$N_{xy\mathbf{z}}$} des couples {$(X=x,Y=y)$} lorsque {$\mathbf{Z}=\mathbf{z}$}. La première dimension de ce tableau (les lignes) correspondent aux différentes valeurs de {$X$} et la deuxième (les colonnes) à celles de {$Y$}. Ainsi, le tableau en haut à droite de resultat est obtenu de la manière suivante: ce tableau correspond à des occurrences de {$(X,Y)$} lorsque {$\mathbf{Z}=0$}. on commence donc par extraire de data le sous-tableau correspondant à la première et à la troisième colonne (les colonnes où X_3=0) et on ne retient que les lignes correspondant à X_0 et X_2 (cf. les paramètres 0 et 2 passés en arguments de create_contingency_table). On obtient donc le sous-tableau:

 np.array ( [ [0, 1],    # instanciations de la variable X_0
              [0, 2]] )  # instanciations de la variable X_2

En utilisant la structure retournée par la fonction create_contingency_table, écrivez une fonction sufficient_statistics: int np.2D-array x dico{string -> int} np.array x int x int x int list -> float qui prend les mêmes arguments que la fonction create_contingency_table et qui renvoie la valeur de {$\chi^2_{X,Y|\mathbf{Z}}$}. Vous pourrez tirer profit du fait que {$N_{x\mathbf{z}} = \sum_{y} N_{xy\mathbf{z}}$} et {$N_{y\mathbf{z}} = \sum_{x} N_{xy\mathbf{z}}$}, ce qui revient à faire des sommes sur chaque ligne ou chaque colonne des tableaux {$T_{X,Y,\mathbf{z}}$}. Attention : il peut arriver que certains {$N_{\mathbf{z}}$} soient égaux à 0. Dans ce cas, vous ne tiendrez pas compte des {$N_{xy\mathbf{z}}$}, {$N_{x\mathbf{z}}$} et {$N_{y\mathbf{z}}$} correspondants dans la formule de {$\chi^2_{X,Y|\mathbf{Z}}$} (car vous feriez des divisions par 0, ce qui est mal).

Vous pourrez vérifier la validité de vos calculs: en utilisant la base de données "2015_tme5_asia.csv", vous devriez obtenir les résultats suivants:

3. Statistique du {$\chi^2$} et degré de liberté

Modifiez votre fonction sufficient_statistics afin qu'elle ne renvoie plus seulement {$\chi^2_{X,Y|\mathbf{Z}}$} mais plutôt un couple ({$\chi^2_{X,Y|\mathbf{Z}}$},DoF), où DoF est le nombre de degrés de liberté de votre statistique. Celui-ci est égal à :

{$$(|X|-1) \times (|Y|-1) \times |\{\mathbf{z} : N_{\mathbf{z}} \neq 0\}| $$}

où {$|X|$} représente le nombre de valeurs possibles que peut prendre la variable {$X$}, autrement dit, c'est la taille de son dictionnaire. Le dernier terme de l'équation est simplement le nombre de {$N_{\mathbf{z}}$} différents de 0.

Vous pourrez vérifier la validité de vos calculs: en utilisant la base de données "2015_tme5_asia.csv", vous devriez obtenir les résultats suivants:

4. Test d'indépendance

En cours, nous avons vu que, pour un risque {$\alpha$} donné, si la statistique {$\chi^2_{X,Y|\mathbf{Z}}$} est inférieure au seuil critique {$c_{\alpha}$} de la loi du {$\chi^2$} à DoF degrés de liberté, alors {$X$} et {$Y$} sont considérés comme indépendants conditionnellement à {$\mathbf{Z}$} ({$X \perp\hspace{-1.7mm}\perp Y | \mathbf{Z}$}). On peut reformuler cette propriété de la manière suivante :

{$$\text{p-value}(\chi^2_{X,Y|\mathbf{Z}}) \geq \alpha \Longleftrightarrow X \perp\hspace{-1.7mm}\perp Y | \mathbf{Z}$$}

La p-value d'un nombre x est l'intégrale de la fonction de densité de la loi du {$\chi^2$} de x à {$+\infty$} (autrement dit, c'est la surface de la partie grisée sur votre table du {$\chi^2$} à partir de l'abscisse x. On a donc p-value{$(c_{\alpha}) = \alpha$}). En statistiques, on considère qu'elle n'a du sens que si les valeurs du tableau de contingence sont toutes supérieures ou égales à 5 (autrement dit, un test d'indépendance du {$\chi^2$} n'est "valide" que si toutes les valeurs du tableau de contingence sont supérieures ou égales à 5). En informatique, on allège souvent cette règle en considérant que le test est valide dès lors que la valeur moyenne des cases est supérieure ou égale à 5. Cet allègement permet de tester la validité du test sans réaliser celui-ci : si le nombre de lignes du CSV est supérieure ou égale à {$d_{min} = 5 \times |X| \times |Y| \times |\mathbf{Z}|$}, le test est considéré comme valide.

Ecrivez une fonction indep_score : int np.2D-array x dico{string -> int} np.array x int x int x int list -> (float,int) qui, étant donné les mêmes paramètres que ceux de la question précédente, vous renvoie un couple contenant la p-value correspondant à {$\chi^2_{X,Y|\mathbf{Z}}$} ainsi que le nombre de degrés de liberté de cette statistique. Vous testerez au préalable si len ( data[0] ), le nombre de lignes/enregistrements de votre CSV, est supérieur ou non à {$d_{min}$} ; si c'est inférieur, vous renverrez le couple (-1,1), qui représente une indépendance. Vous pourrez vous aider de la fonction scipy.stats.chi2.sf ( x, DoF ) qui renvoie la p-value (x) pour une loi à DoF degrés de liberté.

Vous pourrez vérifier la validité de vos calculs: en utilisant la base de données "2015_tme5_asia.csv", vous devriez obtenir les résultats suivants:

Partie optionnelle

5. Meilleur candidat pour être un parent

Ecrivez une fonction best_candidate : int np.2D-array x dico{string -> int} np.array x int x int list x float -> int list qui, étant donné les tableaux data et dico calculés à la question 1, l'index d'une variable aléatoire {$X$}, la liste d'index d'un ensemble de variables aléatoires {$\mathbf{Z}$} et un risque {$\alpha$}, détermine la variable {$Y$} (en fait, l'index de sa colonne dans le CSV), parmi toutes celles à gauche de {$X$} dans le fichier CSV, qui est la plus dépendante de {$X$} conditionnellement à {$\mathbf{Z}$}, autrement dit, celle qui a la plus petite p-value. Si cette p-value est supérieure à {$\alpha$}, cela veut dire que {$\chi^2_{X,Y|\mathbf{Z}}$} est inférieur à {$c_{\alpha}$} et donc que {$Y$} est jugée indépendante de {$X$} conditionnellement à {$\mathbf{Z}$}. Votre fonction renverra une liste vide si {$Y$} est indépendante de {$X$} conditionnellement à {$\mathbf{Z}$}, sinon elle renverra une liste contenant {$Y$}. Vous pourrez tester votre fonction avec {$\alpha$} = 0.05:

6. Création des parents d'un noeud

Ecrivez une fonction create_parents ( data, dico, x, alpha ) qui, étant donné une variable aléatoire x et un niveau de risque alpha, retourne la liste z de ses parents dans le réseau bayésien. L'algorithme est le suivant : partez de z = l'ensemble vide, puis tant que best_candidate ( x, z, alpha ) vous renvoie une liste non vide [y], rajoutez y à z. Lorsque vous sortirez de cette boucle, toutes les autres variables seront indépendantes de x conditionnellement à z.

L'algorithme qui consiste à appliquer, pour chaque noeud/variable aléatoire, votre fonction create_parents correspond, en grande partie, à l'article suivant :

Gregory F. Cooper and Edward Herskovits (1992) "A Bayesian method for the induction of probabilistic networks from data", Machine Learning, Vol. 9, n°4, pp. 309-347.

Vous pourrez tester la validité de votre fonction :

7. Apprentissage de la structure d'un réseau bayésien

Ecrivez une fonction learn_BN_structure ( data, dico, alpha ) qui renvoie un tableau contenant, pour chaque noeud, la liste de ses parents. Ainsi, si votre fonction vous renvoie le tableau ci-dessous,

 array( [ [], [], [], [1, 0], [1], [3, 2], [4, 5], [5] ] )

les noeud correspondant aux 2 premières colonnes du CSV n'ont pas de parents, le noeud de la 3ème colonne a pour parent celui de la 1ère colonne, etc.

Pour visualiser plus aisément votre structure, utilisez la fonction display_BN ci-dessous. Celle-ci prend en paramètres :

1. le tableau des noms des variables aléatoires déterminé à la question 1
2. la structure que vous avez calculée avec votre fonction learn_BN_structure
3. un nom que vous voulez donner à votre réseau
4. un style pour afficher les noeuds

7. Fin de l'apprentissage et calcul probabiliste

Comme précisé au début du TME, apprendre un réseau bayésien consiste à déterminer sa structure graphique et estimer ses paramètres. Vous avez réalisé la première partie. La deuxième, plus simple, peut se faire par maximum de vraisemblance pour chaque table de probabilité des noeuds conditionnellement à leurs parents, comme dans le TME 3. Utilisez la fonction learn_parameters ci-dessous pour effectuer cette tâche. Cette fonction prend en paramètres la structure graphique que vous avez apprise ainsi que le nom du fichier CSV que vous avez utilisé pour votre apprentissage. Elle renvoie un réseau bayésien à la aGrUM. Pour pouvoir utiliser aGrUM, reportez-vous à la question 7 du TME 2.

Vous pouvez maintenant réaliser des calculs probabilistes :

• affichage de la taille du réseau bayésien

• affichage de la table de probabilité conditionnelle d'un noeud du réseau déterminé par son nom (1ère ligne du CSV):

• calcul de la probabilité marginale d'un noeud : P('bronchitis?'):

• affichage graphique d'une distribution de probabilité marginale

• calcul d'une distribution marginale a posteriori : P(bronchitis? | smoking? = true, turberculosis? = false )

8. (Bonus) Autres bases de données

Vous pouvez appliquer vos algorithmes sur des bases un peu plus conséquentes qu'asia: