Pour travailler un peu la manipulation de matrices, dans le cadre du cours d’analyse de données, quelques lignes de code.
La matrice de données est ici
> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",header=TRUE) > M=as.matrix(BD[,2:10]) > colnames(M)=names(BD)[2:10] > rownames(M)=BD$dep > head(M[,1:5]) infra vols eco crim vma D1 44.11 29.51 4.14 2.79 0.08 D2 45.97 26.29 3.97 4.54 0.07 D3 38.83 21.90 4.61 3.33 0.03 D4 49.68 30.11 5.28 4.54 0.07 D5 47.67 27.19 4.94 3.91 0.05 D6 109.21 77.71 6.87 5.25 0.24
La décomposition spectrale de
> round(t(M)%*%M)
infra vols eco crim vma vvi camb
infra 308599 191433 29179 21290 642 7042 30951
vols 191433 119961 17859 12967 404 4427 19335
eco 29179 17859 2946 2042 59 648 2942
crim 21290 12967 2042 1564 42 448 2129
vma 642 404 59 42 2 19 62
vvi 7042 4427 648 448 19 237 665
camb 30951 19335 2942 2129 62 665 3200
roul 55363 34745 5089 3792 113 1217 5612
auto 24841 15747 2214 1676 53 545 2520
(j’ai viré les deux dernières colonnes pour des raisons de présentation) donne ici
> eigen(t(M)%*%M) $values [1] 4.471766e+05 1.239465e+03 3.462936e+02 1.764038e+02 1.303316e+02 6.623592e+01 2.624739e+01 1.761801e+01 1.564292e-01 $vectors
(etc). Pour rappel, les liens entre la diagonalisation et l’analyse spectrale se voit avec l’écriture
> P=eigen(t(M)%*%M)$vectors > P%*%diag(eigen(t(M)%*%M)$values)%*%t(P)
Maintenant, deux remarques. Tout d’abord le lien entre la décomposition spectrale de cette matrice, et la décomposition en valeurs singulières
> sqrt(eigen(t(M)%*%M)$values) [1] 668.7126681 35.2060415 18.6089656 13.2817097 11.4162872 8.1385454 5.1232210 4.1973809 0.3955113 > svd(M) $d [1] 668.7126681 35.2060415 18.6089656 13.2817097 11.4162872 8.1385454 5.1232210 4.1973809 0.3955113 $u [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] -0.08156446 -0.0632514738 -2.596265e-02 0.0960587227 -0.064759395 -0.0430947134 -0.024719699 -0.074347281 0.0550356487 [2,] -0.08087892 0.0541859480 -5.777492e-02 -0.0628249751 -0.099084701 -0.1311683858 0.083484434 0.011872601 -0.0195034904
(etc) et la décomspotion spectrale de l’autre produit matriciel
> sqrt(eigen(M%*%t(M))$values) [1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 1.141629e+01 8.138545e+00 5.123221e+00 4.197381e+00 3.955113e-01 7.546505e-06 6.656212e-06 5.285131e-06 4.270035e-06
(etc).
Maintenant, pour comprendre un peu mieux la recherche de composantes pricinpales, considérons les deux variables
> sM=M[,c(1,3)] > plot(sM)
On a tout intérêt à centrer et réduire les variables (ou changer la métrique, comme vu en cours)
> sMcr=sM > for(j in 1:2) sMcr[,j]=(sMcr[,j]-mean(sMcr[,j]))/sd(sMcr[,j]) > plot(sMcr)
Avant de faire des projections sur un axe, introduisons deux fonctions
> prod_scal=function(x,u){
+ ps=rep(NA,nrow(x))
+ for(i in 1:nrow(x)) ps[i]=sum(x[i,]*u)
+ return(ps)
+ }
> proj=function(x,u){
+ px=x
+ for(j in 1:length(u)){
+ px[,j]=prod_scal(x,u)/sqrt(sum(u^2))*u[j]
+ }
+ return(px)
+ }
Par exemple, si on projette sur l’axe des abscisses,
> points(proj(sMcr,c(1,0)),col="blue")
On peut ensuite chercher la direction de l’axe qui nous donnera le nuage de points qui aura la plus grande inertie
> inertie=function(x) sum(x^2) > Theta=seq(0,3.141592657,length=101) > V=unlist(lapply(Theta,function(theta) inertie(proj(sMcr,c(cos(theta),sin(theta)))))) > plot(Theta,V,type='l')
> (angle=optim(0,function(theta) -inertie(proj(sMcr,c(cos(theta), sin(theta)))))$par) [1] 0.7853516
Visuellement, on obtient
> plot(sMcr) > psMcr=proj(sMcr,c(cos(angle),sin(angle))) > points(psMcr,col="blue")
On notera que l’axe qui donne l’inertie maximale est lié aux vecteurs propres de la décomposition spectrale (celui associé à la plus grande valeur propre)
> c(cos(angle),sin(angle))
[1] 0.7071397 0.7070738
> eigen(t(sMcr)%*%sMcr)
$values
[1] 159.64663 26.35337
$vectors
[,1] [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068 0.7071068
On continuera demain, à manipuler des matrices de données, de faire des projections, avant de commencer l’analyse en composantes principales, à proprement parler, la semaine prochaine
