Commit d2c722b7 authored by vincentvigon's avatar vincentvigon
Browse files

vect alea3

parent 285235d3
No preview for this file type
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
# Vecteurs aléatoires, notamment gaussiens # Vecteurs aléatoires, notamment gaussiens
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
from mpl_toolkits import mplot3d from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
%matplotlib inline %matplotlib inline
np.set_printoptions(precision=2,suppress=True) np.set_printoptions(precision=2,suppress=True)
import scipy.stats as stats import scipy.stats as stats
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Préliminaire: Décomposition en valeurs singulières (SVD) ## Préliminaire: Décomposition en valeurs singulières (SVD)
### Cas général ### Cas général
C'est une factorisation en trois facteurs ayant des propriétés particulières (observez, déduisez). C'est une factorisation en trois facteurs ayant des propriétés particulières (observez, déduisez).
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
A=np.array([[1,2,3],[4,5,6]]) A=np.array([[1,2,3],[4,5,6]])
U,S,V = np.linalg.svd(A) U,S,V = np.linalg.svd(A)
print("A\n",A) print("A\n",A)
print("\nU\n",U) print("\nU\n",U)
print("\nS\n",S) print("\nS\n",S)
print("\nV\n",V) print("\nV\n",V)
diag=np.zeros((2,3)) diag=np.zeros((2,3))
diag[0,0]=S[0] diag[0,0]=S[0]
diag[1,1]=S[1] diag[1,1]=S[1]
print("\nU@diag@V\n",U@diag@V) print("\nU@diag@V\n",U@diag@V)
print("\nU@U.T\n",U@U.T) print("\nU@U.T\n",U@U.T)
print("\nV@V.T\n",V@V.T) print("\nV@V.T\n",V@V.T)
``` ```
%% Output %% Output
A A
[[1 2 3] [[1 2 3]
[4 5 6]] [4 5 6]]
U U
[[-0.39 -0.92] [[-0.39 -0.92]
[-0.92 0.39]] [-0.92 0.39]]
S S
[9.51 0.77] [9.51 0.77]
V V
[[-0.43 -0.57 -0.7 ] [[-0.43 -0.57 -0.7 ]
[ 0.81 0.11 -0.58] [ 0.81 0.11 -0.58]
[ 0.41 -0.82 0.41]] [ 0.41 -0.82 0.41]]
U@diag@V U@diag@V
[[1. 2. 3.] [[1. 2. 3.]
[4. 5. 6.]] [4. 5. 6.]]
U@U.T U@U.T
[[1. 0.] [[1. 0.]
[0. 1.]] [0. 1.]]
V@V.T V@V.T
[[ 1. -0. 0.] [[ 1. -0. 0.]
[-0. 1. -0.] [-0. 1. -0.]
[ 0. -0. 1.]] [ 0. -0. 1.]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Interprétation géométrique: ### Interprétation géométrique:
Considérons application linéaire $T : \mathbb R^n \to \mathbb R^m$, associée à la matrice $A$. On peut trouver une base orthonormale de $\mathbb R^n$ et une base orthonormale pour $\mathbb R^m$ telles que $T$ associe au i-ème vecteur de base de $\mathbb R^n$ un multiple positif du i-ème vecteur de base de $\mathbb R^n$, les vecteurs restants ayant pour image 0. Dans ces bases, l'application $T$ est donc représentée par une matrice diagonale dont les coefficients sont des réels positifs. Considérons application linéaire $T : \mathbb R^n \to \mathbb R^m$, associée à la matrice $A$. On peut trouver une base orthonormale de $\mathbb R^n$ et une base orthonormale pour $\mathbb R^m$ telles que $T$ associe au i-ème vecteur de base de $\mathbb R^n$ un multiple positif du i-ème vecteur de base de $\mathbb R^n$, les vecteurs restants ayant pour image 0. Dans ces bases, l'application $T$ est donc représentée par une matrice diagonale dont les coefficients sont des réels positifs.
***Exo:*** Quel est le lien entre le texte précédent et la formule: ***Exo:*** Quel est le lien entre le texte précédent et la formule:
U,S,V = np.linalg.svd(A) U,S,V = np.linalg.svd(A)
Indiquez précisément quelles sont les bases. Attention, il y a un piège. Indiquez précisément quelles sont les bases. Attention, il y a un piège.
Observons cela avec une application $T : \mathbb R^2 \to \mathbb R^2, x \to ax$. Observons cela avec une application $T : \mathbb R^2 \to \mathbb R^2, x \to ax$.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import matplotlib.patches as patches import matplotlib.patches as patches
def showBasis(vecs,lim=(-2,2),color0=(1,0,0),color1=(0,1,0): def showBasis(vecs,lim=(-2,2),color0=(1,0,0),color1=(0,1,0)):
vec0=vecs[:,0] vec0=vecs[:,0]
vec1=vecs[:,1] vec1=vecs[:,1]
plt.arrow(0,0,vec0[0],vec0[1], head_width=0.1, head_length=0.2,color=color0) plt.arrow(0,0,vec0[0],vec0[1], head_width=0.1, head_length=0.2,color=color0)
plt.arrow(0,0,vec1[0],vec1[1], head_width=0.1, head_length=0.2,color=color1) plt.arrow(0,0,vec1[0],vec1[1], head_width=0.1, head_length=0.2,color=color1)
# plus simplement, on peut faire # plus simplement, on peut faire
# plt.plot([0,vec0[0]],[0,vec1[0]],color=[1,0,blue]) # plt.plot([0,vec0[0]],[0,vec1[0]],color=[1,0,blue])
plt.gca().set_aspect("equal") plt.gca().set_aspect("equal")
# attention avec plt.arrow il faut nécessaire préciser le domaine (sinon c'est moche) # attention avec plt.arrow il faut nécessaire préciser le domaine (sinon c'est moche)
plt.xlim(lim[0],lim[1]) plt.xlim(lim[0],lim[1])
plt.ylim(lim[0],lim[1]) plt.ylim(lim[0],lim[1])
"""matrice de l'application""" """matrice de l'application"""
a=np.array([[1,-1],[-0.7,0.1]]) a=np.array([[1,-1],[-0.7,0.1]])
"""observons comment elle transforme la base canonique """ """observons comment elle transforme la base canonique """
vecs=np.eye(2) vecs=np.eye(2)
plt.subplot(1,2,1) plt.subplot(1,2,1)
showBasis(vecs) showBasis(vecs)
plt.subplot(1,2,2) plt.subplot(1,2,2)
showBasis(a@vecs) showBasis(a@vecs)
``` ```
%% Output %% Output
File "<ipython-input-23-0692074bb260>", line 2
def showBasis(vecs,lim=(-2,2),color0=(1,0,0),color1=(0,1,0):
^
SyntaxError: invalid syntax
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
""" observons maintenant les transformations pas à pas """ """ observons maintenant les transformations pas à pas """
U,S,V = np.linalg.svd(a) U,S,V = np.linalg.svd(a)
print("a\n",a) print("a\n",a)
print("\nU\n",U) print("\nU\n",U)
print("\nS\n",S) print("\nS\n",S)
print("\nV\n",V) print("\nV\n",V)
plt.figure(figsize=(16,4)) plt.figure(figsize=(16,4))
plt.subplot(1,4,1) plt.subplot(1,4,1)
showBasis(vecs) showBasis(vecs)
plt.subplot(1,4,2) plt.subplot(1,4,2)
showBasis(np.eye(2),color0=(0,0,1),color1=(0,0,1)) showBasis(np.eye(2),color0=(0,0,1),color1=(0,0,1))
showBasis(V@vecs) showBasis(V@vecs)
plt.subplot(1,4,3) plt.subplot(1,4,3)
showBasis(np.diag(S),color0=(0,0,1),color1=(0,0,1)) showBasis(np.diag(S),color0=(0,0,1),color1=(0,0,1))
showBasis(np.diag(S)@V@vecs) showBasis(np.diag(S)@V@vecs)
plt.subplot(1,4,4) plt.subplot(1,4,4)
showBasis(U@np.diag(S)@V@vecs) showBasis(U@np.diag(S)@V@vecs)
``` ```
%% Output %% Output
a a
[[ 1. -1. ] [[ 1. -1. ]
[-0.7 0.1]] [-0.7 0.1]]
U U
[[-0.92 0.4 ] [[-0.92 0.4 ]
[ 0.4 0.92]] [ 0.4 0.92]]
S S
[1.53 0.39] [1.53 0.39]
V V
[[-0.78 0.62] [[-0.78 0.62]
[-0.62 -0.78]] [-0.62 -0.78]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***A vous:*** Commentez: à quel moment fait-on une isométrie (laquelle?), une dilatation ? ***A vous:*** Commentez: à quel moment fait-on une isométrie (laquelle?), une dilatation ?
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Cas symétrique défini positif ### Cas symétrique défini positif
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
""" prenons une matrice symétrique, diagonalisons là""" """ prenons une matrice symétrique, diagonalisons là"""
B=np.array([[2,-1],[-1,2]]) B=np.array([[2,-1],[-1,2]])
val_pr,vec_pr=np.linalg.eig(B) val_pr,vec_pr=np.linalg.eig(B)
print("val_pr:\n",val_pr) print("val_pr:\n",val_pr)
print("vec_pr:\n",vec_pr) print("vec_pr:\n",vec_pr)
``` ```
%% Output %% Output
val_pr: val_pr:
[3. 1.] [3. 1.]
vec_pr: vec_pr:
[[ 0.71 0.71] [[ 0.71 0.71]
[-0.71 0.71]] [-0.71 0.71]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Quel est le lien entre la diagonalisation et la svd dans ce cas? Quel est le lien entre la diagonalisation et la svd dans ce cas?
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Vecteurs aléatoires ## Vecteurs aléatoires
Un vecteur aléatoire $X \in \mathbb R^p$, c'est simplement une collection de $p$-variables aléatoires. Mathématiquement, on les note en colonne: Un vecteur aléatoire $X \in \mathbb R^p$, c'est simplement une collection de $p$-variables aléatoires. Mathématiquement, on les note en colonne:
$$ $$
X= \begin{bmatrix} X= \begin{bmatrix}
X_0\\ X_0\\
X_1\\ X_1\\
\vdots\\ \vdots\\
X_{p-1} X_{p-1}
\end{bmatrix} \end{bmatrix}
$$ $$
Mais informatiquement, l'habitude est plutôt de les noter en ligne! Notamment si l'on dispose de plusieurs réalisations $X^{(0)},X^{(1)},...,X^{(i)}$ de $X$, on les disposera ainsi: Mais informatiquement, l'habitude est plutôt de les noter en ligne! Notamment si l'on dispose de plusieurs réalisations $X^{(0)},X^{(1)},...,X^{(i)}$ de $X$, on les disposera ainsi:
$$ $$
\begin{bmatrix} \begin{bmatrix}
X^{(0)} \\ X^{(0)} \\
X^{(1)}\\ X^{(1)}\\
\vdots\\ \vdots\\
X^{(i)}\\ X^{(i)}\\
\vdots \vdots
\end{bmatrix}= \end{bmatrix}=
\begin{bmatrix} \begin{bmatrix}
X^{(0)}_0 & \dots & X^{(0)}_{p-1}\\ X^{(0)}_0 & \dots & X^{(0)}_{p-1}\\
X^{(1)}_0 & \dots & X^{(1)}_{p-1}\\ X^{(1)}_0 & \dots & X^{(1)}_{p-1}\\
\vdots\\ \vdots\\
X^{(i)}_0 & \dots & X^{(i)}_{p-1}\\ X^{(i)}_0 & \dots & X^{(i)}_{p-1}\\
\vdots \vdots
\end{bmatrix} \end{bmatrix}
$$ $$
Cette disposition des données s'appelle une ***dataFrame***. Imaginons que le vecteur aléatoire représente les caractérisitiques d'un individu (poids, taille, QI,...). Alors: Cette disposition des données s'appelle une ***dataFrame***. Imaginons que le vecteur aléatoire représente les caractérisitiques d'un individu (poids, taille, QI,...). Alors:
* chaque ligne de la dataFrame représente un individu * chaque ligne de la dataFrame représente un individu
* chaque colonne représente une caractéristique * chaque colonne représente une caractéristique
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Observons des réalisations d'un vecteur aléatoire `X=(X0,X1)` de dimension 2. Observons des réalisations d'un vecteur aléatoire `X=(X0,X1)` de dimension 2.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
nbData=500 nbData=500
X0=np.random.uniform(low=0,high=10,size=nbData) X0=np.random.uniform(low=0,high=10,size=nbData)
X1=np.random.poisson(lam=X0,size=nbData) X1=np.random.poisson(lam=X0,size=nbData)
""" on crée la dataFrame""" """ on crée la dataFrame"""
X=np.stack([X0,X1],axis=1) X=np.stack([X0,X1],axis=1)
plt.plot(X[:,0],X[:,1],'.') plt.plot(X[:,0],X[:,1],'.')
plt.gca().set_aspect('equal') plt.gca().set_aspect('equal')
``` ```
%% Output %% Output
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***Exo:*** Les deux composantes `X0` et `X1` construites dans le code sont-elle indépendantes? Justifiez votre réponse en observant la fonction `x -> Loi(X1/X0=x)` ***Exo:*** Les deux composantes `X0` et `X1` construites dans le code sont-elle indépendantes? Justifiez votre réponse en observant la fonction `x -> Loi(X1/X0=x)`
***Exo:*** `X0` admet-t-il une densité par rapport à la mesure de Lebesgue? Par rapport à la mesure de comptage? Même quesion pour `X1`, même question pour le couple `X= (X0,X1)`. ***Exo:*** `X0` admet-t-il une densité par rapport à la mesure de Lebesgue? Par rapport à la mesure de comptage? Même quesion pour `X1`, même question pour le couple `X= (X0,X1)`.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Espérance, Covariance, Corrélation ## Espérance, Covariance, Corrélation
### Définissions ### Définissions
* L'espérance de $X$ c'est le vecteur $\mu = \mathbf E[X]$ défini par $\mu_i = \mathbf E[X_i]$. * L'espérance de $X$ c'est le vecteur $\mu = \mathbf E[X]$ défini par $\mu_i = \mathbf E[X_i]$.
* la matrice de covariance de $X$ c'est la matrice $\Sigma^2$ définie par $\Sigma^2_{ij}= \mathrm{cov}(X_i,X_j)$. On peut l'écrire avec une multiplication matricielle: * la matrice de covariance de $X$ c'est la matrice $\Sigma^2$ définie par $\Sigma^2_{ij}= \mathrm{cov}(X_i,X_j)$. On peut l'écrire avec une multiplication matricielle:
$$ $$
\Sigma^2 = \mathbf E[ (X-\mu) (X-\mu)^T ] \Sigma^2 = \mathbf E[ (X-\mu) (X-\mu)^T ]
$$ $$
(où nous supposons, comme pour les vecteurs, que l'espérance d'une matrice aléatoire, c'est la matrice des espérances). (où nous supposons, comme pour les vecteurs, que l'espérance d'une matrice aléatoire, c'est la matrice des espérances).
* Les coefficients de corrélation sont définis par * Les coefficients de corrélation sont définis par
$$ $$
c_{ij} = \frac{\mathrm{cov}(X_i,X_j) }{ \sqrt{\mathbf V(X_i)} \sqrt{\mathbf V(X_j)} } c_{ij} = \frac{\mathrm{cov}(X_i,X_j) }{ \sqrt{\mathbf V(X_i)} \sqrt{\mathbf V(X_j)} }
$$ $$
Quelle fameuse inégalité permet d'affirmer que ces coefficients sont compris entre $-1$ et $+1$? Quelle fameuse inégalité permet d'affirmer que ces coefficients sont compris entre $-1$ et $+1$?
***Exo:*** Soit $a$ une matrice et $b$ un vecteur. On a: ***Exo:*** Soit $a$ une matrice et $b$ un vecteur. On a:
* $\mathbf E[aX+b] =a \mathbf E[X]+b$. * $\mathbf E[aX+b] =a \mathbf E[X]+b$.
* $\mathbf V[aX+b] = a\mathbf V[X]a^T$. * $\mathbf V[aX+b] = a\mathbf V[X]a^T$.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
A partir d'une dataFrame, on peut estimer l'espérance et la matrice de covariance du vecteur aléatoire sous-jacent: A partir d'une dataFrame, on peut estimer l'espérance et la matrice de covariance du vecteur aléatoire sous-jacent:
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
mean=np.mean(X,axis=0) mean=np.mean(X,axis=0)
print("mean:",mean) print("mean:",mean)
X_cen = X - mean X_cen = X - mean
cov = X_cen.T @ X_cen / nbData cov = X_cen.T @ X_cen / nbData
print("covariance à la main:\n",cov) print("covariance à la main:\n",cov)
cov=np.cov(X.T) cov=np.cov(X.T)
print("covariance toute faite:\n",cov) print("covariance toute faite:\n",cov)
``` ```
%% Output %% Output
mean: [5.05 4.88] mean: [5.05 4.88]
covariance à la main: covariance à la main:
[[ 8.33 7.88] [[ 8.33 7.88]
[ 7.88 11.9 ]] [ 7.88 11.9 ]]
covariance toute faites: covariance toute faites:
[[ 8.35 7.89] [[ 8.35 7.89]
[ 7.89 11.92]] [ 7.89 11.92]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Insistons sur le fait que nous n'avons fait que des estimations de l'espérance et la matrice de covariance. D'ailleurs l'estimateur de la covariance serait un peu meilleur si... (Aide: en numpy, l'amélioration se fait avec `np.cov(X.T,ddof=1)`) Insistons sur le fait que nous n'avons fait que des estimations de l'espérance et la matrice de covariance. D'ailleurs l'estimateur de la covariance serait un peu meilleur si... (Aide: en numpy, l'amélioration se fait avec `np.cov(X.T,ddof=1)`)
***Exo:*** Sauriez-vous écrire une formule explicite qui donne les vraies espérances des v.a `X0` et `X1`? Et les variances/covariance? Attention, ce n'est pas facile. Certaines formules comporteront des symboles $\int$ et $\sum$ mélangés. ***Exo:*** Sauriez-vous écrire une formule explicite qui donne les vraies espérances des v.a `X0` et `X1`? Et les variances/covariance? Attention, ce n'est pas facile. Certaines formules comporteront des symboles $\int$ et $\sum$ mélangés.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Interprétation géométrique de la covariance ### Interprétation géométrique de la covariance
La matrice $\Sigma^2$ est symétrique et ses valeurs propres sont positives. Donc on peut l'écrire La matrice $\Sigma^2$ est symétrique et ses valeurs propres sont positives. Donc on peut l'écrire
$$ $$
\Sigma^2=U S^2 U^T \Sigma^2=U S^2 U^T
$$ $$
avec $S^2$ la matrice diagonale formée des valeurs singulières (=valeurs propres) que l'on classe dans l'ordre décroissant $s_{0}^2>s_{1}^2 >... $, et $U$ une matrice orthogonale. avec $S^2$ la matrice diagonale formée des valeurs singulières (=valeurs propres) que l'on classe dans l'ordre décroissant $s_{0}^2>s_{1}^2 >... $, et $U$ une matrice orthogonale.
Considérons un vecteur aléatoire $X$ de matrice de covariance $\Sigma^2= U S^2 U^T$. Notons $U_{i}$ les colonnes de $U$. On note naturellement $s_i$ les racines carrées des $s^2_i=S^2_{i,i}$. Elles sont ordonnées $s_0>s_1 >... $. On note $\mu$ le vecteur espérance de $X$. Si maintenant nous simulons des copies indépendantes de $X$, elles formeront un nuage de point autour de $\mu$, dont la dispersion sera décrite par les $U_i$ et $s_i$ : Considérons un vecteur aléatoire $X$ de matrice de covariance $\Sigma^2= U S^2 U^T$. Notons $U_{i}$ les colonnes de $U$. On note naturellement $s_i$ les racines carrées des $s^2_i=S^2_{i,i}$. Elles sont ordonnées $s_0>s_1 >... $. On note $\mu$ le vecteur espérance de $X$. Si maintenant nous simulons des copies indépendantes de $X$, elles formeront un nuage de point autour de $\mu$, dont la dispersion sera décrite par les $U_i$ et $s_i$ :
***Exo:*** ***Exo:***
* Soit $X$ un vecteur aléatoire de matrice de covariance $\Sigma^2=U S^2 U^T$ et d'espérance $\mu$. Notons $\Sigma = U S U^T.$ Quelles est la matrice de covariance de $\Sigma^{-1} (X-\mu) $ ? Faites le lien avec le fait de centrer-réduire les va. * Soit $X$ un vecteur aléatoire de matrice de covariance $\Sigma^2=U S^2 U^T$ et d'espérance $\mu$. Notons $\Sigma = U S U^T.$ Quelles est la matrice de covariance de $\Sigma^{-1} (X-\mu) $ ? Faites le lien avec le fait de centrer-réduire les va.
* $\Sigma^2$ est symétrique et ses valeurs propres sont positives. La symétrie est évidente. Pour les valeurs propres, il suffit de montrer qu'elle est semi-définie positive. Vérifions-le : $\forall v : v^T \Sigma^2 v = ... \geq 0$. * $\Sigma^2$ est symétrique et ses valeurs propres sont positives. La symétrie est évidente. Pour les valeurs propres, il suffit de montrer qu'elle est semi-définie positive. Vérifions-le : $\forall v : v^T \Sigma^2 v = ... \geq 0$.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
U,S2,_=np.linalg.svd(cov) U,S2,_=np.linalg.svd(cov)
print(U@U.T) print(U@U.T)
S=np.sqrt(S2) S=np.sqrt(S2)
plt.plot(X_cen[:,0],X_cen[:,1],'.'); plt.plot(X_cen[:,0],X_cen[:,1],'.');
showBasis(U@np.diag(S),lim=(-10,10)) showBasis(U@np.diag(S),lim=(-10,10))
``` ```
%% Output %% Output
[[1. 0.] [[1. 0.]
[0. 1.]] [0. 1.]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***Attention 1:*** même en classant les valeurs propres par ordre décroissant, il n'y a pas unicité de la base $U$ (il y a un choix à faire quand des valeurs propres sont égales). ***Attention 1:*** même en classant les valeurs propres par ordre décroissant, il n'y a pas unicité de la base $U$ (il y a un choix à faire quand des valeurs propres sont égales).
***Attention 2 :*** tous les vecteurs aléatoires ne se répartissent pas en patate. Voici des simulations de vecteurs aléatoires de $\mathbb R^2$ qui admettent tous la matrice identité pour matrice de covariance (merci wikipedia). ***Attention 2 :*** tous les vecteurs aléatoires ne se répartissent pas en patate. Voici des simulations de vecteurs aléatoires de $\mathbb R^2$ qui admettent tous la matrice identité pour matrice de covariance (merci wikipedia).
![nuages de point avec covariance identité](img/identite.png) ![nuages de point avec covariance identité](img/identite.png)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Vecteurs aléatoires gaussiens ## Vecteurs aléatoires gaussiens
Un vecteur aléatoire est gaussien, si toutes les combinaisons linéaires qu'on peut en faire sont des v.a gaussiennes. Un vecteur aléatoire est gaussien, si toutes les combinaisons linéaires qu'on peut en faire sont des v.a gaussiennes.
Je trouve cette définition classique très peu constructive. En voici une plus explicite, en deux étapes. Je trouve cette définition classique très peu constructive. En voici une plus explicite, en deux étapes.
### gaussien standart $\to$ gaussien général ### gaussien standart $\to$ gaussien général
* Soit $Y\in \mathbb R^p$. * Soit $Y\in \mathbb R^p$.
On dit que $Y$ est un vecteur gaussien standart, et on note $Y\sim \mathcal N_p(0,I) $ lorsque les composantes $Y_i$ sont des $\mathcal N(0,1)$ indépendantes. Son espérance est $0=(0,...,0)$ et sa matrice de covariance c'est $I$. On dit que $Y$ est un vecteur gaussien standart, et on note $Y\sim \mathcal N_p(0,I) $ lorsque les composantes $Y_i$ sont des $\mathcal N(0,1)$ indépendantes. Son espérance est $0=(0,...,0)$ et sa matrice de covariance c'est $I$.
* Les vecteurs gaussiens sont tous les vecteurs de la forme $a Y+b$ avec $a$ matrice, $b$ vecteur et $Y\sim\mathcal N_p(0,I)$. * Les vecteurs gaussiens sont tous les vecteurs de la forme $a Y+b$ avec $a$ matrice, $b$ vecteur et $Y\sim\mathcal N_p(0,I)$.
Ainsi par construction, la famille des vecteurs gaussiens est stable par combinaison affine. Ainsi par construction, la famille des vecteurs gaussiens est stable par combinaison affine.
### gaussien général $\to$ gaussien standart ### gaussien général $\to$ gaussien standart
Ainsi, par définition, un vecteur gaussien général c'est $X= aY + b$ avec $Y \sim \mathcal N_p(0,I)$. L'espérance de $X$ est alors $\mu = b$ et sa matrice de covariance est alors $\Sigma^2 = a a^T$ (cf. ancien exo). Ainsi, par définition, un vecteur gaussien général c'est $X= aY + b$ avec $Y \sim \mathcal N_p(0,I)$. L'espérance de $X$ est alors $\mu = b$ et sa matrice de covariance est alors $\Sigma^2 = a a^T$ (cf. ancien exo).
Définissons alors $Y' := \Sigma^{-1} (X-\mu)$ (attention, on ne retombe pas forcément sur $Y$). On a que Définissons alors $Y' := \Sigma^{-1} (X-\mu)$ (attention, on ne retombe pas forcément sur $Y$). On a que
$$ $$
X = \Sigma Y' + \mu \qquad \mathrm{ et } \qquad Y' \sim N_p(0,I) X = \Sigma Y' + \mu \qquad \mathrm{ et } \qquad Y' \sim N_p(0,I)
$$ $$
Ci-dessus, on a une écriture "canonique" d'un vecteur Gaussien général (une écriture pas tout à fait unique à cause des choix arbitraires dans la SVD). On voit bien avec cette écriture que la loi de $X$ ne dépend que de sa matrice de covariance $\Sigma^2$ et de son espérance $\mu$. On notera d'ailleurs: Ci-dessus, on a une écriture "canonique" d'un vecteur Gaussien général (une écriture pas tout à fait unique à cause des choix arbitraires dans la SVD). On voit bien avec cette écriture que la loi de $X$ ne dépend que de sa matrice de covariance $\Sigma^2$ et de son espérance $\mu$. On notera d'ailleurs:
$$ $$
X \sim \mathcal N_p(\mu, \Sigma^2) X \sim \mathcal N_p(\mu, \Sigma^2)
$$ $$
Voyons comment on simule un tel vecteur avec `numpy` Voyons comment on simule un tel vecteur avec `numpy`
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
""" choisissons une matrice de covariance """ """ choisissons une matrice de covariance """
Sigma=np.array([[2,-1],[-1,1]]) Sigma=np.array([[2,-1],[-1,1]])
print("cov\n",cov) print("cov\n",cov)
""" un vecteur espérance""" """ un vecteur espérance"""
mu=np.array([3,4]) mu=np.array([3,4])
""" simulons""" """ simulons"""
X= np.random.multivariate_normal(mean=mu,cov=Sigma,size=500) X= np.random.multivariate_normal(mean=mu,cov=Sigma,size=500)
plt.plot(X[:,0],X[:,1],'.'); plt.plot(X[:,0],X[:,1],'.');
``` ```
%% Output %% Output
cov cov
[[ 8.35 7.89] [[ 8.35 7.89]
[ 7.89 11.92]] [ 7.89 11.92]]
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***A vous:*** ***A vous:***
* Pourquoi le nuage de points est-il penché dans ce sens? (et pas selon la première bisectrice par exemple) * Pourquoi le nuage de points est-il penché dans ce sens? (et pas selon la première bisectrice par exemple)
* Que se passe-t-il si l'on prend `cov` une matrice dégénérée? (faites le!) * Que se passe-t-il si l'on prend `cov` une matrice dégénérée? (faites le!)
* en utilisant la théorie, simulez vous-même ce nuage de point sans utiliser `np.random.multivariate_normal` mais en partant simplement de v.a gaussiennes via `np.random.normal`. Aide: il vous faut aussi utiliser la SVD (et relire le paragraphe précédent). * en utilisant la théorie, simulez vous-même ce nuage de point sans utiliser `np.random.multivariate_normal` mais en partant simplement de v.a gaussiennes via `np.random.normal`. Aide: il vous faut aussi utiliser la SVD (et relire le paragraphe précédent).
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
## Densité de vecteur aléatoire ## Densité de vecteur aléatoire
Notre tout premier exemple de vecteur aléatoire n'admettait pas de densité par rapport à la mesure de Lebesgue, ni par rapport à la mesure de comptage. Mais dans ce nombreux cas pratiques, une telle densité existe. Nous traitons ici de la densité par rapport à la mesure de Lebesgue. Notre tout premier exemple de vecteur aléatoire n'admettait pas de densité par rapport à la mesure de Lebesgue, ni par rapport à la mesure de comptage. Mais dans ce nombreux cas pratiques, une telle densité existe. Nous traitons ici de la densité par rapport à la mesure de Lebesgue.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
""" une dataFrame à deux colonnes """ """ une dataFrame à deux colonnes """
def betaBeta(nbData): def betaBeta(nbData):
X0=np.random.beta(a=1,b=3,size=nbData) X0=np.random.beta(a=1,b=3,size=nbData)
X1=np.random.beta(a=2,b=1,size=nbData) X1=np.random.beta(a=2,b=1,size=nbData)
return np.stack([X0,X1],axis=1) return np.stack([X0,X1],axis=1)
""" la densité correspondante """ """ la densité correspondante """
def densityX(x,y): def densityX(x,y):
return stats.beta.pdf(x,a=1,b=3)*stats.beta.pdf(y,a=2,b=1) return stats.beta.pdf(x,a=1,b=3)*stats.beta.pdf(y,a=2,b=1)
lim=1 lim=1
xx=np.linspace(0, lim, 100) xx=np.linspace(0, lim, 100)
yy=np.linspace(0, lim, 100) yy=np.linspace(0, lim, 100)
XX,YY=np.meshgrid(xx,yy) XX,YY=np.meshgrid(xx,yy)
den_square=densityX(XX,YY) den_square=densityX(XX,YY)
plt.imshow(den_square,extent=[0,lim,0,lim],origin="lower",cmap="jet") plt.imshow(den_square,extent=[0,lim,0,lim],origin="lower",cmap="jet")
X=betaBeta(500) X=betaBeta(500)
plt.plot(X[:,0],X[:,1],'.') plt.plot(X[:,0],X[:,1],'.')
plt.gca().set_aspect('equal') plt.gca().set_aspect('equal')
``` ```
%% Output %% Output
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***Exo théorique:*** Les deux composantes `X0` et `X1` sont indépendantes (car on appelle deux fois le générateur aléatoire). Comment est-ce que cela se lit sur la densité? ***Exo théorique:*** Les deux composantes `X0` et `X1` sont indépendantes (car on appelle deux fois le générateur aléatoire). Comment est-ce que cela se lit sur la densité?
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***Exo informatique:*** Refaites un programme similaire choisissant certains paramètres de la loi beta <1. La difficulté c'est qu'on obtient des densités non bornées. Pour effectuer la multiplicaition `f0* f1` il faut alors supprimer les valeurs infinies. Pour cela vous pouvez ***Exo informatique:*** Refaites un programme similaire choisissant certains paramètres de la loi beta <1. La difficulté c'est qu'on obtient des densités non bornées. Pour effectuer la multiplicaition `f0* f1` il faut alors supprimer les valeurs infinies. Pour cela vous pouvez
* soit restreindre le domaine en excluant les bords du carré $[0,1]^2$ * soit restreindre le domaine en excluant les bords du carré $[0,1]^2$
* soit limiter les densités, ex: `f0=np.minimum(f0,1000)` * soit limiter les densités, ex: `f0=np.minimum(f0,1000)`
Mais même en prenant ces précautions, les densités seront moches car les couleurs seront saturées. Astuce: tracez le log de la densité (ce qui revient à changer l'échelle de couleur) Mais même en prenant ces précautions, les densités seront moches car les couleurs seront saturées. Astuce: tracez le log de la densité (ce qui revient à changer l'échelle de couleur)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Cas Gaussien ### Cas Gaussien
La densité des vecteurs gaussiens est fournie dans `scipy.stats` La densité des vecteurs gaussiens est fournie dans `scipy.stats`
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
cov=np.array([[2,-1],[-1,1]]) cov=np.array([[2,-1],[-1,1]])
print("cov\n",cov) print("cov\n",cov)
X= np.random.multivariate_normal(mean=[0,0],cov=cov,size=500) X= np.random.multivariate_normal(mean=[0,0],cov=cov,size=500)
print("X.shape:",X.shape) print("X.shape:",X.shape)
plt.plot(X[:,0],X[:,1],'.') plt.plot(X[:,0],X[:,1],'.')
plt.gca().set_aspect("equal") plt.gca().set_aspect("equal")
lim=3 lim=3
resolution=200 resolution=200
xx=np.linspace(-lim,lim,resolution) xx=np.linspace(-lim,lim,resolution)
XX,YY=np.meshgrid(xx,xx) XX,YY=np.meshgrid(xx,xx)
XY=np.stack([XX,YY],axis=2) XY=np.stack([XX,YY],axis=2)
den=stats.multivariate_normal.pdf(XY,mean=[0,0],cov=cov) den=stats.multivariate_normal.pdf(XY,mean=[0,0],cov=cov)
plt.imshow(den,cmap="jet",extent=[-lim,lim,-lim,lim],origin="lower"); plt.imshow(den,cmap="jet",extent=[-lim,lim,-lim,lim],origin="lower");
``` ```
%% Output %% Output
cov cov
[[ 2 -1] [[ 2 -1]
[-1 1]] [-1 1]]
X.shape: (500, 2) X.shape: (500, 2)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
### Déformation par multiplication matricielle ### Déformation par multiplication matricielle
Soit $X\in \mathbb R^p$ un vecteur aléatoire de densité $f$. Soit $a$ une matrice inversible $p\times p$. Soit $b\in \mathtt R^p$. La densité de $aX+b$ est: Soit $X\in \mathbb R^p$ un vecteur aléatoire de densité $f$. Soit $a$ une matrice inversible $p\times p$. Soit $b\in \mathtt R^p$. La densité de $aX+b$ est:
$$ $$
x\to \frac 1 {| \mathrm{det} \,a |} f \Big( a^{-1} (x-b) \Big ) x\to \frac 1 {| \mathrm{det} \,a |} f \Big( a^{-1} (x-b) \Big )
$$ $$
***Exo:*** Vérifiez cette formule en utilisant la technique de la fonction test. ***Exo:*** Vérifiez cette formule en utilisant la technique de la fonction test.
$$ $$
\mathbf E[\phi(aX+b)] = \int_{\mathbb R^p} \phi(aX+b) \ f(x) \ dx = ... \mathbf E[\phi(aX+b)] = \int_{\mathbb R^p} \phi(aX+b) \ f(x) \ dx = ...
$$ $$
***Exo:*** Retrouvez la densité de $\mathcal N_p(0,\Sigma^2)$ à partir de la densité de $\mathcal N_p(0,I)$. Aide: Avec la svd on sait qu'on peut écrire $\Sigma^2=\Sigma^T \Sigma$. Mais cette racine carrée $\Sigma$(quelque peu arbitraire) ne doit pas apparaître dans l'expression finale. ***Exo:*** Retrouvez la densité de $\mathcal N_p(0,\Sigma^2)$ à partir de la densité de $\mathcal N_p(0,I)$. Aide: Avec la svd on sait qu'on peut écrire $\Sigma^2=\Sigma^T \Sigma$. Mais cette racine carrée $\Sigma$(quelque peu arbitraire) ne doit pas apparaître dans l'expression finale.
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Déformons les réalisations des vecteurs aléatoires `betaBeta`. Déformons les réalisations des vecteurs aléatoires `betaBeta`.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
X= betaBeta(500) X= betaBeta(500)
a=np.array([[1,-1],[-0.7,0.1]]) a=np.array([[1,-1],[-0.7,0.1]])
print("X.shape:",X.shape) print("X.shape:",X.shape)
""" attention: informatiquement les vecteurs aléatoires sont en ligne """ attention: informatiquement les vecteurs aléatoires sont en ligne
et non en colonne comme dans la théorie. et non en colonne comme dans la théorie.
Du coup, la multiplication par la matrice se fait à droite: Du coup, la multiplication par la matrice se fait à droite:
""" """
aX=X@a aX=X@a
plt.plot(aX[:,0],aX[:,1],'.') plt.plot(aX[:,0],aX[:,1],'.')
plt.gca().set_aspect('equal') plt.gca().set_aspect('equal')
``` ```
%% Output %% Output
X.shape: (500, 2) X.shape: (500, 2)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Pour tracer informatiquement la densité déformée il faut bien réfléchir. Pour tracer informatiquement la densité déformée il faut bien réfléchir.
Observons ce que renvoie `meshgrid` sur un petit exemple. Voyons comment nous devons le transformer. Observons ce que renvoie `meshgrid` sur un petit exemple. Voyons comment nous devons le transformer.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
xx=[1,2,3] xx=[1,2,3]
yy=[0,1,2,3] yy=[0,1,2,3]
""" meshgrid répète les abscisses et les ordonnées """ """ meshgrid répète les abscisses et les ordonnées """
XX,YY=np.meshgrid(xx,yy) XX,YY=np.meshgrid(xx,yy)
print("XX\n",XX) print("XX\n",XX)
print("YY\n",YY) print("YY\n",YY)
""" on colle le tout""" """ on colle le tout"""
XY=np.stack([XX,YY]) XY=np.stack([XX,YY])
print("XY.shape:",XY.shape) print("XY.shape:",XY.shape)
""" on transforme les deux matrices en une liste de couple (abscisse,ordonnée)""" """ on transforme les deux matrices en une liste de couple (abscisse,ordonnée)"""
XY_flat=np.reshape(XY,newshape=[2,4*3]) XY_flat=np.reshape(XY,newshape=[2,4*3])
print("XY_flat.shape:",XY_flat.shape) print("XY_flat.shape:",XY_flat.shape)
print("XY_flat\n",XY_flat) print("XY_flat\n",XY_flat)
"""sur chacun de ces couples, on effectue la multipliation matricielle""" """sur chacun de ces couples, on effectue la multipliation matricielle"""
print("a@XY_flat\n",a@XY_flat) print("a@XY_flat\n",a@XY_flat)
``` ```
%% Output %% Output
XX XX
[[1 2 3] [[1 2 3]
[1 2 3] [1 2 3]
[1 2 3] [1 2 3]
[1 2 3]] [1 2 3]]
YY YY
[[0 0 0] [[0 0 0]
[1 1 1] [1 1 1]
[2 2 2] [2 2 2]
[3 3 3]] [3 3 3]]
XY.shape: (2, 4, 3) XY.shape: (2, 4, 3)
XY_flat.shape: (2, 12) XY_flat.shape: (2, 12)
XY_flat XY_flat
[[1 2 3 1 2 3 1 2 3 1 2 3] [[1 2 3 1 2 3 1 2 3 1 2 3]
[0 0 0 1 1 1 2 2 2 3 3 3]] [0 0 0 1 1 1 2 2 2 3 3 3]]
a@XY_flat a@XY_flat
[[ 1. 2. 3. 0. 1. 2. -1. 0. 1. -2. -1. 0. ] [[ 1. 2. 3. 0. 1. 2. -1. 0. 1. -2. -1. 0. ]
[-0.7 -1.4 -2.1 -0.6 -1.3 -2. -0.5 -1.2 -1.9 -0.4 -1.1 -1.8]] [-0.7 -1.4 -2.1 -0.6 -1.3 -2. -0.5 -1.2 -1.9 -0.4 -1.1 -1.8]]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
def deform_density(density, x, y,a): def deform_density(density, x, y,a):
a_inv=np.linalg.inv(a) a_inv=np.linalg.inv(a)
xy=np.stack([x,y]) xy=np.stack([x,y])
axy=a_inv@xy axy=a_inv@xy
ax=axy[0] ax=axy[0]
ay=axy[1] ay=axy[1]
return density(ax,ay) return density(ax,ay)
lim=1 lim=1
resolution=200 resolution=200
xx=np.linspace(-lim,lim,resolution) xx=np.linspace(-lim,lim,resolution)
XX,YY=np.meshgrid(xx,xx) XX,YY=np.meshgrid(xx,xx)
XY=np.stack([XX,YY]) XY=np.stack([XX,YY])
print("XY.shape:",XY.shape) print("XY.shape:",XY.shape)
XY_flat=np.reshape(XY,newshape=[2,resolution*resolution]) XY_flat=np.reshape(XY,newshape=[2,resolution*resolution])
print("XY_flat.shape:",XY_flat.shape) print("XY_flat.shape:",XY_flat.shape)
den_flat=deform_density(densityX, XY_flat[0], XY_flat[1],a) den_flat=deform_density(densityX, XY_flat[0], XY_flat[1],a)
print("den_flat.shape:",den_flat.shape) print("den_flat.shape:",den_flat.shape)
den_square=np.reshape(den_flat,newshape=[resolution,resolution]) den_square=np.reshape(den_flat,newshape=[resolution,resolution])
print("den_square.shape:",den_square.shape) print("den_square.shape:",den_square.shape)
plt.imshow(den_square,cmap="jet",origin="lower",extent=[-lim,lim,-lim,lim]) plt.imshow(den_square,cmap="jet",origin="lower",extent=[-lim,lim,-lim,lim])
plt.plot(aX[:,0],aX[:,1],'.') plt.plot(aX[:,0],aX[:,1],'.')
plt.gca().set_aspect('equal'); plt.gca().set_aspect('equal');
``` ```
%% Output %% Output
XY.shape: (2, 200, 200) XY.shape: (2, 200, 200)
XY_flat.shape: (2, 40000) XY_flat.shape: (2, 40000)
den_flat.shape: (40000,) den_flat.shape: (40000,)
den_square.shape: (200, 200) den_square.shape: (200, 200)
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
***A vous:*** Mince, cela ne colle pas tout à fait. Il y a un bug dans le programme précédent. Aide: la densité est la bonne (à une constante près), ce sont les simulations qui buguent. ***A vous:*** Mince, cela ne colle pas tout à fait. Il y a un bug dans le programme précédent. Aide: la densité est la bonne (à une constante près), ce sont les simulations qui buguent.
***Exo:*** Ajouter les deux directions principales de ce nuage de point, en utilisant `showBasis`. Il y a plusieurs techniques pour les calculer: ***Exo:*** Ajouter les deux directions principales de ce nuage de point, en utilisant `showBasis`. Il y a plusieurs techniques pour les calculer:
* estimer la matrice de covariance avec `np.cov` * estimer la matrice de covariance avec `np.cov`
* calculer exactement cette matrice (ici c'est très facile, il suffit de regarder la formule pour la variance d'une loi beta, que l'on peut aussi avoir avec `stats.beta.std`) * calculer exactement cette matrice (ici c'est très facile, il suffit de regarder la formule pour la variance d'une loi beta, que l'on peut aussi avoir avec `stats.beta.std`)
* ou partir directement de la svd ou de la diagonalisation de matrice `a`. * ou partir directement de la svd ou de la diagonalisation de matrice `a`.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
``` ```
......