#ANALISIS DISCRIMINANTE

####################
#Ejemplo1: B?sico para entender el an?lisis discriminante
####################

#dos especies de mosquitos al que se le mide longitud trompa y longitud de ala

X=matrix(c(138,164,140,170,124,172,136,174,138,182,148,182,154,182,138,190,156,
208,114,178,120,186,118,196,130,196,126,200,128,200),15,2,byrow=T)

Y=c(rep(0,9),rep(1,6))


X=cbind(X,Y)
X

colnames(X)=c("Trompa","Ala","Cat")
X_1=X[1:9,-3]
X_1

S_1=var(X_1)

X_2=X[10:15,-3]
S_2=var(X_2)

S=1/13*(8*S_1+5*S_2)
solve(S)

d=colMeans(X_1)-colMeans(X_2)

w=solve(S)%*%as.matrix(d)

#z=w'x=0.582x_1-0.382x_2

D2=t(as.matrix(d))%*%solve(S)%*%as.matrix(d)
D=sqrt(D2)


##########################################
# Comparaci?n con la funci?n lda de MASS #
##########################################
library(MASS)
l=lda(Cat~.,data=as.data.frame(X))

l$scaling

w

#Encontr? lo mismo:
l$scaling/w

q=qda(Cat~.,data=as.data.frame(X))
q
ls(q)


###########################################
#Clasificamos
dis=t(w)%*%((colMeans(X_1)+colMeans(X_2))/2)

#Asignaremos la clase 0 si z>dis y la clase 1 si z<dis

	#Por ejemplo
	x=as.matrix(c(118,196))
	z=t(w)%*%x
	z
	x=as.matrix(c(138,164))
	z=t(w)%*%x
	z
	x=as.matrix(c(130,190))
	z=t(w)%*%x
	z


##########################
# Ejemplo 2: De J.M Marin, universidad Carlos III, Madrid.
##########################
 
#Esos datos correponden a dos tipos raciales diferentes en los que se
#practicaron diferentes medidas de longitudes achuras de craneo y cara.
#analisis discriminante sobre dos tipos raciales.


source("chap7tibetskull.dat")


##############################
#LDA a pie                   #
##############################
attach(Tibet)
	Tibet1 <- Tibet[Type==1,c(1:5)]
	Tibet2 <- Tibet[Type==2,c(1:5)]

	n1=nrow(Tibet1)
	n2=nrow(Tibet2)
	n=n1+n2
	p=ncol(Tibet1)

# Vector de medias de las dos poblaciones:
	mean1=apply(Tibet1, 2, mean)
	mean2=apply(Tibet2, 2, mean)

# Estimaci?n de la matriz de varianzas covarianzas de toda la poblaci?n:
	S=((n1-1)*var(Tibet1) + (n2-1)*var(Tibet2))/(n-2)
	inv.S=solve(S)

# Eje discrimante:
	w=inv.S %*% (mean1-mean2)


##########################################
# Comparaci?n con la funci?n lda de MASS #
##########################################

# Usamos la funci?n lda de la libreria MASS
library(MASS)
	dis=lda(Type ~ Length + Breadth + Height + Fheight + Fbreadth, 
           data=Tibet, prior=c(0.5,0.5))

# comparo con el vector w encontrado en la parte anterior. 

	# Vector de la funci?n discriminante:
	dis$scaling

# En realidad encontr? un vector colineal:
	dis$scaling / w



###############################################
# Clasificaci?n de una nueva observaci?n a pie#
###############################################


newdata <- rbind(c(171, 140.5, 127, 69.5, 137),
                 c(179, 132, 140, 72, 138.5))

# hallo el punto de corte en la direcci?n discriminante para ambos grupos:
	lda.1 <- mean1 %*% w
	lda.1

	lda.2 <- mean2 %*% w
	lda.2

	puntodecorte=(lda.1 + lda.2)/2
	puntodecorte
	


# Regla de clasificaci?n: ser? de tipo 1 si el score encontrado es mayor que 
#       -30.46349, y de tipo 2 en otro caso.

# calculo el score para newdata
	newdata %*% w

	# la primera observaci?n es de tipo 1
	# la segunda observaci?n es de tipo 2


#######################################################
# Clasificaci?n de una nueva observaci?n usando MASS #
#######################################################

# necesitamos que newdata sea un data.frame
	dimnames(newdata)=list(NULL, c("Length","Breadth","Height",
                                  "Fheight","Fbreadth"))
	newdata.frame=data.frame(newdata)

# usamos la funci?n predict:
	pred=predict(dis, newdata=newdata.frame)
	pred

# clase que se predice:
	pred$class


################################ 
# Performance de la predicci?n #
################################ 

# predicci?n sobre la muestra usando la lda encontrada
	group=predict(dis, method="plug-in")$class
	group

# Quiero comparar con la verdadera etiqueta:
	table(group, Type)

# tasa de error de clasificaci?n:
	(3+3)/n


# un mejor metodo es usando validaci?n cruzada
	predicciones= array(NA, n)
	for (i in 1:n){
   	dat <- Tibet[-i,]
   	dis <- lda(Type ~ Length+Breadth+Height+Fheight+Fbreadth, 
              data=dat, prior=c(0.5,0.5))
   	predicciones[i] <- predict(dis, newdata=Tibet[i,c(1:5)])$class
	}
predicciones

table(predicciones, Type)

# error de clasificaci?n
(6+5)/n



########################
#Ejemplo3: Iris de Fisher
#######################

library(MASS)
attach(iris)

data=data.frame(iris)
#Funci?n discriminante de MASS
	ss=lda(Species~.,data=iris)
#Histogramas de cada grupo
	plot(ss, dimen =1)
#Dibujo en 2D
	plot(ss, dimen =2, col=c("green","red","blue")[unclass(data[,5])], 
	main="LDA plot of Iris data",xlab="Score sobre primer eje discriminante",
	ylab="Score sobre segundo eje discriminante")

#Matriz de confusi?n

	newdata.frame=data.frame(data)
	pred=predict(ss, newdata=newdata.frame)
	pred
	table(newdata.frame[,5],pred$class)
	
	#Con validaci?n cruzada
	ss2=lda(Species~.,data=iris)
	pred=predict(ss2, newdata=newdata.frame)
	pred
	table(newdata.frame[,5],pred$class)

	#Quadratic discriminant analysis
	ss3=qda(Species~.,data=iris)
	pred=predict(ss3, newdata=newdata.frame)
	pred
	table(newdata.frame[,5],pred$class)