```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(FactoMineR)
data(iris)
X <- as.matrix(iris[1:10,1:4])
```
 Centrage de la matrice :
```{r}
X <- sweep(X,2,apply(X,2,mean),FUN="-")  # ou X <- scale(X,scale=FALSE)
n <- nrow(X)
summary(X)
```

## SVD et diagonalisation

Corrélation et calcul matriciel (attention, X doit bien être centrée)
```{r}
COV <- cov(X) * (n-1) / n
round(((t(X)%*%X) / n) - COV , 8)
```

```{r}
diagCOV <- eigen(t(X)%*%X)
diagPS <- eigen(X%*%t(X))
SVD <- svd(X)
```

Les valeurs prores sont identiques et égales au carré des valeurs singulières
```{r}
L1 <- diagCOV$values
L2 <- diagPS$values
L1
L2
SVD$d^2
```

```{r}
U <- diagPS$vectors
V <- diagCOV$vectors
round(t(U)%*%U,8)
round(U%*%t(U),8)
round(t(V)%*%V,8)
round(V%*%t(V),8)
```

Lien entre les vecteurs propres et les résultats de la svd
```{r}
SVD$u[,1:4]/U[,1:4]
SVD$v[,1:4]/V[,1:4]
```

```{r}
round(t(SVD$u)%*%SVD$u,8)
round(SVD$u%*%t(SVD$u),8)  ## ici ce n'est plus diagonale car on restreint au vp non nulles
round(t(SVD$v)%*%SVD$v,8)
round(SVD$v%*%t(SVD$v),8)
```

## Reconstitution des données
```{r}
print(sum((X - SVD$u[,1,drop=FALSE]%*%SVD$d[1]%*%t(SVD$v[,1,drop=FALSE]))^2))  # attention, pour diag Lambda quand 1 seule valeur
for (j in 2:ncol(X)) print(sum((X - SVD$u[,1:j,drop=FALSE]%*%diag(SVD$d[1:j])%*%t(SVD$v[,1:j,drop=FALSE]))^2))
```

## L'ACP normée
```{r}
data(iris)
X <- as.matrix(iris[1:10,1:4])
X <- scale(X,scale=FALSE)
n <- nrow(X)
M <- diag(1/(apply(X,2,var)*(n-1)/n)) ## métrique pour normer
N <- diag(rep(1/n,n))

diagPS <- eigen(X%*%M%*%t(X)%*%N)
round(diagPS$vectors%*%t(diagPS$vectors), 8) ## U'U=Id
round(t(diagPS$vectors)%*%diagPS$vectors, 8) ## UU'=Id

SVD <- svd(N^0.5%*%X %*% M^(0.5))
SVD$u/diagPS$vectors[,1:4]

pca <- PCA(X,gr=F)

pca$eig[,1] ; diagPS$values[1:4]; SVD$d^2

coordInd <- SVD$u[,1:4] %*% diag(SVD$d)
coordInd*sqrt(n)/ pca$ind$coord
coordInd2 <- N^0.5%*%X%*%M^0.5%*%SVD$v
coordInd2 / pca$ind$coord

coordVar <- SVD$v %*% diag(SVD$d) 
coordVar[,1:4]/ pca$var$coord

coordVar2 <- (M^0.5%*%t(X)%*%N^0.5%*%SVD$u)
coordVar2[,1:4]/pca$var$coord
```

## AFC comme une svd particulière
```{r}
data(children)
X <- as.matrix(children[1:8,1:5])
ca <- CA(X, graph=FALSE)

P = X/sum(X)
r = apply(P,1,sum)
c = apply(P,2,sum)
invDr=diag(1/r)
invDc=diag(1/c)

S = invDr^0.5%*%(P-r%*%t(c))%*%invDc^0.5
res <- svd(S)
F <- invDr^0.5%*%res$u[,1:4]%*%diag(res$d[1:4])
F/ca$row$coord

G <- invDc^0.5%*%res$v[,1:4]%*%diag(res$d[1:4])
G/ca$col$coord
```

## ACM comme une SVD pondérée

```{r}
library(FactoMineR)
data(tea)
X <- tea[,10:18]
n <- nrow(X)
Z <- tab.disjonctif(X)
M <- diag(n/apply(Z,2,sum))
Z <- scale(Z, scale=FALSE)
N <- diag(rep(1/n,n))
U <- eigen(Z%*%M%*%t(Z)%*%N)
V <- eigen(t(Z)%*%N%*%Z%*%M)
invDsigma <- solve(M)

S <- svd(N^0.5 %*% Z%*%invDsigma^0.5 / ncol(X))

mca <- MCA(X,graph=FALSE)

mca$eig[1:6,1] * ncol(X) / U$val[1:6]
U$vect[,1:4]%*%diag(sqrt(U$val[1:4]))/mca$ind$coord[,1:4]*sqrt(n/ncol(X)) # F = ZMV = U Lambda
M%*%V$vect[,1:4]%*%diag(sqrt(U$val[1:4]))/mca$var$coord[,1:4] #G = MV Lambda

```

## AFC et ACM
```{r}
data(tea)
don <- tea[,c(14,18)]
mca <- MCA(don)
ca <- CA(table(don))
I <- nrow(table(don))
J <- ncol(table(don))
ncp <- min(5,I-1,J-1)

sweep(mca$var$coord[1:I,1:ncp],2,sqrt(mca$eig[1:ncp,1]),FUN="/")/sweep(ca$row$coord[,1:ncp],2,sqrt(ca$eig[1:ncp,1]),FUN="/")
sweep(mca$var$coord[(I+1):(I+J),1:ncp],2,sqrt(mca$eig[1:ncp,1]),FUN="/")/sweep(ca$col$coord[,1:ncp],2,sqrt(ca$eig[1:ncp,1]),FUN="/")
```

