# The MIT License
# Copyright (c) 2007 The GGobi Foundation
# http://www.ggobi.org/book/code-license.txt


# Read data, and subset
d.crabs<-read.csv("australian-crabs.csv")
d.blue.crabs <- subset(d.crabs,species=="Blue", select=c(sex,FL:BD))

# Load libraries
library(mclust)
library(rggobi)

# Utility functions
f.vc.ellipse <- function(vc, xm, n=500) {
  p<-ncol(vc)
  x<-f.gen.sphere(n,p)

  evc<-eigen(vc)
  vc2<-(evc$vectors)%*%diag(sqrt(evc$values))%*%t(evc$vectors)
  x<-x%*%vc2

  x + matrix(rep(xm, each=n),ncol=p)
}

f.gen.sphere<-function(n=100,p=5) {
  x<-matrix(rnorm(n*p),ncol=p)
  xnew<-t(apply(x,1,f.norm.vec))
  xnew
}

f.norm.vec<-function(x) {
  x<-x/f.norm(x)
  x
}

f.norm<-function(x) { sqrt(sum(x^2)) }

# Model-based clustering
blue.crabBIC <- mclustBIC(subset(d.blue.crabs,select=c(FL,RW)),
    modelNames=c("EEE","EEV","VVV"))
blue.crabBIC
plot(blue.crabBIC)

 # Examine 2 cluster, equal volume, equal shape, variable orientation model
mclst1 <- mclustBIC(subset(d.blue.crabs,select=c(FL,RW)), G=2,
    modelNames="EEV")
smry1 <- mclustModel(subset(d.blue.crabs,select=c(FL,RW)), mclst1,
    G=2, modelNames="EEV")
cl1<-summary(smry1)$classification
table(d.blue.crabs[,1],cl1)
 # Compute the ellipses and plot
vc<-smry1$parameters$variance$sigma[,,1]
xm<-smry1$parameters$mean[,1]
y1<-f.vc.ellipse(vc,xm)
vc<-smry1$parameters$variance$sigma[,,2]
xm<-smry1$parameters$mean[,2]
y2<-f.vc.ellipse(vc,xm)

 # Examine 3 cluster, equal volume, equal shape, variable orientation model
mclst2<-mclustBIC(subset(d.blue.crabs,select=c(FL,RW)),G=3,modelNames="EEV")
smry2<-mclustModel(subset(d.blue.crabs,select=c(FL,RW)),mclst2,G=3,modelNames="EEV")
cl2<-summary(smry2)$classification
table(d.blue.crabs[,1],cl2)
 # Compute the ellipses and plot
vc<-smry2$parameters$variance$sigma[,,1]
xm<-smry2$parameters$mean[,1]
y1<-f.vc.ellipse(vc,xm)
vc<-smry2$parameters$variance$sigma[,,2]
xm<-smry2$parameters$mean[,2]
y2<-f.vc.ellipse(vc,xm)
vc<-smry2$parameters$variance$sigma[,,3]
xm<-smry2$parameters$mean[,3]
y3<-f.vc.ellipse(vc,xm)

 # Examine 2 cluster, variable volume, variable shape, 
 # variable orientation model
mclst3<-mclustBIC(subset(d.blue.crabs,select=c(FL,RW)),G=2,modelNames="VVV")
smry3<-mclustModel(subset(d.blue.crabs,select=c(FL,RW)),mclst3,G=2,modelNames="VVV")
cl3<-summary(smry3)$classification
table(d.blue.crabs[,1],cl3)
 # Compute the ellipses and plot
vc<-smry3$parameters$variance$sigma[,,1]
xm<-smry3$parameters$mean[,1]
y1<-f.vc.ellipse(vc,xm)
vc<-smry3$parameters$variance$sigma[,,2]
xm<-smry3$parameters$mean[,2]
y2<-f.vc.ellipse(vc,xm)

 # All 5 variables
mclst4 <- mclustBIC(subset(d.crabs,select=c(FL:BD)),
    G=1:8, modelNames=c("EEE","EEV","VVV"))
plot(mclst4)

mclst5<-mclustBIC(subset(d.crabs,select=c(FL:BD)),G=4,modelNames="EEV")
smry5<-mclustModel(subset(d.crabs,select=c(FL:BD)),mclst5,G=4,modelNames="EEV")
vc<-smry5$parameters$variance$sigma[,,1]
mn<-smry5$parameters$mean[,1]
y1<-f.vc.ellipse(vc,mn)
vc<-smry5$parameters$variance$sigma[,,2]
mn<-smry5$parameters$mean[,2]
y2<-f.vc.ellipse(vc,mn)
vc<-smry5$parameters$variance$sigma[,,3]
mn<-smry5$parameters$mean[,3]
y3<-f.vc.ellipse(vc,mn)
vc<-smry5$parameters$variance$sigma[,,4]
mn<-smry5$parameters$mean[,4]
y4<-f.vc.ellipse(vc,mn)
mclst5.model<-cbind(matrix(NA,500*4,3),rbind(y1,y2,y3,y4))
colnames(mclst5.model)<-c("species","sex","index","FL","RW","CL","CW","BD")
d.crabs.model<-rbind(d.crabs,mclst5.model)
gd<-ggobi(d.crabs.model)[1]
glyph_color(gd) = c(rep(4,50), rep(1,50), rep(9,50), rep(6,50), rep(8,2000))