#MULTIVARIATE DATA: GEOMETRIC INTERPRETATIONS
#READ DATA:
M=read.table("c:/DATA/Multivariate/Isetosastand.txt",header=F)
M

attach(M)

#PLOT OF VECTORS
plot(V1,V2,asp=1)
segments(0,0,V1,V2)

#TWO SPECIFIC VECTORS:
arrows(0,0,V1[6],V2[6],col='red',lwd=3,code=2,length=0.1)
arrows(0,0,V1[20],V2[20],col='blue',lwd=3,code=2,length=0.1)


#CALCULATE PROJECTION:
#HANDY PROJECTION FUNCTION:

#Projects x onto y

proj <- function(x,y){
 res=((t(x)%*%y)/(t(y)%*%y))*y

 return(res)
}

X=t(M)
score=proj(X[,20],X[,6])
score

#PLOT PROJECTION:
plot(V1,V2,asp=1)
arrows(0,0,V1[6],V2[6],col='red',lwd=3,code=2,length=0.1)
arrows(0,0,V1[20],V2[20],col='blue',lwd=3,code=2,length=0.1)
arrows(V1[20],V2[20],score[1],score[2],col='green',code=2,length=0.1)


#PLOT OF POINTS:
plot(V1,V2,asp=1)
points(V1[6],V2[6],col='red',pch=20)
points(V1[20],V2[20],col='blue',pch=20)
arrows(V1[20],V2[20],V1[6],V2[6],col='purple',lwd=3,code=3,length=0.1)

#FUNCTION make.grid() FINDS POINTS ON A GRID:
make.grid <- function(Xsize,Ysize){

ii=0
jj=0
X=numeric(0)
Y=numeric(0)

for(i in 1:Xsize) {
for(j in 1:Ysize) {
ii=ii+1
X[ii]=i
if(j==Ysize) j=0
jj=jj+1
Y[jj] = j
}
}

X
Y

M=cbind(X-mean(X),Y-mean(Y))
return(M)
}


X=make.grid(20,20)

#LINEAR TRANSFORMATION:
M=matrix(c(1.0,0.3,0.3,1),nrow=2,ncol=2,byrow=T)
Z=M%*%t(X)
Z=t(Z)

#PLOT UNTRANSFORMED AND TRANSFORMED POINTS:
plot(Z,type='n',xlab='X',ylab='Y')
points(X,col='blue')
points(Z,col='red',pch=20)

#PLOT TRANSFORMATION VECTORS:
plot(Z,type='n',xlab='X',ylab='Y')
arrows(X[,1],X[,2],Z[,1],Z[,2],col='purple',lwd=1,code=2,length=0.05)


#WARPING SPACE/DISTANCE AS OBSERVED BY A CIRCLE:
#FUNCTION make.circle() FINDS POINTS ON A CIRCLE AROUND ZERO:
make.circle <- function(size,radius){
X=numeric(0)
Y=numeric(0)

for(k in 1:size) {
theta=(1/size)*2*pi*k
X[k]=radius*cos(theta)
Y[k]=radius*sin(theta)
print(paste("theta,X,Y",theta,X,Y))
}
M=cbind(X,Y)
return(M)
}

C=make.circle(200,2)
C

#LINEAR TRANSFORMATION:
E=M%*%t(C)
E=t(E)


#PLOT UNTRANSFORMED AND TRANSFORMED CIRCLE:
plot(E,asp=1,type='n',xlab='X',ylab='Y')
points(0,0,col='black',pch=3)
points(C,col='blue',pch=20)
points(E,col='red',pch=20)



#FUNCTION FOR SQUARE ROOT OF A SYMMETRIC POSITIVE DEFINITE MATRIX:
matrix.square.root <- function(M) {
M.eig <- eigen(M)
M.square.root <- M.eig$vectors %*% diag(sqrt(M.eig$values)) %*% solve(M.eig$vectors)
return(M.square.root)
}


#MAHALANOBIS CORRECTED PLOT:
M=read.table("c:/DATA/Multivariate/Isetosastand.txt",header=F)
attach(M)
S=cov(M)
S

S.sqr=matrix.square.root(S)
S.sqr

M.mah=solve(S.sqr)%*%t(M)
M.mah=t(M.mah)

#PLOT OF POINTS:
plot(M.mah[,1],M.mah[,2],type='n',asp=1,xlab="V1",ylab="V2")
points(V1,V2,col='blue',pch=20)
points(M.mah[,1],M.mah[,2],col='red',pch=20)