# PCA and non-normality
pdf("pcanormal.pdf",width=9,height=4.5)
par(mfrow=c(1,2))
par(mar=c(4,4,2,2))
library(mvtnorm)
set.seed(122)
x=scale(rmvnorm(40,c(0,0),matrix(c(1,-.9,-.9,1),ncol=2)))
xtr=scale(cbind(exp(x[,1]),exp(x[,2])))
plot(xtr,xlim=c(-3,3),ylim=c(-3,3),xlab="x1",ylab="x2",cex.lab=1.2)
x.pc=princomp(xtr)
l1=-x.pc$loa[,1]
l2=-x.pc$loa[,2]
c1=3.5;c2=3;arrows(-l1[1]*c1,-l1[2]*c1,l1[1]*c2,l1[2]*c2,length=0.20,angle=20)
#c1=6;c2=6;arrows(-l2[1]*c1,-l2[2]*c1,l2[1]*c2,l2[2]*c2,length=0.20,angle=20)
abline(h=0,lty=3)
abline(v=0,lty=3)
text(2.7,-2.7,expression(u[1]),cex=1.4)
plot(x,xlim=c(-3,3),ylim=c(-3,3),xlab="x1",ylab="log(x2)",cex.lab=1.2)
x.pc=princomp(x)
l1=-x.pc$loa[,1]
l2=-x.pc$loa[,2]
c1=3.5;c2=3;arrows(-l1[1]*c1,-l1[2]*c1,l1[1]*c2,l1[2]*c2,length=0.20,angle=20)
#c1=6;c2=6;arrows(-l2[1]*c1,-l2[2]*c1,l2[1]*c2,l2[2]*c2,length=0.20,angle=20)
abline(h=0,lty=3)
abline(v=0,lty=3)
text(2.7,-2.7,expression(u[1]),cex=1.4)
dev.off()
pdf("pcaout.pdf",width=9,height=4.5)
par(mfrow=c(1,2))
par(mar=c(4,4,2,2))
library(mvtnorm)
set.seed(124)
xin=rmvnorm(25,c(0,0),matrix(c(1,0.9,0.9,1),ncol=2))
xout1=rmvnorm(2,c(-0.5,2),0.2*matrix(c(1,0,0,1),ncol=2))
xout2=rmvnorm(3,c(1,-1.9),0.2*matrix(c(1,0,0,1),ncol=2))
x=rbind(xin,xout1,xout2)
plot(x,xlim=c(-3,3),ylim=c(-3,3),xlab="x1",ylab="x2",cex.lab=1.2)
x.pc=princomp(x,cor=T)
l1=-x.pc$loa[,1]
l2=x.pc$loa[,2]
c1=3;c2=3;arrows(-l1[1]*c1,-l1[2]*c1,l1[1]*c2,l1[2]*c2,length=0.20,angle=20)
#c1=6;c2=6;arrows(-l2[1]*c1,-l2[2]*c1,l2[1]*c2,l2[2]*c2,length=0.20,angle=20)
abline(h=0,lty=3)
abline(v=0,lty=3)
text(2.7,-2.7,expression(u[1]),cex=1.4)
library(ellipse)
lines(ellipse(0.5,scale=c(.12,.22),centre=c(-0.65,1.8)))
lines(ellipse(0.5,scale=c(.12,.22),centre=c(0.9,-2.1)))
library(robustbase)
x.mcd=covMcd(x)
plot(x,xlim=c(-3,3),ylim=c(-3,3),xlab="x1",ylab="x2",cex.lab=1.2)
x.rpc=princomp(x,covmat=x.mcd,cor=T)
l1=-x.rpc$loa[,1]
l2=x.rpc$loa[,2]
c1=3;c2=3;arrows(-l1[1]*c1,-l1[2]*c1,l1[1]*c2,l1[2]*c2,length=0.20,angle=20)
abline(h=0,lty=3)
abline(v=0,lty=3)
text(2.5,2.5,expression(u[1]),cex=1.4)
lines(ellipse(0.5,scale=c(.12,.22),centre=c(-0.65,1.8)))
lines(ellipse(0.5,scale=c(.12,.22),centre=c(0.9,-2.1)))
dev.off()