Commit 74341c64 by Eric Coissac

Ecriture du modèle

parent f96cd20c
---
title: "Untitled"
author: "Aurlie & Eric"
author: "Christelle & Eric"
date: "27/03/2018"
output: html_document
output:
html_document: default
---
```{r setup, include=FALSE}
......@@ -33,21 +34,21 @@ env = sweep(env,MARGIN = 2,colMeans(env),'-')
env = sweep(env,MARGIN = 2,apply(env,2,sd),'/')
```
```{r eval=FALSE, include=FALSE}
```{r}
euk=sqrt(euk)
euk = sweep(euk,MARGIN = 2,colMeans(euk),'-')
euk = sweep(euk,MARGIN = 2,apply(euk,2,sd),'/')
#euk = sweep(euk,MARGIN = 2,colMeans(euk),'-')
#euk = sweep(euk,MARGIN = 2,apply(euk,2,sd),'/')
```
```{r eval=FALSE, include=FALSE}
```{r}
bac=sqrt(bac)
bac = sweep(bac,MARGIN = 2,colMeans(bac),'-')
bac = sweep(bac,MARGIN = 2,apply(bac,2,sd),'/')
#bac = sweep(bac,MARGIN = 2,colMeans(bac),'-')
#bac = sweep(bac,MARGIN = 2,apply(bac,2,sd),'/')
```
```{r}
euk.dist=quasieuclid(vegdist(euk,method = "jaccard"))
euk.dist=vegdist(euk,method = "jaccard")
euk.pco =dudi.pco(euk.dist,full = TRUE)
euk.pco.li = euk.pco$li
......@@ -57,7 +58,7 @@ text(euk.pco.li[,1:2],labels = rownames(euk.pco.li),cex=0.6)
```
```{r}
bac.dist=quasieuclid(vegdist(bac,method = "jaccard"))
bac.dist=vegdist(bac,method = "jaccard")
bac.pco =dudi.pco(bac.dist,full = TRUE)
bac.pco.li = bac.pco$li
......@@ -83,21 +84,233 @@ y = procrustes(bac.pco.li,euk.pco.li,scale = TRUE)
```{r}
SCR.E.B = sum((euk.pco.li-x$Yrot)^2)
SCT.E = sum((euk.pco.li- colMeans(euk.pco.li))^2)
SCR.E.B/SCT.E
1-SCR.E.B/SCT.E
```
```{r}
SCR.B.E = sum((bac.pco.li-y$Yrot)^2)
SCT.B = sum((bac.pco.li- colMeans(bac.pco.li))^2)
SCR.B.E/SCT.B
1-SCR.B.E/SCT.B
```
```{r}
x = procrustes(sweep(euk.pco.li,2,colMeans(euk.pco.li),"-"),
sweep(bac.pco.li,2,colMeans(bac.pco.li),"-"),scale = TRUE)
y = procrustes(sweep(bac.pco.li,2,colMeans(bac.pco.li),"-"),
sweep(euk.pco.li,2,colMeans(euk.pco.li),"-"),scale = TRUE)
myproc = function (X, Y, scale = TRUE, symmetric = FALSE, scores = "sites",
...)
{
X <- scores(X, display = scores, ...)
Y <- scores(Y, display = scores, ...)
if (nrow(X) != nrow(Y))
stop("Matrices have different number of rows: ", nrow(X),
" and ", nrow(Y))
if (ncol(X) < ncol(Y)) {
warning("X has fewer axes than Y: X adjusted to comform Y\n")
addcols <- ncol(Y) - ncol(X)
for (i in 1:addcols) X <- cbind(X, 0)
}
ctrace <- function(MAT) sum(MAT^2)
c <- 1
if (symmetric) {
X <- scale(X, scale = FALSE)
Y <- scale(Y, scale = FALSE)
X <- X/sqrt(ctrace(X))
Y <- Y/sqrt(ctrace(Y))
}
xmean <- apply(X, 2, mean)
ymean <- apply(Y, 2, mean)
if (!symmetric) {
X <- scale(X, scale = FALSE)
Y <- scale(Y, scale = FALSE)
}
XY <- crossprod(X, Y)
sol <- svd(XY)
A <- sol$v %*% t(sol$u)
print(A)
if (scale) {
c <- sum(sol$d)/ctrace(Y)
}
Yrot <- c * Y %*% A
b <- xmean - c * ymean %*% A
R2 <- ctrace(X) + c * c * ctrace(Y) - 2 * c * sum(sol$d)
reslt <- list(Yrot = Yrot, X = X, ss = R2, rotation = A,
translation = b, scale = c, xmean = xmean, symmetric = symmetric,
call = match.call())
reslt$svd <- sol
class(reslt) <- "procrustes"
reslt
}
```
```{r}
myproc3 = function (X, Y, Z)
{
if (nrow(X) != nrow(Y))
stop("Matrices have different number of rows: ", nrow(X),
" and ", nrow(Y))
if (nrow(X) != nrow(Z))
stop("Matrices have different number of rows: ", nrow(X),
" and ", nrow(Y))
ctrace <- function(MAT) sum(MAT^2)
c <- c(1,1)
X <- scale(X, scale = FALSE)
Y <- scale(Y, scale = FALSE)
Z <- scale(Z, scale = FALSE)
xmean <- apply(X, 2, mean)
ymean <- apply(Y, 2, mean)
x$scale
1/y$scale
XY <- crossprod(X, Y)
sol1 <- svd(XY)
A1 <- sol1$v %*% t(sol1$u)
XZ <- crossprod(X, Z)
sol2 <- svd(XZ)
A2 <- sol2$v %*% t(sol2$u)
Yrot <- Y %*% A1
Zrot <- Z %*% A2
CovXY = crossprod(X, Yrot)
CovXY = crossprod(X, Zrot)
CovYZ = crossprod(Yrot, Zrot)
b <- xmean - c * ymean %*% A
R2 <- ctrace(X) + c * c * ctrace(Y) - 2 * c * sum(sol$d)
reslt <- list(Yrot = Yrot, X = X, ss = R2, rotation = A,
translation = b, scale = c, xmean = xmean, symmetric = symmetric,
call = match.call())
reslt$svd <- sol
class(reslt) <- "procrustes"
reslt
}
```
## Estimation des coefficients lors d'une rgression linaire multiple (2 variables)
$$ \hat{y}=a_1 x_1 + a_2 x_2 + b$$
On cherche minimiser l'erreur commise entre la prdiction et la vraie valeur (mthode des moindres carrs)
$$
\left\{
\begin{array}{rcr}
\sum_i{(y_i- \hat{y_i})^2}=0 \\
\sum_i{(y_i- (a_1 x_{1,i} + a_2 x_{2,i} + b))^2}=0
\end{array}
\right.
$$
Pour trouver les paramtres qui minimise cette fonction, on calcule les trois drives partielles et on les annule :
$$
\left\{
\begin{array}{rcr}
-2 \sum_i{x_{1,i}(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0 \\
-2 \sum_i{x_{2,i}(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0 \\
-2 \sum_i{(y_i- a_1 x_{1,i} + a_2 x_{2,i} + b)}=0
\end{array}
\right.
$$
A partir de la dernire quation, on obtient :
$$
\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}-b=0
\\
\Leftrightarrow \hat{b}=\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}
$$
En remplaant dans la premire quation on obtient :
$$
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-b\frac{\sum_i{x_{1,i}}}{n}=0
\\
\Leftrightarrow
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-(\overline{y}-a_1\overline{x_1}-a_2\overline{x_2}) \overline{x_1}=0
\\
\Leftrightarrow
\frac{\sum_i{x_{1,i}y_i}}{n}-a_1\frac{\sum_i{x_{1,i}^2}}{n}-a_2\frac{\sum_i{x_{1,i}x_{2,i}}}{n}-(\overline{y}-a_1\overline{x_1}-a_2\overline{x_2} ) \overline{x_1}=0
\\
\Leftrightarrow
S_{xy}-a_1S^{2}_{x_{1}}-a_2S_{x_1x_2}=0
\\
\Leftrightarrow
\hat{a}_1=\frac{S_{x_1y}-a_2S_{x_1x_2}}{S^{2}_{x_{1}}}
$$
De la mme manire, on obtient :
$$
\hat{a}_2=\frac{S_{x_2y}-a_1S_{x_1x_2}}{S^{2}_{x_{2}}}
$$
## Calcul matriciel :
Soit X la matrice expliquer et Y la matrice explicative.
Par procuste on a trois oprations :
* translation (revient aligner les barycentres : centrer )
* chelle : homothtie
* rotation trouver l'angle qui minimise : mthode des moindres carrs
Seule contrainte de la mthode, la rotation doit tre effectue en dernier.
On cherche : $$ Min|X-(1c'+\rho Y A)| $$
On se retrouve minimiser comme en relation linaire l'cart entre X notre tableau de donnes explicative et la transformation de la matrice Y (qui reprsente l'approximation en rgression linaire) :
$$1c'+\rho Y A $$
avec $c'$ translation (centrage), $\rho$ chelle (normalisation), $A$ rotation
D'un point de vue matriciel : il faut que A'A=I orthogonal, diag(A'A)=1
Les matrices centres de X et Y sont notes :
$$\tilde{Y}, \tilde{X}$$
Dcomposition en valeur singulire : avec U'U=V'V=I
$$\tilde{X}' \tilde{Y}=U \Lambda V'$$
La matrice de rotation est donc estime par : $$\hat{A}=VU'$$
Paramtre d'chelle qui correspond c dans procruste :
$$\hat{\rho}=\frac{trace(\hat{A}\tilde{X}'\tilde{Y})}{trace(\tilde{Y}'\tilde{Y})}$$
Si on fait l'analogie avec la fonction programme dans procrustre (quivalence : translation <=> b, chelle <=> c et rotation A) :
$$Yrot=\hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
$$
? j'aurai plus tt crit, la translation a dj du tre faite dans procruste avant cette tape??
$$Yrot=1c'+ \hat{\rho} Y A
\\
b=\overline{X}-\hat{\rho}\overline{Y} A
$$
? Yrot est la nouvelle estimation de Y qui a t transform pour "tre reprsent dans la mme base" que X. ???
Et du coup la diffrence entre X et Yrot reprsente la part non explique par Y du tableau de donnes X
$$
\sum(X-Y_{rot})^2=SRC_{residuel}
$$
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment