library(dplyr)
library(ggplot2)
library(plotly)
library(gridExtra)

En este laboratorio se quiere aprender un clasificador que pueda determinar si una flor, de acuerdo a las medidas de pétalos y sépalos, pertenece o no a una determinada familia. Por ello trabajaremos con un conjunto de datos (dataset) iris.data que ha sido generado por un experto en el tema. Para ello utilizaremos un modelo de regresión lineal.

Funcion que lee N datapoints del dataset iris.data y lo devuelve

read_dataset<-function(){
  dataset.location <- "./data/iris.data"
  XY <- read.csv(dataset.location,header=TRUE,stringsAsFactors=F)
  return(XY)
}

Función que conviernte el dataset a una matriz numerica

pre_process <- function()
{  
  XY <- read_dataset()
  
  #### REMOVE BLOCK CODE
  
  XY$Y <- as.character(XY$Y)
  XY <- XY %>% mutate(Y=ifelse(Y=="Iris-virginica",1,-1))
  XY <- XY %>%mutate(id=row_number()) #%>% select(id,X1,X2,X3,X4,Y) # Mantenemos un id para luego buscar por clase
  XY <- apply(as.matrix(XY),2,as.numeric) # requerido para utilizar el operador %*%
  XY
  #### END
  return(XY)
}

Función para realizar el aprendizaje del discriminador

Recibe una matriz numerica con los datos preprocesados y devuelve una lista list(XEntrenamiento, YEntrenamiento, XTesteo, YTesteo, w) donde:

  1. XEntrenamiento: es la matriz X correspondiente a datos de entrenamiento
  2. YEntrenamiento: es la matriz Y correspondiente a datos de entrenamiento
  3. XTesteo: es la matriz X correspondiente a datos de testeo
  4. YTesteo: es la matriz Y correspondiente a datos de testeo
  5. w: es la matriz de pesos w obtenida en (4)
learn <- function(data,N)
{  
  #set.seed(N)
  finEntr  <- N * 0.7
  XYDesordenada <- data[sample(nrow(data),N),]
  XEntrenamiento <- NULL
  YEntrenamiento <- NULL
  XTesteo <- NULL
  YTesteo <- NULL

  #### REMOVE BLOCK CODE
  XDesordenada <- XYDesordenada[,1:4]
  YDesordenada <- XYDesordenada[,5]
  IdDesordenada <- XYDesordenada[,6]
  XEntrenamiento <- XDesordenada[1:finEntr,]
  YEntrenamiento <- YDesordenada[1:finEntr]
  IdEntrenamiento <- IdDesordenada[1:finEntr]
  
  XTesteo <- XDesordenada[(finEntr+1):N,]
  YTesteo <- YDesordenada[(finEntr+1):N]
  IdTesteo <- IdDesordenada[(finEntr+1):N]
  
  # Se agrega uno para el BIAS
  XEntrenamiento <- cbind(matrix(1,ncol=1,nrow=nrow(XEntrenamiento)),XEntrenamiento)
  XTesteo <- cbind(matrix(1,ncol=1,nrow=nrow(XTesteo)),XTesteo)
  
  w <- solve(t(XEntrenamiento)%*%XEntrenamiento)%*%t(XEntrenamiento)%*%YEntrenamiento
  
  #### END
  ret <- list(XEntrenamiento=XEntrenamiento,
              YEntrenamiento=YEntrenamiento,
              XTesteo=XTesteo,
              YTesteo=YTesteo,
              IdEntrenamiento=IdEntrenamiento,
              IdTesteo=IdTesteo,
              w=w)
  return(ret)
}

Función para realizar la evaluacion del desempeño del discriminador

  1. recibe una matriz numerica con los datos preprocesados
  2. Evaluación el discriminador calculado su error
  3. Muestra por pantalla y guarda en el archivo outputs.raws los resultados (punto c)
  4. Descripción de campos:
  5. nombre de dataset
  6. N
  7. número de repetición
  8. error de clasificación
  9. Recall de las clases positivas (classes positivas que fueron correctamente detectadas como positivas)
  10. Recall de las clases negativas (classes negativas que fueron correctamente detectadas como negativas)
  11. Vector de pesos.

Version utilizando dyplr

evaluate <- function(data)
{  

  N <- c(50,100,150)
  salida <- c()
  for(n in N)
  {
    for(i in 1:10){
      exitos<-0;
      recall<-0;
      l <- learn(data,n)      
      #### REMOVE BLOCK CODE
      XTesteo <- l$XTesteo
      XTesteo <-cbind(XTesteo,Ypred=c(ifelse(XTesteo %*% l$w>0,1,-1)))
      XTesteo <-cbind(XTesteo,Y=l$YTesteo)
      XTesteo<-as.data.frame(XTesteo)
      XTesteo$Ypred<-as.factor(XTesteo$Ypred)
      XTesteo$Y<-as.factor(XTesteo$Y)
      error<-XTesteo %>% summarise(error=sum(Y!=Ypred)/n()) %>% select(error)
      posrecall<-XTesteo %>%  summarise(posrecall=sum(Y==Ypred& Y==1)/sum(Y==1)) %>% select(posrecall)
      negrecall<-XTesteo %>% summarise(negrecall=sum(Y==Ypred& Y==-1)/sum(Y==-1)) %>% select(negrecall)
      #### END
      print(l$w)
      salida<-rbind(salida,cbind("iris.data",n,i,error,posrecall,negrecall,t(l$w)))
    }
}
  ref<-c("nombre","N","rep","error","recall_pos","recall_neg","bias","W1","W2","W3","W4");
  colnames(salida)<-ref
  print(salida)
  write.csv(salida,file="output.csv",quote = F);
}

Esta sección presenta los resultados de los experimentops

Analisis de los resultados

  1. Distribucion del error por N
  2. Promedio del error por N
  3. Distribucion del recall negativo (classes positivas que fueron correctamente detectadas como positivas)
  4. Distribucion del recall positivo
  5. Analisis del error por clase
  6. Analisis de la influencia de las variables (pesos w)
csv.data <- read.csv("output.csv",header=T)
csv.data$N<-as.factor(csv.data$N)

error=ggplot(csv.data,aes(x=N,y=error))+
  geom_boxplot(colour='skyblue')+
  stat_summary(fun.y=mean, colour="red", geom="point", shape=19, size=2)+
  theme_classic()+
  ggtitle("Dist. del error ")


recall_pos=ggplot(csv.data,aes(x=N,y=recall_pos))+
  geom_boxplot(colour='skyblue')+
  stat_summary(fun.y=mean, colour="red", geom="point", shape=19, size=2)+
  theme_classic()+
  ggtitle("Dist. del recall_pos relativo")

recall_neg=ggplot(csv.data,aes(x=N,y=recall_neg))+
  geom_boxplot(colour='skyblue')+
  stat_summary(fun.y=mean, colour="red", geom="point", shape=19, size=2)+
  theme_classic()+
  ggtitle("Dist. del recall_neg relativo")

data_weights<-csv.data %>% select(N,bias,W1,W2,W3,W4) %>% reshape2::melt(id.vars="N",variable.name="weigths") 
grid.arrange(error,recall_neg,recall_pos,ncol=3)

Analisis de los pesos

Analisis del Bias

El bias, puede interpretarse como el Valor Esperado cuando las variables \(X1,X2,X3,X4\) toman valor \(0\). En este caso se observa que el Bias es negativo, lo que significa que el valor esperado sera negativo. Esto ultimo resulta razonable, ya las clase \(-1\) es mayoritaria (2/3) y por lo tanto se esperaria observarla con mayor frecuencia.

Analisis de los coeficientes \(W\)

Las variables \(X2\) \(X4\) correspondiente a “Sepal.Width” y “Petal.Width” en el dataset iris, presentan el mayor peso. Ya Las restantes variables tienen valores muy cercanos a \(0\). Se puede interpretar como el ancho “width”, es el atributo mas importante para construccion del modelo.

weights=ggplot(data_weights,aes(x=weigths,y=value))+
  geom_boxplot(colour='skyblue')+
  ylab("valor")+xlab("pesos")+
  theme_bw()+
  ggtitle("Dist. de los pesos")+
  facet_wrap(~N)
weights

names(iris)
## [1] "Sepal.Length" "Sepal.Width"  "Petal.Length" "Petal.Width" 
## [5] "Species"

Proyeccion en 2D de las predicciones de un caso con N=100

PCA (Principal Component Analisys)

(tomado textual de wikipedia) El PCA construye una transformación lineal que escoge un nuevo sistema de coordenadas para el conjunto original de datos en el cual la varianza de mayor tamaño del conjunto de datos es capturada en el primer eje (llamado el Primer Componente Principal), la segunda varianza más grande es el segundo eje, y así sucesivamente. Para construir esta transformación lineal debe construirse primero la matriz de covarianza o matriz de coeficientes de correlación. Debido a la simetría de esta matriz existe una base completa de vectores propios de la misma. La transformación que lleva de las antiguas coordenadas a las coordenadas de la nueva base es precisamente la transformación lineal necesaria para reducir la dimensionalidad de datos. Además las coordenadas en la nueva base dan la composición en factores subyacentes de los datos iniciales.

Una de las ventajas del PCA para reducir la dimensionalidad de un grupo de datos, es que retiene aquellas características del conjunto de datos que contribuyen más a su varianza, manteniendo un orden de bajo nivel de los componentes principales e ignorando los de alto nivel. El objetivo es que esos componentes de bajo orden a veces contienen el aspecto “más importante” de esa información.

l<-learn(pre_process(),100)
YObtenido <- l$XTesteo %*% l$w
YObtenido <- data.frame(Y=YObtenido) 
YObtenido<-YObtenido %>% mutate(Ypred=ifelse(Y>0,1,-1)) %>% select(Ypred)

pca<-prcomp(l$XTesteo[,2:(ncol(l$XTesteo))], center = TRUE, scale. = TRUE) 

iris_data<-read_dataset()
pca_labeled<-data.frame(pca$x,Y=l$YTesteo,Ypred=YObtenido,especie=iris_data[l$IdTesteo,5]
)

pca_labeled$Y<-as.factor(pca_labeled$Y)
pca_labeled$Ypred<-as.factor(pca_labeled$Ypred)
pca_labeled<-pca_labeled %>% mutate(asignacion=ifelse(Ypred!=Y,"Incorrecto","Correcto"))

pca_labeled$asignacion<-as.factor(pca_labeled$asignacion)


ggplot(pca_labeled,aes(x=pca_labeled$PC1,y=pca_labeled$PC2,col=Y,text=especie)) +
 # geom_point()+
  geom_point(aes(shape=asignacion),size=3)+
  ylab("PC1")+xlab("PC2")+
  theme_classic()+
scale_shape_manual(values=c(8,6))

Version dinamica que mantiene informacion original sobre las clases

ggplotly()