Ligne de niveau d'une gaussienne

Bonjour à tous
Dans le code Python ci dessous :
def f2(X1,X2,beta, Sigma):
    assert (np.linalg.det(Sigma) > 0),"Sigma is not invertible" 
    X = np.array([X1,X2])
    first = np.dot( np.linalg.inv(Sigma), (X-beta) )
    quadratic = np.dot( np.transpose(first) , X-beta) 
    den = ( (2*np.pi)**(d/2) )*( (np.linalg.det(Sigma))**(1/2) )
    return np.exp(-(quadratic)/2)/den

#Paramètres
d = 2
Sigma = np.array([5,1,1,1]).reshape(2,2)
n = 10**3
beta = [0,0]

# Génération d'un n échantillon
sample =[]
for i in range(n):
    X = stats.norm(0,1).rvs(2)
    Y = np.dot(Sigma, X) + beta
    sample.append(Y)

# Affichage pour Sigma inversible
plt.figure(figsize = [20,10])

for i in range(n):
    plt.scatter(sample[ i][0], sample[ i][1])

# Affichage des lignes de niveau
x = np.linspace(-5, 5, 10**2)
y = np.linspace(-5, 5, 10**2)

X,Y = np.meshgrid(x,y)
Z = np.zeros( [len(X),len(Y)] )
for i in range(len(x)):
    for j in range(len(y)):
        Z[ i][j] = f2(x[ i],y[j],beta, Sigma)

plt.contour(X,Y,Z, cmap='RdGy')    
plt.colorbar()

plt.show()
Le but est de tracer les lignes de niveau d'une gaussienne en dimension 2. Vous ne trouvez pas bizarre la figure obtenue ? Les lignes de niveaux sont des ellipses mais elles ne vont pas dans le sens du nuage de points. C'est normal ?

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]

Réponses

  • Bonjour,

    Non ce n'est pas normal. Au minimum l’échantillonnage me semble suspect, avec une multiplication par la matrice de covariance et non une racine.

    Je te propose à la va-vite une variante, sans erreur j’espère... Deux suggestions:
    * En numpy (matlab, etc.), utilise du code vectorisé. Fuis les boucles for, qui sont moins rapides et moins lisibles. D'ailleurs la gymnastique mentale qu'impose l’écriture tensorielle est intéressante.
    * Quand apparaissent des formes symétriques, écris du code qui numériquement exploite et/ou préserve la symétrie.
    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt
    
    def f2(X, mu, sigma_svd):
        """
        X: (N,2)
        mu: (2,)
        sigma_svd: output of np.linalg.svd on (2,2) matrix
        """
        d = len(mu)
        u, s, _ = sigma_svd
        isqrt_sigma = u / np.sqrt(s)
        Z = np.matmul(X-mu, isqrt_sigma)
        E = np.sum(Z**2, axis=-1) + np.log(s).sum() + np.log(2*np.pi)*d
        return np.exp(-E/2)
    
    #Paramètres
    Sigma = np.array([5,1,1,3]).reshape(2,2)
    N = 10**3
    mu = np.zeros(2)
    
    # Square root
    u, s, vh = sigma_svd = np.linalg.svd(Sigma, full_matrices=True, hermitian=True)
    sqrt_sigma = vh*np.sqrt(s[:,np.newaxis])
    
    # Génération d'un n échantillon
    Z = np.random.normal(size=N*2).reshape(-1,2)
    sample = np.matmul(Z, sqrt_sigma) + mu
    
    # Affichage pour Sigma inversible
    plt.figure(figsize = [10,5])
    plt.scatter(sample[:,0], sample[:,1], zorder=10, alpha=.25, color='grey', edgecolors='none')
    
    # Affichage des lignes de niveau
    res = 10**2
    x = np.linspace(-5, 5, res)
    y = np.linspace(-5, 5, res)
    dx2 = (x[-1]-x[0])/res*(y[-1]-y[0])/res
    
    X,Y = np.meshgrid(x,y) # res*res vectors
    p = f2(np.stack([X.reshape(-1),Y.reshape(-1)], axis=1), 
           mu, sigma_svd).reshape(res,res) # (res*res,) into (res,res) 
    
    plt.contourf(X,Y,p*dx2, 100, cmap='bone_r')    
    plt.colorbar()
    ax = plt.gca()
    ax.set_xlim(-5,5)
    ax.set_ylim(-5,5)
    
    plt.show()
    

    Correction. sqrt_sigma = vh*np.sqrt(s[:,np.newaxis]) et non u*s
  • Bonjour Talbon,

    Alors je vais vous répondre en deux temps, tout d'abord une question de maths puis peut être une question de code.
    Pourquoi je dis peut être ? C'est simplement que je suis intéressé par vos suggestions pour améliorer mon code :
    * En numpy (matlab, etc.), utilise du code vectorisé. Fuis les boucles for, qui sont moins rapides et moins lisibles. D'ailleurs la gymnastique mentale qu'impose l’écriture tensorielle est intéressante.
    * Quand apparaissent des formes symétriques, écris du code qui numériquement exploite et/ou préserve la symétrie.

    Et peut être que j'aurais une question lors de cette deuxième étape qui sait :)o. Mais je ne le fais pas maintenant, pourquoi ? Parce que d'abord je vais essayer de comprendre un peu mathématiquement le problème.

    Déjà vous avez mis le doigt dessus, j'ai mal simulé mon échantillon, voici une correction :
    # Ecriture de Sigma = A*A^t
    vp, vecp = np.linalg.eig(Sigma)
    P = vecp
    racD = np.zeros([2,2], dtype = float)
    racD[0][0] = np.sqrt(vp[0])
    racD[1][1] = np.sqrt(vp[1])
    A = np.dot(P,racD)
    
    # Génération d'un n échantillon
    sample =[]
    for i in range(n):
        X = stats.norm(0,1).rvs(2)
        Y = np.dot(A, X) + beta
        sample.append(Y)
    

    Je remarque que le code n'est pas optimal, on pourra supprimer la boucle for lors de la deuxième étape.

    La densité de la loi $N(\beta, \Sigma)$ est
    $$
    \boxed{
    f(x)= \frac{\exp\left( \frac{-1}{2}\langle (\Sigma)^{-1}(x-\beta), (x-\beta) \rangle \right)}{(2\pi)^{\frac{d}{2}} (\operatorname{det}(\Sigma))^{\frac{1}{2}} }
    }
    $$
    Soit $x\in \mathbb{R}^{d}$ appartient à la ligne de niveau en $c$ est
    $$
    c = \frac{\exp\left( \frac{-1}{2}\langle (\Sigma)^{-1}(x-\beta), (x-\beta) \rangle \right)}{(2\pi)^{\frac{d}{2}} (\operatorname{det}(\Sigma))^{\frac{1}{2}} }
    \Leftrightarrow \exists c' \in \mathbb{R} ; c' = \langle (\Sigma)^{-1}(x-\beta), (x-\beta) \rangle \Leftrightarrow x\in \operatorname{Fr} \left( \mathcal{E}(c') \right)
    $$
    avec $\mathcal{E}(k)$ l'intérieur d'une ellipsoide centrée en $\beta$. Et le truc qui est fou c'est que pour $k$ à partir d'un certain rang on arrive à avoir 95% des points dans l'ellipsoïde.

    Il suffit d'observer le code suivant :
    #Paramètres
    d = 2
    Sigma = np.array([5,1,1,1]).reshape(2,2)
    n = 10**3
    beta = [0,0]
    
    # Génération d'un n échantillon
    sample = stats.multivariate_normal(beta, Sigma).rvs(size = n)
    
    # Affichage pour Sigma inversible
    plt.figure(figsize = [20,10])
    for i in range(n):
        plt.scatter(sample[ i][0], sample[ i][1])
    
    # Affichage de l'ellipse
    k = stats.chi2(2).ppf(0.95)
    
    # Test
    test = True
    if test :
        X = np.linspace(-10, 10, 10**2)
        Y = np.linspace(-10, 10, 10**2)
        for x in X:
            for y in Y:
                Sigma_inv = np.linalg.inv(Sigma)
                u = [x,y]
                if np.dot( np.transpose(np.dot(Sigma_inv,u)),  u) <= k:
                    plt.scatter(x,y, color = "black")
    
    plt.show()
    
    

    J'admet que ce n'est pas la bonne méthode pour tracer une ellipse, il vaut mieux donner une équation paramétrique mais ici c'était juste pour montrer qu'on a bien une ellipse.

    Comme on vient de le voir, il n'y a pas de doute les lignes de niveaux devraient être des ellipses centrées en beta. Le fait que quand on utilise plt.contour on n'obtient des ellipses de confiance disons pas dans le bon sens, est toujours un mystère pour moi mathématiquement.
  • Oui, les lignes de niveau sont bien des ellipses / ellipsoïdes centrés sur la moyenne, avec le théorème spectral qui te donne une information géométrique (orientation spatiale et demi-axes).

    Si les ellipses obtenues via plt.contour ne correspondent pas au nuage de points, c'est que le code contient une autre erreur en amont (ça marche dans la variante que j'ai proposée ; pour t'en convaincre tu peux remplacer contourf par contour, ou bien réduire le nombre de niveaux de $100$ à $20$ ou $10$).

    Tiens, d'ailleurs je vois peut-être une erreur. Le dernier bloc de code avec $X,Y,Z$ et $x,y$ n'est peut-être pas cohérent, quant à l'ordre dans lequel tu remplis $Z$ et l'ordre dans lequel $X,Y$ parcourent la grille.
  • Oui bien vu ! Tu as trouvé l'erreur. Le code suivant fonctionne mais n'est pas optimal.
    #Paramètres
    d = 2
    Sigma = np.array([5,1,1,1]).reshape(2,2)
    n = 10**3
    beta = [5,5]
    
    # Ecriture de Sigma = A*A^t
    vp, vecp = np.linalg.eig(Sigma)
    P = vecp
    racD = np.zeros([2,2], dtype = float)
    racD[0][0] = np.sqrt(vp[0])
    racD[1][1] = np.sqrt(vp[1])
    A = np.dot(P,racD)
    
    # Génération d'un n échantillon
    sample =[]
    for i in range(n):
        X = stats.norm(0,1).rvs(2)
        Y = np.dot(A, X) + beta
        sample.append(Y)
    
    # Affichage pour Sigma inversible
    plt.figure(figsize = [20,10])
    
    for i in range(n):
        plt.scatter(sample[ i][0], sample[ i][1])
    
    # Affichage des lignes de niveau
    x = np.linspace(beta[0] - 5 ,beta[1] + 5, 10**2)
    y = np.linspace(beta[0] - 5 ,beta[1] + 5, 10**2)
    
    X,Y = np.meshgrid(x,y)
    Z = np.zeros( [len(X),len(Y)] )
    for i in range(len(x)):
        for j in range(len(y)):
            Z[ i][j] = stats.multivariate_normal(beta, Sigma).pdf([X[ i,j],Y[ i,j]])
    
    plt.contour(X,Y,Z,levels=14, cmap="RdBu_r")
    plt.colorbar()
    plt.show()
    
    En fait ça venait de ma mauvaise compréhension de plt.meshgrid ! Je vais regarder sur internet mais j'ai du mal à comprendre comment ça marche ce truc :)o.
  • Pour commencer l'étape 2 je vais essayer de virer les boucles for :
    - * En numpy (matlab, etc.), utilise du code vectorisé. Fuis les boucles for, qui sont moins rapides et moins lisibles. D'ailleurs la gymnastique mentale qu'impose l’écriture tensorielle est intéressante.
  • A quoi correspond le -1 dans reshape(-1,2) s'il vous plaît ?
  • Comment marche np.matmul(Z, Sigma) ? Déjà je trouve bizarre qu'on ne fasse pas Sigma,Z après Sigma est une matrice 2x2 et Z est un ensemble de vecteur de taille 2. Donc je pense qu'on remplace chaque vecteur de taille 2 par un nouveau vecteur de taille 2 qui est le vecteur obtenu par multiplication par Z.
  • Le $-1$ dans reshape est remplacé par l'entier positif tel que le nombre d'éléments du vecteur soit le même avant et après redimensionnement. Si $v$ est $3\times 5\times 6$ et que j’écris v.reshape(-1,9), c'est équivalent à v.reshape(10,9). v.reshape(-1,7) renverrait une erreur, car $7$ ne divise pas le produit des dimensions. v.reshape(-1,-1,5) n'est pas non plus autorisé, car c'est intrinsèquement ambigu. v.reshape(2,-1,5) marche et est équivalent à v.reshape(2,9,5). (Souvent reshape ne coûte rien, c'est une autre façon de voir le même tableau en mémoire. Bon, ici ce n'est pas important.)

    Pour le fonctionnement général de np.matmul, je te renvoie à la documentation. Pour le cas particulier np.matmul($Z, C$) où les deux arguments sont $2$D, c'est simplement un produit matriciel. Par exemple si $Z$ est $N\times 2$ et $C$ est $2\times 2$, cela renvoie $ZC$ qui est encore $N\times 2$. Chaque ligne $(ZC)_i$ (de taille $1\times 2$) correspond au produit de la ligne $Z_i$ ($1\times 2$) de $Z$ par la matrice $C$ ($2\times2$).

    Dans ma version j'ai un $X-\mu$ qui est $N\times 2$ : chaque ligne $\bar{x}_i\triangleq (X-\mu)_i$ est un point de donnée centré. On veut les $\bar{x}_i^T \sqrt{\Lambda}\sqrt{\Lambda}^T \bar{x}_i$ qui sont $N$ scalaires. Ils s’écrivent $\Vert \sqrt{\Lambda}^T \bar{x}_i \Vert_{\mathbb{R}^{2}}^2$, mais $\Vert \bar{x}_i^T \sqrt{\Lambda}\Vert_{\mathbb{R}^{1\times 2}}^2$ est tout aussi vrai. C'est cette deuxième forme que j'utilise car elle tombe tout cuit, pour tous les points avec np.matmul($X-\mu, \sqrt{\Lambda})$, et ce sans avoir à faire la moindre transposition.
  • D'accord j'ai compris pour le -1 dans reshape merci.
  • Pour le produit matriciel aussi j'ai compris vous avez transposé le produit AX. Ça demande effectivement une gymnastique parce qu'il faut en plus ajouter des lignes à X transposé et on récupère un vecteur dont les lignes sont le résultat.
  • Pouvez vous m'expliquer s'il vous plaît comment vous vectorisez
    for i in range(len(x)):
        for j in range(len(y)):
            Z[ i][j] = stats.multivariate_normal(beta, Sigma).pdf([X[ i,j],Y[ i,j]])
    
  • édit. J'avais mal lu la question. :)o Tu as une réponse partielle dans le dernier paragraphe de mon précédent message. On veut les $E_i = -\log{p(x_i)} = \frac{1}{2}\left((x_i-\beta)^T \Sigma^{-1} (x_i-\beta) + \log{|2\pi \Sigma|}\right)$. $E_i\mapsto \exp(-E_i)$ donne les probabilités si nécessaire.

    On exploite le théorème spectral : $\Sigma^{-1} = U D^{-1} U^T = (U D^{-1/2}) (U D^{-1/2})^T$.
    Pour $U$ et $D$, faire par exemple : u, s, _ = np.linalg.svd(Sigma, hermitian=True). Alors $U \equiv$ u et $D\equiv$ s.

    Remplace $\sqrt{\Lambda}$ dans le précédent post par $U D^{-1/2}$ ($\equiv$ U/np.sqrt(s)) pour vectoriser les $(x_i-\beta)^T \Sigma^{-1} (x_i-\beta)$.

    Le deuxième terme $\log{|2\pi \Sigma|}$ est une constante de $x$.
    On peut calculer ce scalaire comme $\sum_i \log{D_{ii}} + d \log{2\pi}$ ($\equiv$ np.log(s).sum()+np.log(2*np.pi)*d).
Connectez-vous ou Inscrivez-vous pour répondre.