from numpy import *
# This line configures matplotlib to show figures embedded in the notebook,
%matplotlib inline
import matplotlib.pyplot as plt
import pylab as pl # import MATLAB-like API 

N = 400
X = random.randn(N,2) # distribution normale (0,1) avec tirages indépendants

T = array([[-1,1],[4,2]]) # matrice quelconque 2x2 
X = dot(X,T) # on transforme X linéairement X*T

Xc = zeros(X.shape)
for k in range(0,2):
    Xc[:,k] = X[:,k] - mean(X[:,k])

pl.plot([-30, 30],[0, 0],'k')
pl.plot([0, 0],[-20, 20],'k')
pl.plot(Xc[:,0],Xc[:,1],'b.')
pl.plot(Xc[1,0],Xc[1,1],'ro') # X[1,:] tagué rouge 
pl.plot(Xc[8,0],Xc[8,1],'co') # X[8,:] tagué cyan 
pl.axis('equal')
pl.axis([-10,10,-10,10])

U, s, Vs = linalg.svd(Xc, full_matrices=True)

s, Vs

P = dot(Xc,Vs.T)

pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge 
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan 
pl.axis('equal')
pl.axis([-10,10,-10,10])

len1 = s[0]*2/sqrt(N-1) # 2 écarts types le long du semi-axe u1
len2 = s[1]*2/sqrt(N-1) # 2 écarts types le long du semi-axe u2

pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot([0,len1],[0,0],'g',linewidth=4)
pl.plot([0,0],[0,len2],'g',linewidth=4)
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge 
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan 
pl.axis('equal')
pl.axis([-10,10,-10,10])
pl.xlabel('premiere composante')
pl.ylabel('deuxieme composante')

theta = linspace(0,2*pi,101)
circle = zeros((101,2))
circle[:,0] = cos(theta)
circle[:,1] = sin(theta)

ellipse = dot(circle*s*2/sqrt(N-1),Vs)

pl.plot([-30, 30],[0, 0],'k') # axe x
pl.plot([0, 0],[-20, 20],'k') # axe y
pl.plot(Xc[:,0],Xc[:,1],'b.') # nuage de points
pl.plot(ellipse[:,0],ellipse[:,1],'g',linewidth=4) # ellipse avec semi-axes orientés par V et de longueurs 2*sigma/sqrt(N-1)
pl.plot([0,Vs[0,0]*len1],[0,Vs[0,1]*len1],'g',linewidth=4) # semi-axes orientés par V1 et de longueur 2*sigma1/sqrt(N-1)
pl.plot([0,Vs[1,0]*len2],[0,Vs[1,1]*len2],'g',linewidth=4) # semi-axes orientés par V2 et de longueur 2*sigma2/sqrt(N-1)
pl.plot(Xc[1,0],Xc[1,1],'ro') # X[1,:] tagué rouge 
pl.plot(Xc[8,0],Xc[8,1],'co') # X[8,:] tagué cyan 
pl.axis('equal') # ratio d'aspect des axes 1:1 important pour l'interprétation géométrique 
pl.axis([-10,10,-10,10])
pl.xlabel('X1')
pl.ylabel('X2')

pl.plot(P[:,0],P[:,1],'b.')
pl.plot(array([-30,30]),array([0, 0]),'k')
pl.plot(array([0, 0]),array([-30,30]),'k')
pl.plot(circle[:,0]*len1,circle[:,1]*len2,'g',linewidth=4) # ellipse avec semi-axes le long des composantes principales
pl.plot(P[1,0],P[1,1],'ro') # X[1,:] tagué rouge 
pl.plot(P[8,0],P[8,1],'co') # X[8,:] tagué cyan 
pl.axis('equal')
pl.axis([-10,10,-10,10])
pl.xlabel('premiere composante')
pl.ylabel('deuxieme composante')