setwd("~/Dropbox/Teaching/MVA/AAMA/computing/")

# install.packages("car")
# install.packages("klaR")
library(MASS)
library(car)

# iris data
# 4 dime, n = 150.
data(iris)
attach(iris) 
iris$Species <- factor(iris$Species)
scatterplotMatrix(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width | Species , diagonal="density",smooth=FALSE,reg.line=FALSE)
 

irisv = iris[51:150,]
irisv$Species <- factor(irisv$Species)
attach(irisv)
scatterplotMatrix(~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width | Species , diagonal="density",smooth=FALSE,reg.line=FALSE)

gpr <- prcomp(irisv[,1:4])
Z <- gpr$x
scatterplotMatrix(~ Z| Species , diagonal="density",smooth=FALSE,reg.line=FALSE)

scatterplot(PC1 ~ PC2 | Species ,Z, smooth=FALSE,reg.line=FALSE)
fit <- lda(Species ~ Z[,1:2]) 
fit # show results
abline(0,-fit$scaling[2]/fit$scaling[1],pch=5)
#treat data as matrix
z = as.matrix(Z)
fit <- qda(Species ~ Z[,1:2])
fit # show results
library(klaR)
partimat(Species ~ Z[,1:2],method="lda",prec=100)
partimat(Species ~ Z[,1:2],method="qda",prec=100)

## naive methods - Mean Difference
mu1 <-colMeans(Z[1:50,])
mu2 <-colMeans(Z[51:100,])
overalmu = (mu1+mu2)/2
scatterplot(PC1 ~ PC2 | Species ,Z, smooth=FALSE,reg.line=FALSE)
#abline(0,-fit$scaling[2]/fit$scaling[1],pch=5)
points(mu1[2],mu1[1],pch=17,cex=1)
points(mu2[2],mu2[1],pch=17,cex=1,col=2)
md = mu2-mu1
abline(0,-md[1]/md[2],pch=5)

## naive methods - Naive Bayes
Z1 = Z[,1]/sd(Z[,1])
Z2 = Z[,2]/sd(Z[,2])
scatterplot( Z2 ~ Z1| Species, smooth=FALSE,reg.line=FALSE)
mu1 <- c(mean(Z1[1:50]),mean(Z2[1:50]))
mu2 <- c(mean(Z1[51:100]),mean(Z2[51:100]))
points(mu2[2],mu2[1],pch=17,cex=1,col=1)
points(mu1[2],mu1[1],pch=17,cex=1,col=2)
md = mu2-mu1
abline(0,-md[1]/md[2],pch=5)



## 3-way classification
attach(iris)
gpr <- prcomp(iris[,1:4])
Z <- gpr$x
scatterplotMatrix(~ Z| Species , diagonal="density",smooth=FALSE,reg.line=FALSE)

scatterplot(PC1 ~ PC2 | Species ,Z, smooth=FALSE,reg.line=FALSE)

fit <- lda(as.matrix(Z[,1:2]),Species)
fit # show results
prediction<-predict(fit, Z[,1:2])
prediction$class
fit <- qda(Species ~ Z[,1:2])
fit # show results
library(klaR)
partimat(Species ~ Z[,1:2],method="lda",prec=100)
partimat(Species ~ Z[,1:2],method="qda",prec=100)