# Este es un ejemplo del calculo de la variable canonica de Fisher
# El codigo asume que se abre R en el directorio en el cual estan grabados los datos
# para evitar tramitar con los "paths"
# Datos y ejemplo del libro "An R and S-plus companion to Multivariate Analysis" por B. Everitt
# Para obtener los datos, bajar el archivo zip disponible en 
# "http://biostatistics.iop.kcl.ac.uk/publications/everitt/"
# 



Tibet<-source("chap7tibetskull.dat")$value; 
print(Tibet);
names(Tibet);
# Leyendo los datos de cada muestra:
attach(Tibet);

X1 <- Tibet[Type==1,-6];
X2 <- Tibet[Type==2,-6];
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
par(mfrow=c(p,p));
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);
		points(Tibet[Type==2,i],Tibet[Type==2,j],xlab="",ylab="",col="blue",pch=16);
		}
	}
}






# Medias y matrix var-cov S_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/coef(disc.an) #(De donde viene esta constante?)

# Standarized canonical variate
std.ys <- (ys-rep(mean(ys,n1+n2)))/sqrt((t(a.vec)%*%S.pool%*%a.vec));

# Standarized canonical variate from R lda ftn.:
std.ys.2 <- (ys.2-rep(mean(ys.2,n1+n2)))/sqrt((t(lmax)%*%S.pool%*%lmax))

# Comparen esto (la varianza de las coordenadas de los puntos sobre a.vec)
sqrt((t(a.vec)%*%S.pool%*%a.vec))

# contra esto
sqrt((t(lmax)%*%S.pool%*%lmax))

sd.avec.por.X <- sqrt((t(a.vec)%*%S.pool%*%a.vec));
sd.avec.por.X


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.varavec <- S.pool%*%a.vec;
var.can1    <- t(a.vec)%*%S.pool%*%a.vec;
can1.covs   <- rbind(cov.varavec,var.can1);
var.cov.mat <- rbind(cbind(S.pool,cov.varavec), 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  

#> var.cov.mat
#             Length    Breadth    Height   Fheight  Fbreadth
#Length    59.013464  9.0080719 17.218889 20.119575 20.110294 -10.9098039
#Breadth    9.008072 48.2605229  1.077222  4.339183 30.046078   0.6196078
#Height    17.218889  1.0772222 36.197778  4.837778  4.108333  -2.7666667
#Fheight   20.119575  4.3391830  4.837778 18.306797 12.985294  -6.6431373
#Fbreadth  20.110294 30.0460784  4.108333 12.985294 43.696078  -7.1470588
#         -10.909804  0.6196078 -2.766667 -6.643137 -7.147059   3.5014415
#> cov.varavec <- S.pool%*%a.vec;
#> var.can1    <- t(a.vec)%*%S.pool%*%a.vec;
#> can1.covs   <- rbind(cov.varavec,var.can1);
#> var.cov.mat <- rbind(cbind(S.pool,cov.varavec), t(can1.covs));
#> 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
#           [,1]       [,2]       [,3]       [,4]       [,5]        [,6]
#[1,]  1.0000000 0.16879527  0.3725535  0.6121207  0.3960241 -0.75895837
#[2,]  0.1687953 1.00000000  0.0257732  0.1459842  0.6542904  0.04766477
#[3,]  0.3725535 0.02577320  1.0000000  0.1879312  0.1033008 -0.24574939
#[4,]  0.6121207 0.14598417  0.1879312  1.0000000  0.4591181 -0.82974301
#[5,]  0.3960241 0.65429041  0.1033008  0.4591181  1.0000000 -0.57780685
#[6,] -0.7589584 0.04766477 -0.2457494 -0.8297430 -0.5778069  1.00000000
























pcomp.trial <- princomp(~Length+Breadth+Height+Fheight+Fbreadth);



#manova.examp <- manova(cbind(Length,Breadth,Height,Fheight,Fbreadth)~Type, data=Tibet);