# Ilustracja skalowania wielowymiarowego (MDS) metoda najmniejszych kwadratow
# i porownanie jego z PCA. Samce oznaczone sa trojkatami. Rezultatem programu 
# jest rysunek w .eps. Funkcja system.time podaje czas obliczen: na procesorze 
# PIII 1GHz trwaly one ponad 20 min.

s=function(Z){ Z=matrix(Z,ncol=2); sum( (DX-dist(Z))^2 ) }

print(system.time({

  library(MASS); X=as.matrix(crabs[,4:8])
  DX<<-dist(X)
  HX=t(t(X)-apply(X,2,mean))
  Z1=HX%*%svd(HX)$v[,1:2]
  Z2=matrix( nlm(s,as.vector(Z1))$estimate, ncol=2)
  Z3=matrix( nlm(s,as.vector(HX[,1:2]))$estimate, ncol=2)

}))

cl=c(rep(1,50),rep(2,50),rep(3,50),rep(4,50))
sex= c(rep(2,50),rep(1,50),rep(2,50),rep(1,50))
postscript("mds.eps",pointsize=10)
par(mfrow=c(1,3),pty="s",cex=1.2,las=1)
plot(Z1,col=cl,xlab="pc1",ylab="pc2",pch=sex,cex.main=.8,main="PCA")
plot(Z2,col=cl,xlab="z1",ylab="z2",pch=sex,cex.main=.8,
        main="Least Squares MDS started from PCA")
plot(Z3,col=cl,xlab="z1",ylab="z2",pch=sex,cex.main=.8,
        main="Least Squares MDS started from HX[,1:2]")
graphics.off()