#
# Aufgabe 1
#

S11 <- matrix(c(1, 0.49, 0.53, 0.49, 0.51, 
                0.49, 1, 0.57, 0.46, 0.53, 
                0.53, 0.57, 1, 0.48, 0.57, 
                0.49, 0.46, 0.48, 1, 0.57, 
                0.51, 0.53, 0.57, 0.57, 1), 5, 5, byrow=F)


S21 <- matrix(c(0.33, 0.30, 0.31, 0.24, 0.38,  
                0.32, 0.21, 0.23, 0.22, 0.32,  
                0.20, 0.16, 0.14, 0.12, 0.17,  
                0.19, 0.08, 0.07, 0.19, 0.23,  
                0.30, 0.27, 0.24, 0.21, 0.32,  
                0.37, 0.35, 0.37, 0.29, 0.36,  
                0.21, 0.20, 0.18, 0.16, 0.27), 7, 5, byrow=T)
                
S12 <- t(S21)

S22 <- matrix(c(1.00, 0.43, 0.27, 0.24, 0.34, 0.37, 0.40,  
                0.43, 1.00, 0.33, 0.26, 0.54, 0.32, 0.58,  
                0.27, 0.33, 1.00, 0.25, 0.46, 0.29, 0.45,  
                0.24, 0.26, 0.25, 1.00, 0.28, 0.30, 0.27,   
                0.34, 0.54, 0.46, 0.28, 1.00, 0.35, 0.59, 
                0.37, 0.32, 0.29, 0.30, 0.35, 1.00, 0.31,  
                0.40, 0.58, 0.45, 0.27, 0.59, 0.31, 1.00), 7, 7, byrow=T)


ei <- eigen(solve(sqrt.m(S11))%*%S12%*%solve(S22)%*%S21%*%solve(sqrt.m(S11)))

fi <- eigen(solve(sqrt.m(S22))%*%S21%*%solve(S11)%*%S12%*%solve(sqrt.m(S22)))

ei$values                                # canonical correlations

t(ei$vectors[,1])%*%solve(sqrt.m(S11))   # a_i

t(fi$vectors[,1])%*%solve(sqrt.m(S22))   # b_i

sqrt.m<-function(vc)  
{ 
  if(dim(vc)[1] != dim(vc)[2]) stop("Matrix is not square") 
  eg <- eigen(vc) 
  if(min(eg$values) < 0 || max(eg$values)/min(eg$values) > 1000)  
    stop("Variance-covariance matrix is singular") 
  eg$values <- sqrt(eg$values) 
  sqrtvc <- eg$vector %*% diag(eg$values) %*% t(eg$vector) 
  sqrtvc 
}  


#
# Aufgabe 2 (a und b)
#
# siehe auch pp155-157
#

crabs<-read.table("Desktop/crabs.txt",T)

crabsB<-subset(crabs, sp==1)

attach(crabsB)

plot(FL, RW, col=sex)

plot.contours2d(FL[sex==1],RW[sex==1], clevels=0.95, label=F)
plot.contours2d(FL[sex==2],RW[sex==2], clevels=0.95, label=F)

abline(28,-1.26)

plot(FL, RW, col=sex, xlim=c(8,21), ylim=c(6, 19))

rot <- as.matrix(cbind(FL, RW))%*%c(2.93, -3.68)

plot(density(rot[sex==1], bw=2), xlim=c(-15, 10), ylim=c(0,0.15))
lines(density(rot[sex==2], bw=2))

abline(v=-2.6255)

# Aufgabe 3

plot(FL, RW, col=sex, pch=sex, xlim=c(6,23), ylim=c(6, 19))
plot.contours2d(FL[sex==1],RW[sex==1], clevels=0.95, label=F)
plot.contours2d(FL[sex==2],RW[sex==2], clevels=0.95, label=F)

vc<-(var(cbind(FL[sex==1],RW[sex==1]))+var(cbind(FL[sex==2],RW[sex==2])))/2
plot.contours2d(FL[sex==1],RW[sex==1], clevels=0.95, label=F,vc=vc,col="red")
plot.contours2d(FL[sex==2],RW[sex==2], clevels=0.95, label=F,vc=vc,col="red")

plot.contours2d(FL,RW, clevels=0.95, label=F, col="green")

# Aufgabe 4

olives<-read.table("Desktop/olives.txt",T)
attach(olives)

Region.int<-unclass(Region)
plot(stearic,oleic,col=Region.int)
for (i in 1:3)
  plot.contours2d(stearic[Region.int==i],oleic[Region.int==i], clevels=0.95, label=F,col=i)

table(predict(lda(Region~stearic+oleic))$class,Region)