#install.packages("mvtnorm")
#install.packages("scatterplot3d")
library(mvtnorm)
set.seed(100)
n <- 200
mu = c(0,0)
rho = 0.5
Sig = matrix(c(1,rho,rho,1),nrow=2)
x <- rmvnorm(n,mu,Sig)

# plot simulated ata 
dev.off()
par(mfrow=c(1,1)) 
plot(x[,1],x[,2],,xlab="v_1",ylab="v_2",
     main="n (= 200) data points in 2D",asp = 1,
     xlim = c(-3,3), ylim=c(-3,3))
dev.copy(pdf,'PCA_ex1.pdf')
dev.off()

# add directions and see their projections

## directions of interest
d1 = c(cos(0),sin(0))
d2 = c(cos(pi/4),sin(pi/4))
d3 = c(cos(pi/2),sin(pi/2))
d11 <- cbind(d1,-d1)
d22 <- cbind(d2,-d2)
d33 <- cbind(d3,-d3)

par(mfrow=c(1,1)) 
plot(x[,1],x[,2],,xlab="v_1",ylab="v_2",
     main="n (= 200) data points in 2D",asp = 1,
     xlim = c(-3,3), ylim=c(-3,3))
lines(3*d11[1,],3*d11[2,],col="red",lwd=3)
dev.copy(pdf,'PCA_ex1_dir1.pdf')
proj_scores = (x%*%d1)
proj = proj_scores %*% d1 
var_proj = round(var(proj_scores),digits = 2)
points(proj[,1],proj[,2],col="red")
dev.copy(pdf,'PCA_ex1_dir2.pdf')
dev.off()

par(mfrow=c(1,2)) 
plot(x[,1],x[,2],,xlab="v_1",ylab="v_2",
     main="n (= 200) data points in 2D",asp = 1,
     xlim = c(-3,3), ylim=c(-3,3),cex.main=0.8)
lines(3*d11[1,],3*d11[2,],col="red",lwd=3)  
points(proj[,1],proj[,2],col="red") 
plot.density(density(proj_scores), main = 'Projection Scores',cex.main=0.8) 
points(proj_scores,jitter(array(0.2,c(n,1)), factor = 5, amount = NULL),col="red")
mtext("variance",side=1,line=2,at=-2,font=2)
mtext(var_proj,side=1,line=2,at=1,font=2)
dev.copy(pdf,'PCA_ex1_dir3.pdf')
dev.off()

## now plot all 
par(mfrow=c(2,2)) 
plot(x[,1],x[,2],,xlab="v_1",ylab="v_2",
     main="n (= 200) data points in 2D",asp = 1,
     xlim = c(-3,3), ylim=c(-3,3),cex.main=0.8)
lines(3*d11[1,],3*d11[2,],col="red",lwd=3)   
lines(3*d22[1,],3*d22[2,],col="blue",lwd=3)   
lines(3*d33[1,],3*d33[2,],col="green",lwd=3)   

plot.density(density(proj_scores), main = 'Projection Scores',cex.main=0.8) 
points(proj_scores,jitter(array(0.2,c(n,1)), factor = 5, amount = NULL),col="red")
mtext("variance",side=1,line=2,at=-2,font=2)
mtext(var_proj,side=1,line=2,at=1,font=2)

proj_scores = (x%*%d2) 
var_proj = round(var(proj_scores),digits = 2)

plot.density(density(proj_scores), main = 'Projection Scores',cex.main=0.8) 
points(proj_scores,jitter(array(0.2,c(n,1)), factor = 5, amount = NULL),col="blue")
mtext("variance",side=1,line=2,at=-2,font=2)
mtext(var_proj,side=1,line=2,at=1,font=2)


proj_scores = (x%*%d3) 
var_proj = round(var(proj_scores),digits = 2)
plot.density(density(proj_scores), main = 'Projection Scores',cex.main=0.8) 
points(proj_scores,jitter(array(0.2,c(n,1)), factor = 5, amount = NULL),col="green")
mtext("variance",side=1,line=2,at=-2,font=2)
mtext(var_proj,side=1,line=2,at=1,font=2)
dev.copy(pdf,'PCA_ex1_dir4.pdf')
dev.off()

### 3D examples

n <- 200
mu = c(0,0,0)
rho = 0.5
Sig = matrix(c(1.2,rho,0,rho,1,0,0,0,0.4),nrow=3)
xx <- rmvnorm(n,mu,Sig)

par(mfrow=c(1,1)) 
library(scatterplot3d) 
scatterplot3d(xx[,1],xx[,2],xx[,3], main="3D Scatterplot",xlim=c(-3,3),ylim=c(-3,3),zlim=c(-3,3))