Ligne de niveau d'une gaussienne

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 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.

A quoi correspond le -1 dans reshape(-1,2) s'il vous plaît ?

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.