# Objetivo: Calcular la variable canonica de Fisher para los datos de medidas de craneos humanos de 2 poblaciones del Tibet
# Datos y ejemplo del libro "An R and S-plus companion to Multivariate Analysis" por B. Everitt

# Lectura de datos:
datos <- read.table("http://www.cimat.mx/~ponciano/STAT90MOES02/TibetMatrix.txt", header=T);
Tibet <- data.frame(datos);
print(Tibet); # imprimir datos
names(Tibet); # nombres de las variables en Tibet

# cargar los datos como un data.frame:
attach(Tibet);
X1 <- Tibet[Type==1,-6]; # los datos del 'tipo 1', exceptuando la sexta columna (es la variable 'Type')
X2 <- Tibet[Type==2,-6]; # los datos del 'tipo 2', exceptuando la sexta columna (es la variable 'Type')
n1 <- dim(X1)[1];
n2 <- dim(X2)[1];
p  <- dim(X1)[2];


# Graficas exploratorias:
pairs(Tibet[,-6],col="red");
# Ahora lo mismo pero coloreando rojo los datos de tipo 1 y azul los datos de tipo 2:
# Preguntarme o averiguar usando 'help(par)' que es cada argumento dentro de estas dos bucles 'for'.
# La funcion 'par' determina los parametros para graficar lo que uno quiere. 'mfrow' significa 'multi-figure, by row'
par(mfrow=c(p,p),mai=c(0.25,0.25,0.25,0.25),oma=c(1.45,1.5,0.5,0.55));
for(i in 1:p){
	for(j in 1:p){
		if(j==i){plot(0,0,pch=".",main="",xlab="",ylab="",xaxt="n",yaxt="n");
			     text(0,0,names(Tibet)[i], cex=2);
			     }else{ 
		plot(Tibet[Type==1,i],Tibet[Type==1,j],xlab="",ylab="",col="red",pch=16, cex.axis=1.5);
		points(Tibet[Type==2,i],Tibet[Type==2,j],xlab="",ylab="",col="blue",pch=16, cex.axis=1.5);
		}
	}
}

# Medias y matriz de var-cov conjunta:

X1.means <- apply(X1,2,mean);
X2.means <- apply(X2,2,mean);
S.pool  <- (1/(n1+n2-2))*((n1-1)*var(X1) + (n2-1)*var(X2));


# Calculando los coeficientes de la var. canon. de Fisher:
d.vec <- X1.means -X2.means;
Sm1   <- solve(S.pool);
a.vec <- Sm1%*%d.vec; # Este es el vector de coefs.

# Coord. de los puntos en la var. canon. :

y1s <- t(a.vec)%*%t(as.matrix(X1));
y2s <- t(a.vec)%*%t(as.matrix(X2));
ys  <- c(y1s,y2s);

min.ys <- min(ys); max.ys <- max(ys);
hist(y1s,xlim=c(min.ys,max.ys));X11();hist(y2s,xlim=c(min.ys,max.ys));


#----------- OK: comparando con el output de la funcion 'lda' de la libreria MASS----#

library(MASS)

disc.an <- lda(Type~Length+Breadth+Height+Fheight+Fbreadth, data=Tibet, prior=c(0.5,0.5));

lmax <- matrix(coef(disc.an),nrow=p,ncol=1); # Este es el vector 'a' calculado por la funcion 'lda'
y1s.2 <- t(lmax)%*%t(as.matrix(X1));
y2s.2 <- t(lmax)%*%t(as.matrix(X2));
ys.2  <- c(y1s.2,y2s.2);


# Lmax difiere de a.vec por una constante:
a.vec/lmax #(De donde viene esta constante?)

# es nada menos que la varianza estimada de a'X: 
sqrt((t(a.vec)%*%S.pool%*%a.vec))

# Entonces si estandarizamos las ys:
std.ys <- (ys-rep(mean(ys,n1+n2)))/sqrt((t(a.vec)%*%S.pool%*%a.vec));
std.ys

# Encontramos lo mismo que si estandarizamos las ys.2 que entrega "lda":
std.ys.2 <- (ys.2-rep(mean(ys.2,n1+n2)))/sqrt((t(lmax)%*%S.pool%*%lmax))


# Noten que la varianza de las y's entregadas por "lda" es 1: "lda" divide las ys por la varianza
# pero no las centra!!!!
sqrt((t(lmax)%*%S.pool%*%lmax))


# Para verificar, lmax multiplicado por la varianza de a'X es igual a a.vec
sd.avec.por.X <- sqrt((t(a.vec)%*%S.pool%*%a.vec));
lmax*as.numeric(sd.avec.por.X)
a.vec
# y entonces puedo encontrar las ys originales asi:
t(lmax*as.numeric(sd.avec.por.X))%*%t(as.matrix(rbind(X1,X2)))
ys

# Computing the correlation matrix between variables and a.vec:
cov.yapX    <- S.pool%*%a.vec;
var.can1    <- t(a.vec)%*%S.pool%*%a.vec;
can1.covs   <- rbind(cov.yapX,var.can1);
var.cov.mat <- rbind(cbind(S.pool,cov.yapX), t(can1.covs));
var.cov.mat
SDm1        <- matrix(rep(0,(p+1)*(p+1)),nrow=p+1,ncol=p+1);
diag(SDm1)  <- 1/sqrt(diag(var.cov.mat));
Rmat        <- SDm1%*%var.cov.mat%*%SDm1;
Rmat