<jupyter_start><jupyter_text>MAPSI - TME  - Rappels de Proba/stats  I- La planche de Galton ( obligatoire)  I.1- Lo<NAME> BernouilliÉcrire une fonction `bernouilli: float ->int` qui prend en argument la paramètre $p \in [0,1]$ et qui renvoie aléatoirement $0$ (avec la probabilité $1-p$) ou $1$ (avec la probabilité $p$).<jupyter_code>def bernouilli(p):
    if(np.random.random(1)>p):
        return 1
    else:
        return 0<jupyter_output><empty_output><jupyter_text>I.2- Loi binomialeÉcrire une fonction `binomiale: int , float -> int` qui prend en argument un entier $n$ et $p \in [0,1]$ et qui renvoie aléatoirement un nimbre tiré selon la distribution ${\cal B}(n,p)$.<jupyter_code>def binomiale(n, p):
    acc=0
    for i in range(n):
        acc += bernouilli(p)
    return acc<jupyter_output><empty_output><jupyter_text>I.3- Histogramme de la loi binomialeDans 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-contre.  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)$. Écrire un script qui crée un tableau 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, tracer l'histogramme de ce tableau. Vous pourrez utiliser la fonction hist de matplotlib.pyplot:<jupyter_code>import matplotlib.pyplot as plt

nb_billes = 10000
nb_etape = 100
proba = 0.5
tab = np.array([])

for i in range(nb_billes):
    tab = np.append(tab, binomiale(nb_etape, proba))

plt.hist(tab, nb_etape, width=1.1)
plt.title("Histogramme répartition des billets dans la planche de Galton")
plt.show()<jupyter_output><empty_output><jupyter_text>Pour le nombre de bins, calculez le nombre de valeurs différentes dans votre tableau.<jupyter_code>nb_unique = np.size(np.unique(tab))
print("Nombre de barres de l'histogramme:", nb_unique)<jupyter_output>Nombre de barres de l'histogramme: 37<jupyter_text>II- Visualisation d'indépendances ( obligatoire)  II.1- Loi normale centrée réduiteOn 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σ$ à $2σ$, les $y_i$ correspondant à la valeur de la fonction de densité de la loi normale centrée de variance $σ^2$, autrement dit ${\cal N}(0,σ^2)$.Écrire une fonction `normale : int , 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.<jupyter_code>import random
import math

def normale(k, sigma):
    if(k%2==0):
        raise ValueError("le nombre k doit etre impair")
    else:
        res = np.zeros(shape=(k, 2))
        x = np.linspace(-2*sigma, 2*sigma, k)
        for i in range(k):
            x_val = x[i]
            y = (np.exp(-0.5*(x_val/sigma)**2))/(sigma*np.sqrt(2*np.pi))
#             print("x:",x)
#             print("y:",y)
            res[i][0] = x_val
            res[i][1] = y
        return res<jupyter_output><empty_output><jupyter_text>Vérfier la validité de votre fonction en affichant grâce à la fonction plot les points générés dans une figure.<jupyter_code>nb_points = 2001
sigma = 10
points = normale(nb_points, sigma)

plt.plot(points[:,0], points[:,1])
plt.title("Loi normale entre -2 sigma et +2 sigma")
plt.show()<jupyter_output><empty_output><jupyter_text>II.2- Distribution de probabilité affineDans 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-contre. Écrire une fonction `proba_affine : int , 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 vérifiera ici aussi 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 à $\frac{1}{k}$ (afin que $\sum y_i=1$). Si la pente est différente de $0$, il suffit de choisir, $\forall i=0,...,k−1$,$$y_i=\frac{1}{k}+(i−\frac{k−1}{2})×slope$$Vous pourrez aisément vérifier que, ici aussi, $\sum y_i=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.<jupyter_code>def proba_affine(k, slope):
    if k % 2 == 0:
        raise ValueError("le nombre k doit etre impair")
    if abs ( slope  ) > 2. / ( k * k ):
        raise ValueError("la pente est trop raide : pente max = " +
        str ( 2. / ( k * k ) ) )
    res = np.zeros(shape=(k, 2))
    for i in range(k):
        res[i][0] = i
        res[i][1] = (1/k) + (i - ((k-1)/2)) * slope
    return res
nb_points = 101
slope = 1e-4
points = proba_affine(nb_points, slope)
plt.plot(points[:,0], points[:,1])
plt.title("Courbe distribution de probabilité affine")
plt.show()<jupyter_output><empty_output><jupyter_text>II.3- Distribution jointeÉcrire une fonction `Pxy : float np.array , 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:<jupyter_code>PA = np.array ( [0.2, 0.7, 0.1] )
PB = np.array ( [0.4, 0.4, 0.2] )<jupyter_output><empty_output><jupyter_text>alors `Pxy(A,B)` renverra le tableau :```np.array([[ 0.08,  0.08,  0.04],          [ 0.28,  0.28,  0.14],          [ 0.04,  0.04,  0.02]])```<jupyter_code>def Pxy(x,y):
    n1 = np.size(x)
    n2 = np.size(y)
    res = np.zeros(shape=(n1, n2))
    for i in range(n1):
        for j in range(n2):
            res[i][j] = x[i]*y[j]
    return res
print("Distribution jointe de PA PB:\n",Pxy(PA, PB))<jupyter_output>Distribution jointe de PA PB:
 [[0.08 0.08 0.04]
 [0.28 0.28 0.14]
 [0.04 0.04 0.02]]<jupyter_text>II.4- Affichage de la distribution jointeLe 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 la commande `%matplotlib notebook` fonctione, vous pouvez interagir avec la courbe. 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.<jupyter_code>from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline  
# essayer `%matplotib notebook` pour interagir avec la visualisation 3D

def dessine ( P_jointe ):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    x = np.linspace ( -3, 3, P_jointe.shape[0] )
    y = np.linspace ( -3, 3, P_jointe.shape[1] )
    X, Y = np.meshgrid(x, y)
    ax.plot_surface(X, Y, P_jointe, rstride=1, cstride=1 )
    ax.set_xlabel('A')
    ax.set_ylabel('B')
    ax.set_zlabel('P(A) * P(B)')
    plt.show ()
dessine(Pxy(PA,PB))
k = 101
sigma = 10
slope = 1e-5
dessine(Pxy(normale(k, sigma)[1], proba_affine(k, slope)[1]))<jupyter_output><empty_output><jupyter_text>III- Indépendances conditionnelles ( obligatoire)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 :<jupyter_code># creation de P(X,Y,Z,T)                     
P_XYZT = np.array([[[[ 0.0192,  0.1728],
                     [ 0.0384,  0.0096]],    
                                             
                    [[ 0.0768,  0.0512],     
                     [ 0.016 ,  0.016 ]]],   
                                             
                   [[[ 0.0144,  0.1296],  
                     [ 0.0288,  0.0072]],    
                                             
                    [[ 0.2016,  0.1344],     
                     [ 0.042 ,  0.042 ]]]])<jupyter_output><empty_output><jupyter_text>Ainsi, $\forall (x,y,z,t) \in \{0,1\}^4$, `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)$.   III.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)$. Il s'agit donc de vérifer que dans la loi $P$, $$P(X,T|Y,Z)=P(X|Y,Z)\cdot P(T|Y,Z)$$Pour cela, tout d'abord, calculer à 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)s, vous devez obtenir le tableau suivant : ```np.array([[ 0.336,  0.084],          [ 0.464,  0.116]])```Ainsi $P(Y=0,Z=1)=$ `P_YZ[0][1]` $=0.084$<jupyter_code>P_YZ = np.zeros((2,2))
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            for t in range(P_XYZT.shape[3]):
                P_YZ[y][z] += P_XYZT[x][y][z][t]

print(P_YZ)<jupyter_output>[[0.336 0.084]
 [0.464 0.116]]<jupyter_text>Ensuite, calculer 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)=\frac{P(X,Y,Z,T)}{P(Y,Z)}$$<jupyter_code>P_XTcondYZ = np.zeros((2,2,2,2))
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            for t in range(P_XYZT.shape[3]):
                P_XTcondYZ[x][y][z][t] += P_XYZT[x][y][z][t] / P_YZ[y][z]

print(P_XTcondYZ)<jupyter_output>[[[[0.05714286 0.51428571]
   [0.45714286 0.11428571]]

  [[0.16551724 0.11034483]
   [0.13793103 0.13793103]]]


 [[[0.04285714 0.38571429]
   [0.34285714 0.08571429]]

  [[0.43448276 0.28965517]
   [0.36206897 0.36206897]]]]<jupyter_text>Calculer à 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)=∑_Y P(X,T|Y,Z)$$<jupyter_code>P_XcondYZ = np.zeros((2,2,2))
P_TcondYZ = np.zeros((2,2,2))
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            for t in range(P_XYZT.shape[3]):
                P_XcondYZ[x][y][z] += P_XTcondYZ[x][y][z][t]
                P_TcondYZ[y][z][t] += P_XTcondYZ[x][y][z][t]

print(P_XcondYZ)
print(P_TcondYZ)<jupyter_output>[[[0.57142857 0.57142857]
  [0.27586207 0.27586207]]

 [[0.42857143 0.42857143]
  [0.72413793 0.72413793]]]
[[[0.1 0.9]
  [0.8 0.2]]

 [[0.6 0.4]
  [0.5 0.5]]]<jupyter_text>Enfin, tester 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)×P(T|Y,Z)$$<jupyter_code># A = np.round(P_XcondYZ * P_TcondYZ, 5) # Erreur !
A = np.zeros((2, 2, 2, 2))
ecart = 0
epsilon = 1e-3
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            for t in range(P_XYZT.shape[3]):
#                 A[x][y][z][t] = P_XcondYZ[x][y][z] * P_TcondYZ[y][z][t]
                ecart += np.abs((P_XcondYZ[x][y][z] * P_TcondYZ[y][z][t]) - P_XTcondYZ[x][y][z][t]) 
    
if ecart < epsilon:
    print("indépendant")
else:
    print("pas indépendant")
    
# A = np.round(A, 5)
# B = np.round(P_XTcondYZ, 5)
# print(A.shape)
# print("=======")
# print(B.shape)
# print(A == B)<jupyter_output>indépendant<jupyter_text>III.2- Indépendance de X et (Y,Z)On souhaite maintenant déterminer si $X$ et $(Y,Z)$ sont indépendantes. Pour cela, commencer par calculer à partir de `P_XYZT` le tableau `P_XYZ` représentant la distribution $P(X,Y,Z)$. Ensuite, calculer à 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)=∑_Y∑_Z P(X,Y,Z)$$Si vous ne vous êtes pas trompé(e), P_X doit être égal au tableau suivant : ```np.array([ 0.4,  0.6])```<jupyter_code>P_XYZ = np.zeros((2,2,2))
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            for t in range(P_XYZT.shape[3]):
                P_XYZ[x][y][z] += P_XYZT[x][y][z][t]
print(P_XYZ)

# P_X = np.zeros((2))
# for x in range(P_XYZT.shape[0]):
#     for y in range(P_XYZT.shape[1]):
#         for z in range(P_XYZT.shape[2]):
#             for t in range(P_XYZT.shape[3]):
#                 P_X[x] += P_XYZT[x][y][z][t]
# print(P_X)

P_X = np.zeros((2))
for x in range(P_XYZT.shape[0]):
    for y in range(P_XYZT.shape[1]):
        for z in range(P_XYZT.shape[2]):
            P_X[x] += P_XYZ[x][y][z]
print(P_X)<jupyter_output>[[[0.192 0.048]
  [0.128 0.032]]

 [[0.144 0.036]
  [0.336 0.084]]]
[0.4 0.6]<jupyter_text>Enfin, si $X$ et $(Y,Z)$ sont bien indépendantes, on doit avoir $$P(X,Y,Z)=P(X)×P(Y,Z)$$<jupyter_code>A = np.zeros((2, 2, 2))
ecart = 0
epsilon = 1e-3
for x in range(P_XYZ.shape[0]):
    for y in range(P_XYZ.shape[1]):
        for z in range(P_XYZ.shape[2]):
            A[x][y][z] = P_X[x] * P_YZ[y][z]
            ecart += np.abs(P_X[x] * P_YZ[y][z] - P_XYZ[x][y][z])

if ecart < epsilon:
    print("indépendant")
else:
    print("pas indépendant")
            
# A = np.round(A, 5)
# B = np.round(P_XYZ, 5)
# print(A.shape)
# print("=======")
# print(B.shape)
# print(A == B)<jupyter_output>pas indépendant<jupyter_text>IV- Indépendances conditionnelles et consommation mémoire  ( obligatoire)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,…,X_n$ de variables aléatoires binaires (elles ne peuvent prendre que 2 valeurs : 0 et 1).  Simplification du code : utilisation de pyAgrumManipuler des probabilités et des opérations sur des probabilités complexes est difficiles avec les outils classiques. La difficulté principale est certainement le problème du mapping entre axe et variable aléatoire. `pyAgrum` propose une gestion de `Potential` qui sont des tableaux multidimensionnels dont les axes sont caractérisés par des variables et sont donc non ambigüs.Par exemple, après l'initiation du `Potential PABCD` :<jupyter_code>import pyAgrum as gum
import pyAgrum.lib.notebook as gnb

X,Y,Z,T=[gum.LabelizedVariable(x,x,2) for x in "XYZT"]
pXYZT=gum.Potential().add(T).add(Z).add(Y).add(X)
pXYZT[:]=[[[[ 0.0192,  0.1728],
            [ 0.0384,  0.0096]],
           [[ 0.0768,  0.0512],
            [ 0.016 ,  0.016 ]]],
          [[[ 0.0144,  0.1296],
            [ 0.0288,  0.0072]],
           [[ 0.2016,  0.1344],
            [ 0.042 ,  0.042 ]]]]<jupyter_output><empty_output><jupyter_text>On peut alors utiliser la méthode `margSumOut` qui supprime les variables par sommations: `p.margSumOut(['X','Y'])` correspond à calculer $\sum_{X,Y} p$La réponse a question III.1 se calcule donc ainsi :<jupyter_code>pXT_YZ=pXYZT/pXYZT.margSumOut(['X','T'])
pX_YZ=pXT_YZ.margSumOut(['T'])
pT_YZ=pXT_YZ.margSumOut(['X'])

if pXT_YZ==pX_YZ*pT_YZ:
    print("=> X et T sont indépendants conditionnellemnt à Y et Z")
else:
    print("=> pas d'indépendance trouvée")<jupyter_output>=> X et T sont indépendants conditionnellemnt à Y et Z<jupyter_text>La réponse à la question III.2 se calcule ainsi :<jupyter_code>pXYZ=pXYZT.margSumOut("T")
pYZ=pXYZ.margSumOut("X")
pX=pXYZ.margSumOut(["Y","Z"])
if pXYZ==pX*pYZ:
    print("=> X et YZ sont indépendants")
else:
    print("=> pas d'indépendance trouvée")
gnb.sideBySide(pXYZ,pX,pYZ,pX*pYZ,
               captions=['$P(X,Y,Z)$','$P(X)$','$P(Y,Z)$','$P(X)\cdot P(Y,Z)$'])<jupyter_output><empty_output><jupyter_text>`asia.txt` contient la description d'une probabilité jointe sur un ensemble de $8$ variables aléatoires binaires (256 paramètres). Le fichier est produit à partir du site web suivant `http://www.bnlearn.com/bnrepository/`.Le code suivant permet de lire ce fichier et d'en récupérer la probabilité jointe (sous forme d'une `gum.Potential`) qu'il contient :<jupyter_code>def read_file ( filename ):
    """
    Renvoie les variables aléatoires et la probabilité contenues dans le
    fichier dont le nom est passé en argument.
    """
    Pres = gum.Potential ()
    vars=[]

    with open ( filename, 'r' ) as fic:
        # on rajoute les variables dans le potentiel
        nb_vars = int ( fic.readline () )
        for i in range ( nb_vars ):
            name, domsize = fic.readline ().split ()
            vars.append(name)
            variable = gum.LabelizedVariable(name,name,int (domsize))
            Pres.add(variable)

        # on rajoute les valeurs de proba dans le potentiel
        cpt = []
        for line in fic:
            cpt.append ( float(line) )
        Pres.fillWith( cpt )
    return vars,Pres

vars,Pjointe=read_file('asia.txt')
# afficher Pjointe est un peu délicat (retire le commentaire de la ligne suivante)
# Pjointe

print('Les variables : '+str(vars))
# Noter qu'il existe une fonction margSumIn qui, à l'inverse de MargSumOut, élimine 
# toutes les variables qui ne sont pas dans les arguments
Pjointe.margSumIn(['tuberculosis?','lung_cancer?'])<jupyter_output><empty_output><jupyter_text>IV.1- test d'indépendance conditionnelleEn utilisant la méthode `margSumIn` (voir juste au dessus), écrire une fonction `conditional_indep: Potential,str,str,list[str]->bool` qui rend vrai si dans le `Potential`, on peut lire l'indépendance conditionnelle.Par exemple, l'appel`conditional_indep(Pjointe,'bronchitis?', 'positive_Xray?',['tuberculosis?','lung_cancer?'])` vérifie si bronchitis est indépendant de `posititve_Xray` conditionnellement à `tuberculosis?` et `lung_cancer?`D'un point de vue général, on vérifie que $X$ et $Y$ sont indépendants conditionnellement à $Z_1,\cdots,Z_d$ par l'égalité :$$P(X,Y|Z_1,\cdots,Z_d)=P(X|Z_1,\cdot,Z_d)\cdot P(Y|Z_1,\cdots,Z_d)$$Ces trois probabilités sont calculables à partir de la loi jointe de $P(X,Y,Z_1,\cdots,Z_d)$.Remarque Vérifier l'égalité `P==Q` de 2 `Potential` peut être problématique si les 2 sont des résultats de calcul : il peut exister une petite variation. Un meilleur test est de vérifier `(P-Q).abs().max()<epsilon` avec `epsilon` assez petit.<jupyter_code>epsilon = 1e-3

# def conditional_indep(P,X,Y,Zs): # VERIFIER
#     if len(Zs) != 0:
#         pX_Zs = P.margSumIn(X) / P.margSumIn(Zs) 
#         pY_Zs = P.margSumIn(Y) / P.margSumIn(Zs)
#         PXY_Zs = P.margSumIn([X, Y]) / P.margSumIn(Zs) 
#         A = pX_Zs * pY_Zs
#         B = PXY_Zs
#         print("A:",A)
#         print("B:",B)
#         print("A-B",A-B)
#         test = (A-B).abs().max() < epsilon
#         if(test):
#             return "=> X et Y sont indépendants conditionnellemnt à Zs"
#         else:
#             return "=> X et Y ne sont pas indépendants conditionnellemnt à Zs"
#     else:
#         pX = P.margSumIn(X)
#         pY = P.margSumIn(Y)
#         PXY = P.margSumIn([X, Y]) 
#         A = pX * pY
#         B = PXY
#         print("A:",A)
#         print("B:",B)
#         print("A-B",A-B)
#         test = (A-B).abs().max() < epsilon
#         if(test):
#             return "=> X et Y sont indépendants conditionnellemnt à Zs"
#         else:
#             return "=> X et Y ne sont pas indépendants conditionnellemnt à Zs"
      
def conditional_indep(P,X,Y,Zs):
    variables = [X, Y]
    for var in Zs:
        variables.append(var)
    if len(Zs) != 0:
        pXYZs = P.margSumIn(variables) # Besoin de cette valeur pour calculer PXY_ZS
        pXY_Zs = pXYZs/pXYZs.margSumOut([X, Y]) # On enleve X Y pour avoir X et Y conditionnellement à Zs
        pX_Zs = pXY_Zs.margSumOut([Y])
        pY_Zs = pXY_Zs.margSumOut([X])
        A = pX_Zs * pY_Zs
        B = pXY_Zs
#         print("A:",A)
#         print("B:",B)
#         print("A-B",A-B)
        test = (A-B).abs().max() < epsilon
        if(test):
            return 1 # Independant
        else:
            return 0 # Pas independant
    else:
        pX = P.margSumIn(X)
        pY = P.margSumIn(Y)
        PXY = P.margSumIn([X, Y]) 
        A = pX * pY
        B = PXY
#         print("A:",A)
#         print("B:",B)
#         print("A-B",A-B)
        test = (A-B).abs().max() < epsilon
        if(test):
            return 1
        else:
            return 0   
conditional_indep(Pjointe,
                  'bronchitis?', 
                  'positive_Xray?',
                  ['tuberculosis?','lung_cancer?'])
conditional_indep(Pjointe,
                  'bronchitis?', 
                  'visit_to_Asia?',
                  [])<jupyter_output><empty_output><jupyter_text>IV.2- Factorisation compacte de loi jointeOn 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. Écrire une fonction `compact_conditional_proba: Potential,str-> Potential` qui, étant donné une probabilité jointe $P(X_{i_0},\ldots,X_{i_n})$, une variable aléatoire $X_{i_n}$, retourne cette probabilité conditionnelle $P(X_{i_n} | {\cal K})$. Pour cela, nous vous proposons l'algorithme itératif suivant:```K=SPour tout X in K:  Si X indépendante de Xin conditionnellement à K\{X) alors    Supprimer X de Kretourner P(Xin|K)$```Trois petites aides :1- La fonction precédente `conditional_indep` devrait vous servir...2- Obtenir la liste des noms des variables dans un `Potential` se fait par l'attribut  ```P.var_names```3- Afin que l'affichage soit plus facile à comprendre, il peut ê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 : ```proba = proba.putFirst(Xin)```<jupyter_code>def compact_conditional_proba(P, X):
    K = P.var_names # Ensemble des variables dans la Pjointe
    K.remove(X)
    
    for k in K:
        tmp_K = K.copy()
        tmp_K.remove(k)
        if conditional_indep(P, k, X, K): # Check si les variables dans K sont independantes
            K.remove(k) # Supprime de la liste si independantes
        
    var = []
    var.append(X)
    for k in K:
        var.append(k)
    PX_K = P.margSumIn(var) / P.margSumIn(K) 
    PX_K = PX_K.putFirst(X)
    
    return PX_K  
compact_conditional_proba(Pjointe,"visit_to_Asia?")
compact_conditional_proba(Pjointe,"dyspnoea?")<jupyter_output><empty_output><jupyter_text>IV.3- Création d'un réseau bayésienUn 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_conditional_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.Écrire une fonction `create_bayesian_network : Potential -> Potential list` qui, étant donné une probabilité jointe, vous renvoie la liste des $P(X_i | {\cal K_i})$. Pour cela, il vous suffit d'appliquer l'algorithme suivant:```liste = [] P = P(X_0,...,X_n)Pour i de n à 0 faire:  calculer Q = compact_conditional_proba(P,X_i)  afficher la liste des variables de Q  rajouter Q à liste  supprimer X_i de P par marginalisationretourner 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.<jupyter_code>def create_bayesian_network(P):
    res = []
    liste_names = P.var_names
    for X_i in reversed(liste_names):
        print(X_i)
        Q = compact_conditional_proba(P, X_i)
#         print(Q)
        res.append(Q)
        P = P.margSumOut(X_i)
    return res
create_bayesian_network(Pjointe)<jupyter_output>visit_to_Asia?
tuberculosis?
smoking?
lung_cancer?
tuberculosis_or_lung_cancer?
bronchitis?
positive_Xray?
dyspnoea?<jupyter_text>IV.4- Gain en compressionOn 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`.<jupyter_code># votre code<jupyter_output><empty_output><jupyter_text>V- Applications pratiques (optionnelle) 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 "dim" est la taille de la probabilité sous forme compacte (somme des tailles des probabilités conditionnelles compactes).<jupyter_code># chargement de pyAgrum
import pyAgrum as gum
import pyAgrum.lib.notebook as gnb

# chargement du fichier bif ou dsl
bn = gum.loadBN ( "asia.bif" )

# affichage de la taille des probabilités jointes compacte et non compacte
print(bn)

# affichage graphique du réseau bayésien
bn<jupyter_output>BN{nodes: 8, arcs: 8, domainSize: 256, dim: 36}
{'fraction_non_alphanumeric': 0.13117406119849928, 'fraction_numerical': 0.028017760644339654, 'mean_word_length': 3.6242240967690593, 'pattern_counts': {'":': 0, '<': 102, '<?xml version=': 0, '>': 114, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '384299', 'license': 'MIT', 'n_tokens_mistral': 11616, 'n_tokens_neox': 11064, 'n_words': 3169, 'path': 'S1/MAPSI/TME/TME2/enonce/.ipynb_checkpoints/TME2Durand-checkpoint.ipynb', 'repo_name': 'hanzopgp/MasterBioInformaticsPTs', 'size': '29051'}