#==============  Carga de los datos ========================================================================

rm(list=ls(all=TRUE))

jardin <- read.csv("E:\\Máster\\Asignaturas\\Análisis Exploratorio de Datos\\Profesor\\Examen 2008\\jardin.txt", header = TRUE, sep = ",")
names(jardin)

jardin=as.matrix(jardin)
class(jardin)
jardin

#plot(datos,col = "blue")

summary(jardin)
library(MASS)



# =============  Puesta a punto  de los datos =============================================================

#==Cada clase es  una variable.
#==El salario de el sera la clase 1, representa la primera columna de la matriz
#==El salario de ella sera la clase 2, representa la segunda columna de la matriz
x1 =as.matrix(jardin[ ,1])  
x2 =as.matrix(jardin[ ,2]) 
mean (x1)
#==Generando una nueva matriz con los datos a clasificar
x=cbind(x1,x2)
class(x)
#x

#==============  Clasificando los datos ==================================================================

# Asignemos la codificacion a las clases {0,1} (0,1 porque asi es mi respuesta la variable jardin)
# 0 representa la clase 1 (sueldo el) y 1 representa la clase 2 (sueldo ella).
y = as.matrix(jardin[ ,4])



#==============  Graficamos los datos =====================================================================

windows(6,5.5,xpos=0,ypos=0)
# pch es el código del tipo de símbolo
# col es el color que se aplica al símbolo
plot(x, pch=y+2, col = y+3, xlab="SalarioEl", ylab = "SalarioElla")  
# NO se debe cerrar la grafica resultante
 page(y)

#==============  Generando las muestras de entrenamiento y prueba ==========================================

# Extraemos una muestra aleatoria simple de tamaño 825
isample = sample(1:825) 
# se podrían sacar dos isamples de 1 a 100 para representar bien las dos clases en la muestra


  # Extraemos una muestra para entrenamiento de tamaño 625
  itrain = isample[1:625]
  xtrain = x[itrain]
  ytrain = y[itrain]


  # El resto de los datos queda como muestra de test 200
  xtest = x[-itrain]
  ytest = y[-itrain]


# Densidades condicionadas
#p(x/Y=0)  
fx0 = 0.5 * dmvnorm(x1, mean(x1))
fx1 = 0.5 * dmvnorm(x2, mean(x2))
 fx = 0.5 * fx0 + 0.5 * fx1 

 py1x = 0.5 * fx1 / fx





#==============  Contruimos la cuadricula   ====================================================================

# Generamos la secuencia de valores en cada eje
rx1 = seq(0,8,0.1)
rx1
rx2 = seq(0,6,0.1)
rx2

# Ponemos los puntos (pares) en fila en una matriz con dos columnas
# cogiéndolos en un orden tal que varía primero la x2 y luego la x1.
# Es decir, los cogemos en este orden:
#           (0.0,0.0), (0.0,0.1), (0.0,0.2), (0.0,0.3)....(0.0,6.0)
#           (0.1,0.0), (0.1,0.1), (0.1,0.2), (0.1,0.3)....(0.1,6.0)
#           ....
#           ....
#           (8.0,0.0), (8.0,0.1), (8.0,0.2), (8.0,0.3)....(8.0,6.0)

# El primer rep usa como patrón de repetición rep(length(rx2),length(rx1)) 
# que sirve para generar una estructura de length(rx2) huecos de tamaño length(rx1)
# que acojeran las repeticiones del vector completo rx1.
# El segundo rep simplemente repite el vector rx2 completo tantas veces
# como la longitud length(rx1)

xgrid = cbind(rep(rx1,rep(length(rx2),length(rx1))),rep(rx2,length(rx1)))
xgrid 



# ======================= REGRESIÓN LOGÍSTICA =============================================================

# Utilizamos la función glm del paquete stats que viene con R para realizar la regresión logística


#=====  Primero preparamos el data frame que necesita la función glm como argumento  =========================

#Longitud de la matriz x y la de y de entrenamiento y prueba
ngrid = length(xgrid[ ,1])
ngrid 
ntrain = length(ytrain)
ntrain
ntest = length(ytest)
ntest

#creo el data frame con las variables
xy.df = as.data.frame(cbind(x,y))
names(xy.df) = c("x1","x2","y")
xy.df

#ceo el dataframe con la cuadricula
xygrid.df = data.frame(cbind(xgrid,rep(0,ngrid)))
names(xygrid.df) = names(xy.df)
xygrid.df

#junto los dos data frame
xy.all.df = rbind(xy.df,xygrid.df)
xy.all.df 

# ========  Realizamos la regresión logística con el paquete glm que viene en el paquete stats que viene con R


    # Transformamos la respuesta a códigos {0,1} NO esto no se hace porque ya la tengo en cero uno
    #xy01.all.df = xy.all.df
    #xy01.all.df$y = (xy.all.df$y + 1)/2
    #logistic_glm = glm(y~x1+x2, binomial(link = "logit"), data=xy01.all.df, subset = itrain)


logistic_glm = glm(y~x1+x2, binomial(link = "logit"), data=xy.all.df, subset = itrain)
summary(logistic_glm)


    #Graficas del modelo glm, grafica de normalidad, atipicos d. cook... No los necesito
    #windows(6,5.5)
    #plot(logistic_glm)

    #phat_glm_train = predict(logistic_glm,xy01.all.df[itrain, ], type = "response") # con type="response" da las probabilidades
    #phat_glm_rest = predict(logistic_glm,xy01.all.df[-itrain, ], type = "response")



phat_glm_train = predict(logistic_glm,xy.all.df[itrain, ], type = "response") # con type="response" da las probabilidades

phat_glm_rest = predict(logistic_glm,xy.all.df[-itrain, ], type = "response")

phat_glm_test = phat_glm_rest[1:ntest]
phat_glm_grid = phat_glm_rest[(ntest+1):(ntest+ngrid)]



Error_glm_train = sum((2*(phat_glm_train>1/2)-1)!=ytrain)/ntrain
Error_glm_test = sum((2*(phat_glm_test>1/2)-1)!=ytest)/ntest
Error_glm_test



# matrix(x,nx1,nx2) rellena por columnas,
# pero queremos llenar por filas para que quede como outer
# genero una matriz traspuesta (nfilas y ncolumnas intercambiadas) y traspongo después
nfilas = length(rx1)
ncolumnas = length(rx2)


Prob_glm_grid <- t(matrix(phat_glm_grid,ncolumnas,nfilas))
Prob_glm_grid



# ==============  Graficamos la frontera de la logística===================================
dev.set(2)
contour(rx1,rx2,Prob_glm_grid, levels=0.5, lwd = 2, col = "blue", add = TRUE)

# Graficamos la estimación de las probabilidades a posteriori para la clase 1
# realizada por la regresión logística
windows(6,5.5)
persp(x=rx1,y=rx2,Prob_glm_grid,theta=-135,phi=50,  xlab="x1",ylab="x2",zlab="Probabilidad a Posteriori")




#============ Intentando el modelo CART =====================================

# ABRIR EL PAQUETE 
library (rpart)


mcart <- rpart(y~x1+x2, data=xy.all.df, subset = itrain)
mcart2 <- rpart(y~x1+x2, data=xy.all.df, subset = itrain,parms=list(prior=c(.65,.35), split='information'))
mcart3 <- rpart(y~x1+x2, data=xy.all.df, subset = itrain, control=rpart.control(cp=.05))
par(mfrow=c(1,2), xpd=NA) # otherwise on some devices the text is clipped
plot(mcart)
text(mcart, use.n=TRUE)
plot(mcart)
text(mcart2, use.n=TRUE)

mcart
mcart3 
mcart2


#phat_cart_train = predict(mcart3,xy.all.df[itrain, ], type = "response") # con type="response" da las probabilidades
#phat_cart_rest = predict(mcart3,xy.all.df[-itrain, ], type = "response")
#phat_cart_test = phat_cart_rest[1:ntest]
#phat_cart_grid = phat_cart_rest[(ntest+1):(ntest+ngrid)]


#Error_cart_train = sum((2*(phat_cart_train>1/2)-1)!=ytrain)/ntrain
#Error_cart_test = sum((2*(phat_cart_test>1/2)-1)!=ytest)/ntest
#Error_cart_test