Ligne de niveau d'une gaussienne
dans Statistiques
Bonjour à tous
Dans le code Python ci dessous :
[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]
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]
Connectez-vous ou Inscrivez-vous pour répondre.
Réponses
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
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.
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.
#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.- * 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.
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.
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]])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).