Rappels de Probas / stats

Partie obligatoire

La planche de Galton

1. Loi de Bernoulli

Ecrivez une fonction bernoulli : float -> int qui prend en argument un paramètre réel p compris entre 0 et 1, et qui renvoie soit 1 (succès) avec la probabilité p, soit 0 (échec) avec la probabilité 1-p.

2. Loi binomiale

Ecrivez une fonction binomiale : int x float -> int qui prend en argument un paramètre entier n et un paramètre réel p et qui renvoie un nombre tiré aléatoirement selon la distribution binomiale {$\cal{B}(n,p)$}.

3. Histogramme de la loi binomiale

Dans cette question, on considère une planche de Galton de hauteur n. On rappelle que des bâtons horizontaux (oranges) sont cloués à cette planche comme le montre la figure ci-dessous. Des billes bleues tombent du haut de la planche et, à chaque niveau, se retrouvent à la verticale d'un des bâtons. Elles vont alors tomber soit à gauche, soit à droite du bâton, jusqu'à atteindre le bas de la planche. Ce dernier est constitué de petites boites dont les bords sont symbolisés par les lignes verticales grises. Chaque boite renferme des billes qui sont passées exactement le même nombre de fois à droite des bâtons oranges. Par exemple, la boite la plus à gauche renferme les billes qui ne sont jamais passées à droite d'un bâton, celle juste à sa droite renferme les billes passées une seule fois à droite d'un bâton et toutes les autres fois à gauche, et ainsi de suite.

La répartition des billes dans les boites suit donc une loi binomiale {$\cal{B}(n,0.5)$}. Ecrivez un script qui crée un array numpy de 1000 cases dont le contenu correspond à 1000 instanciations de la loi binomiale {$\cal{B}(n,0.5)$}. Afin de voir la répartition des billes dans la planche de Galton, tracez l'histogramme de ce tableau. Vous pourrez utiliser la fonction hist de matplotlib.pyplot:

Pour le nombre de bins, calculez le nombre de valeurs différentes dans votre tableau.

Visualisation d'indépendances

1. Loi normale centrée

On souhaite visualiser la fonction de densité de la loi normale. Pour cela, on va créer un ensemble de k points (x_i,y_i), pour des x_i équi-espacés variant de -2{$\sigma$} à 2{$\sigma$}, les y_i correspondant à la valeur de la fonction de densité de la loi normale centrée de variance {$\sigma^2$}, autrement dit {$\cal{N}(0,\sigma^2)$}.

Ecrivez une fonction normale : int x float -> float np.array qui, étant donné un paramètre entier k impair et un paramètre réel sigma renvoie l'array numpy des k valeurs y_i. Afin que l'array numpy soit bien symmétrique, on lèvera une exception si k est pair:

Vous vérifierez la validité de votre fonction en affichant grâce à la fonction plot les points générés dans une figure.

2. Distribution de probabilité affine

Dans cette question, on considère une généralisation de la distribution uniforme: une distribution affine, c'est-à-dire que la fonction de densité est une droite, mais pas forcément horizontale, comme le montre la figure ci-dessous:

Ecrivez une fonction proba_affine : int x float -> float np.array qui, comme dans la question précédente, va générer un ensemble de {$k$} points {$y_i$}, {$i = 0,...,k-1$}, représentant cette distribution (paramétrée par sa pente slope). On supposera que l'entier {$k$} est impair. Si la pente est égale à 0, c'est-à-dire si la distribution est uniforme, chaque point {$y_i$} devrait être égal à {$1/k$} (afin que la somme des {$y_i$} soit égale à 1). Si la pente est différente de 0, il suffit de choisir, pour tout {$i = 0,...,k-1,$}

{$$y_i = \frac{1}{k} + \left( i - \frac{k - 1}{2} \right) \times slope$$}

Vous pourrez aisément vérifier que la somme des {$y_i$} est alors éagle à 1. Afin que la distribution soit toujours positive (c'est quand même un minimum pour une distribution de probabilité), il faut que la pente slope ne soit ni trop grande ni trop petite. Le bout de code ci-dessous lèvera une exception si la pente est trop élevée et indiquera la pente maximale possible.

3. Distribution jointe

Ecrivez une fonction Pxy : float np.array x float np.array -> float np.2D-array qui, étant donné deux tableaux numpy de nombres réels à 1 dimension générés par les fonctions des questions précédentes et représentant deux distributions de probabilités {$P(A)$} et {$P(B)$}, renvoie la distribution jointe {$P(A,B)$} sous forme d'un tableau numpy à 2 dimensions de nombres réels, en supposant que {$A$} et {$B$} sont des variables aléatoires indépendantes. Par exemple, si:

alors Pxy ( PA,PB ) renverra le tableau :

4. Affichage de la distribution jointe

Le code ci-dessous permet d'afficher en 3D une probabilité jointe générée par la fonction précédente. Exécutez-le avec une probabilité jointe résultant de la combinaison d'une loi normale et d'une distribution affine. Si le contenu de la fenêtre est vide, redimensionnez celle-ci et le contenu devrait apparaître. Cliquez à la souris à l'intérieur de la fenêtre et bougez la souris en gardant le bouton appuyé afin de faire pivoter la courbe. Observez sous différents angles cette courbe. Refaites l'expérience avec une probaiblité jointe résultant de deux lois normales. Essayez de comprendre ce que signifie, visuellement, l'indépendance probabiliste. Vous pouvez également recommencer l'expérience avec le logarithme des lois jointes.

Indépendances conditionnelles

Dans cet exercice, on considère quatre variables aléatoires booléennes {$X$}, {$Y$}, {$Z$} et {$T$} ainsi que leur distribution jointe {$P(X,Y,Z,T)$} encodée en python de la manière suivante :

Ainsi, pour tout x,y,z,t {$\in \{0,1\}$}, P_XYZT[x][y][z][t] correspond à {$P(X=x,Y=y,Z=z,T=t)$} ou, en version abrégée, à {$P(x,y,z,t)$}.

1. Indépendance de {$X$} et {$T$} conditionnellement à {$(Y,Z)$}

On souhaite tester si les variables aléatoires {$X$} et {$T$} sont indépendantes conditionnellement à {$(Y,Z)$}. Pour cela, calculez à partir de P_XYZT le tableau P_YZ représentant la distribution {$P(Y,Z)$}. On rappelle que {$P(Y,Z) = \sum_{X,T} P(X,Y,Z,T)$}. Le tableau P_YZ est donc un tableau à deux dimensions, dont la première correspond à {$Y$} et la deuxième à {$Z$}. Si vous ne vous êtes pas trompé(e), vous devez obtenir le tableau suivant :

Ainsi, P_YZ[0][1] =0.084 = {$P(Y=0, Z=1)$}.

Ensuite, calculez le tableau P_XTcondYZ représentant la distribution {$P(X,T|Y,Z)$}. Ce tableau a donc 4 dimensions, chacune correspondant à une des variables aléatoires. De plus, les valeurs de P_XTcondYZ sont obtenues en utilisant la formule des probabilités conditionnelles: {$P(X,T|Y,Z) = P(X,Y,Z,T) / P(Y,Z)$}.

Calculez à partir de P_XTcondYZ les tableaux à 3 dimensions P_XcondYZ et P_TcondYZ représentant respectivement les distributions {$P(X|Y,Z)$} et {$P(T|Y,Z)$}. On rappelle que {$P(X|Y,Z) = \sum_T P(X,T|Y,Z)$}.

Enfin, testez si {$X$} et {$T$} sont indépendantes conditionnellement à {$(Y,Z)$}: si c'est bien le cas, on doit avoir {$P(X,T|Y,Z) = P(X|Y,Z) \times P(T|Y,Z)$}.

2. Indépendance de {$X$} et {$(Y,Z)$}

On souhaite maintenant déterminer si {$X$} et {$(Y,Z)$} sont indépendantes. Pour cela, commencez par calculer à partir de P_XYZT le tableau P_XYZ représentant la distribution {$P(X,Y,Z)$}. Ensuite, on calcule à partir de P_XYZ les tableaux P_X et P_YZ représentant respectivement les distributions {$P(X)$} et {$P(Y,Z)$}. On rappelle que {$P(X) = \sum_Y \sum_Z P(X,Y,Z)$}. Si vous ne vous êtes pas trompé(e), P_X doit être égal au tableau suivant :

Enfin, si {$X$} et {$(Y,Z)$} sont bien indépendantes, on doit avoir {$P(X,Y,Z) = P(X) \times P(Y,Z)$}.

Indépendances conditionnelles et consommation mémoire

Le but de cet exercice est d'exploiter les probabilités conditionnelles et les indépendances conditionnelles afin de décomposer une probabilité jointe en un produit de "petites probabilités conditionnelles". Cela permet de stocker des probabilités jointes de grandes tailles sur des ordinateurs "standards". Au cours de l'exercice, vous allez donc partir d'une probabilité jointe et, progressivement, construire un programme qui identifie ces indépendances conditionnelles.

Pour simplifier, dans la suite de cet exercice, nous allons considérer un ensemble {$X_0,\ldots,X_n$} de variables aléatoires binaires (elles ne peuvent prendre que 2 valeurs : 0 et 1).

1. Package de manipulation de probabilités

Comme vous avez pu le constater dans l'exercice précédent, manipuler des probabilités (conditionnelles) multivariées peut s'avérer complexe du fait que l'on doit utiliser des boucles for dépendant explicitement des variables aléatoires, ce qui rend difficile l'écriture d'un programme générique s'appliquant à toute distribution de probabilité. Pour pallier ce problème, dans cet exercice, vous allez exploiter pyAgrum, qui est un module python3 spécifiquement conçu pour manipuler des distributions de probabilité multivariées. Pour l'importer sur les machines de la PPTI, il suffit d'effectuer les instructions suivantes :

Chez vous, vous pouvez l'installer facilement (cf. http://agrum.gitlab.io/pages/pyagrum.html). Après importation, vous pouvez créer des variables aléatoires et des tables de probabilité de la manière suivante :

Les Potential de pyAgrum sont donc des tables de probabilité faciles à manipuler. Notamment, ils possèdent les opérateurs et méthodes suivants:

• Les opérateurs standards *, /, + et - peuvent être appliqués sur des Potential. Ainsi, si P_AB et P_B sont des Potential représentant respectivement les distributions {$P(A,B)$} et {$P(B)$}, alors P_AB / P_B est un Potential représentant {$P(A|B)$}.
• La méthode margSumOut prend en argument une liste de noms de variables (string list) et supprime par sommation toutes ces variables de la distribution de probabilité. Par exemple, si P_ABC est un Potential représentant {$P(A,B,C)$}, alors P_ABC.margSumOut(['B','C']) est un Potential représentant {$P(A) = \sum_B \sum_C P(A,B,C)$}
• La méthode margSumIn est similaire à la méthode précédente excepté que l'on supprime toutes les variables sauf celles en arguments. Ainsi P_ABC.margSumIn(['B','C']) est un Potential représentant {$P(B,C) = \sum_A P(A,B,C)$}.
• La méthode variablesSequence() renvoie la liste de toutes les LabelizedVariable sur lesquelles porte la distribution de probabilité. Ainsi, P_ABC.variablesSequence() renverra une liste contenant les LabelizedVariable correspondant aux variables {$A$}, {$B$} et {$C$}.

2. Téléchargement d'une distribution jointe

Téléchargez le fichier asia.txt et sauvegardez le sous le nom de asia.txt. Celui-ci contient une probabilité jointe sur un ensemble de 8 variables aléatoires. Le fichier est produit à partir du site web suivant :

http://www.bnlearn.com/bnrepository/

Le code suivant vous permettra de lire ce fichier et d'en récupérer un tableau numpy des LabelizedVariable ainsi que la probabilité jointe (sous forme d'un Potential) qu'il contient :

3. Test d'indépendance conditionnelle

Ecrivez une fonction conditional_indep : Potential x LabelizedVariable x LabelizedVariable x LabelizedVariable List x float -> bool qui prend en arguments une distribution de probabilité jointe {$P(X_1,\ldots,X_n)$}, deux variables distinctes {$X$} et {$Y$} ainsi qu'une liste de variables {$\mathbf{Z}$} dont on supposera que {$X$}, {$Y$} et {$\mathbf{Z}$} sont disjoints et que leur union est incluse dans {$\{X_1,\ldots,X_n\}$}, et enfin un nombre réel {$\epsilon$}. La fonction renverra un booléen indiquant si {$X$} est indépendant de {$Y$} conditionnellement à {$\mathbf{Z}$}. La démarche la plus simple pour déterminer cela est sans doute la suivante :

1. Calculez à partir de la probabilité jointe {$P(X_1,\ldots,X_n)$} la distribution {$P(X,Y,\mathbf{Z})$} en supprimant par sommation toutes les variables de {$P$} autres que {$X$}, {$Y$} et {$\mathbf{Z}$}. Dans un Potential, la liste des variables peut être obtenue grâce à la méthode variablesSequence(). Notez également que le nom d'une LabelizedVariable peut être obtenu en lui appliquant la méthode name().
2. Calculez à partir de {$P(X,Y,\mathbf{Z})$} les distributions {$P(X|\mathbf{Z})$} et {$P(Y|\mathbf{Z})$}.
3. Calculez {$Q(X,Y|\mathbf{Z}) = P(X,Y|\mathbf{Z}) - P(X|\mathbf{Z}) \times P(Y|\mathbf{Z})$}.
4. Soit Q le Potential représentant la fonction {$Q(X,Y|\mathbf{Z})$}. S'il y a indépendance entre {$X$} et {$Y$} conditionnellement à {$\mathbf{Z}$}, Q devrait être une table ne contenant que des 0. Malheureusement, du fait de l'imprécision des calculs en virgule flotante, il se peut que Q contienne des valeurs légèrement différentes de 0. Aussi, pour tester s'il y a bien indépendance, il est plus judicieux de vérifier si les valeurs absolues des éléments de la table Q (ou, de manière équivalente, leur max) sont inférieures à l'{$\epsilon$} passé en argument de votre fonction conditional_indep. Vous pouvez tester cela simplement en appelant Q.abs().max().

4. Compactage de probabilités conditionnelles

On sait que si un ensemble de variables aléatoires {${\cal S} = \{X_{i_0},\ldots,X_{i_{n-1}}\}$} peut être partitionné en deux sous-ensembles {$\cal K$} et {$\cal L$} (c'est-à-dire tels que {${\cal K} \cap {\cal L} = \emptyset$} et {${\cal K} \cup {\cal L} = \{X_{i_0},\ldots,X_{i_{n-1}}\}$}) tels qu'une variable {$X_{i_n}$} est indépendante de {${\cal L}$} conditionnellement à {${\cal K}$}, alors:

{$P(X_{i_n}|X_{i_0},\ldots,X_{i_{n-1}}) = P(X_{i_n} | {\cal K},{\cal L}) = P(X_{i_n} | {\cal K}).$}

C'est ce que nous avons vu au cours n°2 (cf. définition des probabilités conditionnelles). Cette formule est intéressante car elle permet de réduire la taille mémoire consommée pour stocker {$P(X_{i_n}|X_{i_0},\ldots,X_{i_{n-1}})$}: il suffit en effet de stocker uniquement {$P(X_{i_n} | {\cal K})$} pour obtenir la même information.

Ecrivez une fonction compact_conditional_proba : Potential x LabelizedVariable x float -> bool qui, étant donné une probabilité jointe {$P(X_{i_0},\ldots,X_{i_n})$}, une variable aléatoire {$X_{i_n}$} et un nombre réel {$\epsilon$}, retourne cette probabilité conditionnelle {$P(X_{i_n} | {\cal K})$}: pour cela, nous vous proposons l'algorithme itératif suivant:

{${\cal K} = {\cal S}$}
Pour tout {$X_{i_j} \in {\cal K}$} faire:

Fait
Retourner {$P(X_{i_n} | {\cal K})$}

Nous vous suggérons d'exploiter votre fonction conditional_indep afin de tester les indépendances conditionnelles.

Si vous souhaitez observer les valeurs des distributions de probabilité conditionnelles que vous avez créées, vous pouvez utiliser le code suivant, où proba est un Potential.

Afin que l'affichage soit plus facile à comprendre, il pourrait être judicieux de placer la variable {$X_{i_n}$} en premier dans la liste des variables du Potential, ce que l'on peut faire avec le code suivant :

5. Création d'un réseau bayésien

Un réseau bayésien est simplement la décomposition d'une distribution de probabilité jointe en un produit de probabilités conditionnelles: vous avez vu en cours que {$P(A,B) = P(A|B)P(B)$}, et ce quel que soient les ensembles de variables aléatoires disjoints {$A$} et {$B$}. En posant {$A = X_n$} et {$B = \{X_0,\ldots,X_{n-1}\}$}, on obtient donc:

{$$P(X_0,\ldots,X_n) = P(X_n | X_0,\ldots,X_{n-1}) P(X_0,\ldots,X_{n-1}).$$}

On peut réitérer cette opération pour le terme de droite en posant {$A = X_{n-1}$} et {$B=\{X_0,\ldots,X_{n-2}\}$}, et ainsi de suite. Donc, par récurrence, on a:

{$$P(X_0,\ldots,X_n) = P(X_0) \times \prod_{i=1}^n P(X_i | X_0,\ldots,X_{i-1} ).$$}

Si on applique à chaque terme {$P(X_i | X_0,\ldots,X_{i-1} )$} la fonction compact_proba, on obtient une décomposition:

{$$P(X_0,\ldots,X_n) = P(X_0) \times \prod_{i=1}^n P(X_i | {\cal K_i}),$$}

avec {$K_i \subseteq \{X_0,\ldots,X_{i-1}\}$}. Cette décomposition est dite compacte car son stockage nécessite en pratique beaucoup moins de mémoire que celui de la distribution jointe. C'est ce que l'on appelle un réseau bayésien.

Ecrivez une fonction create_bayesian_network : Potential x float -> Potential list qui, étant donné une probabilité jointe et un nombre réel {$\epsilon$}, vous renvoie la liste des {$P(X_i | {\cal K_i})$}. Pour cela, il vous suffit d'appliquer l'algorithme suivant:

liste = {$\emptyset$}
P = {$P(X_0,\ldots,X_n)$}
Pour {$i$} variant de manière décroissante de {$n$} à 0 faire:

fait
retourner liste

Il est intéressant ici de noter les affichages des variables de Q: comme toutes les variables sont binaires, Q nécessite uniquement (2 puissance le nombre de ces variables) nombres réels. Ainsi une probabilité sur 3 variables ne nécessite que {$2^3=8$} nombres réels.

Testez votre fonction sur la distribution jointe que vous avez lue à partir du fichier asia.txt.

6. Gain du compactage

On souhaite observer le gain en termes de consommation mémoire obtenu par votre décomposition. Si P est un Potential, alors P.toarray().size est égal à la taille (le nombre de paramètres) de la table P. Calculez donc le nombre de paramètres nécessaires pour stocker la probabilité jointe lue dans le fichier asia.txt ainsi que la somme des nombres de paramètres des tables que vous avez créées grâce à votre fonction create_bayesian_network.

Partie optionnelle

7. Applications pratiques

La technique de décomposition que vous avez vue est effectivement utilisée en pratique. Vous pouvez voir le gain que l'on peut obtenir sur différentes distributions de probabilité du site :

http://www.bnlearn.com/bnrepository/

Cliquez sur le nom du dataset que vous voulez visualiser et téléchargez son .bif ou .dsl. Afin de visualiser le contenu du fichier, vous allez utiliser pyAgrum. Le code suivant vous permettra alors de visualiser votre dataset: la valeur indiquée après "domainSize" est la taille de la probabilité jointe d'origine (en nombre de paramètres) et celle après "parameters" est la taille de la probabilité sous forme compacte (somme des tailles des probabilités conditionnelles compactes).