Manipulation matricielle de données

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")

https://freakonometrics.hypotheses.org/files/2016/01/proj1.png

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


Leave a Reply

Your email address will not be published. Required fields are marked *