#---------------------------------------------------------------#
# un dataset ficticio, creaciónn de las matrices                     #

rm(list=ls())

                            
g <- function(x){s <- round(abs(sin(x)*1000 ))
                 set.seed(s)                 # randomize
                 # pero garantiza que g(constante) es una constante
                 x*(1- 0.5*rnorm(1)*sin(x))} # para añadir perturbación

set.seed(2); x1 <-  matrix(rnorm(10000))
set.seed(3); x2 <-  matrix(rnorm(10000))
y     <- matrix(2*x1-7*x2+1)
y     <- apply(y,1,g)
Z     <- data.frame(x1,x2,y)


#---------------------------------------------------------------#
# set entrenamiento y validación set                           #
                            
names(Z)                    # X = [x1,x2] ; y = [y]

require(caret) # nos da una partición elegante
set.seed(1) # iniciar generador de números aleatorios
index <- createDataPartition(y=Z$y, # extraer índice
                             p= 60/100, # que es una secuencia de
                             list=FALSE) # 6000 numss aleatorios

str(index) # comprobar que es correcto

Ze <- Z[ index,] # obtener 6000 filas
Zv <- Z[-index,] # extracciones de 4000 In complementarios

head(index) # verificaciones - primeros 5 índices
head(Ze) # primeras 5 filas en estos índices
dim (Ze) # tenemos 6000 líneas
head(Zv) # 5 primeras In en otros índices
dim(Zv) # tenemos 4000 lineas

y <- matrix(Z$y) # las y
ye <- matrix(Ze$y) # las y's de entrenamiento
yv <- matriz(Zv$y) # las y's de validación

dim(y) # vector 10000 coord
dim(ye) # vector 6000 coord
dim(yv) # vector 4000 coord

#---------------------------------------------------------------#
# visor de datos interactivo                                  #
library(iplots)

iplot(Ze$x1,Ze$x2)
ihist(Ze$y)


#---------------------------------------------------------------#
# test Breusch Pagan                                            #

library(lmtest) # paquete con pruebas para lm
bptest(y ~ x1 + x2, data = Ze) # la prueba
? bptest # ayuda de prueba


#------------------------------------------------- --------------#
# creación de modelo lineal en conjunto de entrenamiento #
# aplicación en conjunto de validación #

f <- lm(y ~ x1+x2,data=Ze) # creación del modelo y=f(X)
yv_ <- predict(f,newdata = Zv) # f(Xv)
yv_ <- matriz(yv_)

# almacenar nuestros resultados
# para uso futuro
prediccion_lm <- data.frame(Zv[,-3], # eliminar y
                            prediccion_y = yv_,
                            row.names = NULL)
tail(prediccion_lm)
head(prediccion_lm)

f_lm <- f # memorisación del modelo para
# comparación futura


require(ggplot2) # para tener qplot()
# visualización del resultado
qplot(yv,yv_, # comparación y e y predichas
      geom = c("point", "smooth"),
      method="lm",
      alfa = I(1/2)
) # aspecto satisfactorio

#---------------------------------------------------------------#
# profundizar en el estudio del resultado                           #

e <- yv - yv_ # vector de errores (= residuos)

qplot(yv_,e, # función de dispersión del error de yv
      geom = c("punto", "suave"),
      method="lm",
      alfa = I(1/2)
)
# porcentaje
d <- max(yv_)-min(yv_) # desviación máxima en la predicción
p <- e/d *100 # proporción de error

qplot(yv_,p, # estudio proporción error
      geom = c("point", "smooth"),
      method="lm",
      ylab = "porcentaje de error",
      alfa = I(1/2)
)


qplot(yv_, # densidad de predicciones
      geom = "density",
      método="lm",
      size = I(1.1), # grosor de línea
      color = I ("steelblue") # color de línea
)

#---------------------------------------------------------------#
# cálculo manual de los indicadores actuales                        #

library(matrixcalc)         # para hacer cálculos matricieles

n  <- length(yv)            # num elementos de yv
e  <- yv - yv_              # vector de errores (= residus)
n  <- length(yv)            # num elementos de yv

SSE  <- entrywise.norm( e, 2 )**2        # SSE
MSE  <- SSE / n    # MSE
RMSE <- sqrt(MSE)                        # RMSE
MAE <- entrywise.norm( e, 1 ) / n        # MAE
ME <- mean(e)                            # ME
NRMSE <- RMSE/(max(yv) - min(yv))        # NRMSE
CV_RMSE <- RMSE/mean(yv)                 # CV_RMSE
SRD <- sqrt(mean(e**2) - mean(e)**2)     # SRD

(ME**2 + SRD**2) - MSE                   # verificación

MEAN   <- mean(yv)                       # media yv
CENTER <- yv-MEAN                        # centrado yv
SST  <- entrywise.norm( CENTER, 2 )**2   # SST 
R2 <- 1- SSE/SST                         # R2                      
 
id1 <- data.frame(n,SSE,MSE,RMSE,MAE,ME,NRMSE,CV_RMSE,SRD,SST,R2)
id1


#---------------------------------------------------------------#
# trabajamos directamente en el modelo de predicción generado    #
summary(f)                # detalle del modelo

a <- f$coefficients[2]    # extracción de los coeff del modelo
b <- f$coefficients[3]
c <- f$coefficients[1]

f_ <- function(x1,x2){    # re creación de la función del 
         a*x1+b*x2+c      # modelo
         }

                                     # visualización del modelo
x1_ <- seq(min(x1), max(x1), by=0.2) # grid x1
x2_ <- seq(min(x2), max(x2), by=0.2) # grid x2
y_ <- outer(x1_, x2_, f_)            # un valor por parej

persp(x1_,                           # perspectiva 3D
      x2_,
      y_,
      phi=40,                        # ángulos presentación
      theta=60,
      col= "green",
      ticktype = "detailed",         # eje detalaldo
      shade=.0001,
      cex.lab=1.0 )                  # atención zoomer para
                                     # visualizarlo

#---------------------------------------------------------------#
paleta <- colorRampPaleta(c("red","yellow"))(100)
diag <- persp(x1_,                           # perspectiva 3D
              x2_,
              y_,
              phi=40,                        # ángulo presentación
              theta=60,
              col=paleta[round((y_+max(y_))/(max(y_)-min(y_))* 100)],          
              ticktype = "detailed",         # eje detallado
              shade=.0001,
              cex.lab=1.5 )             # atención zoomer para
                                        # visualizarlo

index_ <- createDataPartition(y=Z$y,      # extracción índice
                              p= 2/100,   # seguido de
                              list=FALSE) # 6000 num aleatorios

u <- as.vector(matrix(Z[index_,1]))
v <- as.vector(matrix(Z[index_,2]))
w <- as.vector(matrix(Z[index_,3]))

nube <- trans3d(u,                       # nube de puntos
                 v,
                 v,
                 pmat = diag)

points(nube,                            # inserción en diag
       pch = v + max(v)+1,      # en forma  de símbolos
       cex = v/3.5 + max(v/3),  # juega con el tamaño de los símbolos 
       col = v + max(v)+1       # juega con el color de los símbolos
       )

#---------------------------------------------------------------#
# líneas de nivel                                              #
contour(x1_,x2_,y_,
        ,method = "edge",
        xlab="x1",
        ylab="x2",
        col=round(abs(y_))%%5+1) # jugar con los colores

#---------------------------------------------------------------#
# prueba con el algoritmo randomForest #
# aplicación al conjunto de validación #

require(randomForest) # árboles de decisión del bosque
f <- randomForest(y ~ x1+x2,
                      data=Ze,
                      ntree = 30) # creación del modelo y=f(X)
yv_ <- predict(f,newdata = Zv) # f(X)
yv_ <- matrix(yv_)


# visualización del resultado
qplot(yv,yv_, # comparación y e y predichas
      geom = c("point", "smooth"),
      method="lm",
      alfa = I(1/2)
) # apariencia satisfactoria

#------------------------------------------------- --------------#
# estudio en profundidad del resultado #

e <- yv - yv_ # vector de errores (= residuos)

qplot(yv_,e, # función de dispersión del error de yv
      geom = c("point", "smooth"),
      method="lm",
      alfa = I(1/2)
)
# porcentaje
d <- max(yv_)-min(yv_) # desviación máxima en la predicción
p <- e/d # proporción de error

qplot(yv_,p, # estudio proporción error
      geom = c("point", "smooth"),
      method="lm",
      alfa = I(1/2)
)


qplot(yv_, # densidad de predicciones
      geom = "density",
      method="lm",
      alfa = I(1/2)
)


e <- yv - yv_ # vector de errores (= residuos)
n <- length(yv) # num elementos de yv

SSE  <- entrywise.norm( e, 2 )**2        # SSE
MSE  <- SSE / n    # MSE
RMSE <- sqrt(MSE)                        # RMSE
MAE <- entrywise.norm( e, 1 ) / n        # MAE
ME <- mean(e)                            # ME
NRMSE <- RMSE/(max(yv) - min(yv))        # NRMSE
CV_RMSE <- RMSE/mean(yv)                 # CV_RMSE
SRD <- sqrt(mean(e**2) - mean(e)**2)     # SRD

(ME**2 + SRD**2) - MSE                   # verificación

MEAN   <- mean(yv)                       # media yv
CENTER <- yv-MEAN                        # centrado yv
SST  <- entrywise.norm( CENTER, 2 )**2   # SST 
R2 <- 1- SSE/SST                         # R2                      
R2 
id2 <- data.frame(n,SSE,MSE,RMSE,MAE,ME,NRMSE,CV_RMSE,SRD,SST,R2)

#---------------------------------------------------------------#
# trabajamos directamente en el modelo de predicción generado    #
summary(f)                # detalle del modelo

names(Z)
f_ <- function(x1,x2){    # creación función predicción
      x1 <- matrix(c(x1))
      x2 <- matrix(c(x2))
      y<- matrix(c(1))
      linea_X <- data.frame(x1,x2,y)
      predict(f,newdata = linea_X )
      }

# visualización del modelo
x1_ <- seq(min(x1), max(x1), by=0.2) # grid x1
x2_ <- seq(min(x2), max(x2), by=0.2) # grid x2
y_ <- outer(x1_, x2_, f_)            # un valor por pareja

paleta <- colorRampPaleta(c("red","yellow"))(100)

diag <- persp(x1_,                           # perspectiva 3D
      x2_,
      y_,
      phi=40,                        # ángulos presentación
      theta=60,
      col=paleta[round((y_+max(y_)+10)/(max(y_)-min(y_))* 100)],
      ticktype = "detailed",         # eje detallado
      shade=.0001,
      cex.lab=1.5 )                  # atención zoomer para
# el visualizador
index_ <- createDataPartition(y=Z$y,      # extracción índice
                              p= 2/100,   # que es una sucesión de
                              list=FALSE) # 6000 num aleatorios

u <- as.vector(matrix(Z[index_,1]))
v <- as.vector(matrix(Z[index_,2]))
w <- as.vector(matrix(Z[index_,3]))

nube <- trans3d(u,
                 v,
                 v,
                 pmat = diag)

points(nube,                            # inserción en diag
       pch = v + max(v)+1,      # jugar sobre la forma de los símbolos
       cex = v/4 + max(v/3.5),  # jugar sobre el tamaño de los símbolos
       col = v + max(v)+1       # jugar con el color de los símbolos
)


#---------------------------------------------------------------#
# líneas de nivel                                              #
contour(x1_,x2_,y_,
        ,method = "edge",
        xlab="x1",
        ylab="x2",
        col=round(abs(y_))%%5+1) # jugar con los colores

#---------------------------------------------------------------#
# comparación de los indicadores
id <- rbind(id1,id2)
id

