######################################################
#            REGRESIÓN Y ANÁLISIS DE DATOS           #
######################################################

data.frame(airquality)
names(airquality)
airquality$Ozone
ind<-which(is.na(airquality$Ozone))
airquality$Temp[-ind]

REGRESION<-lm(airquality$Ozone[-ind]~airquality$Temp[-ind])
plot(REGRESION)
plot(airquality$Temp[-ind],airquality$Ozone[-ind])
abline(lm(airquality$Ozone[-ind]~airquality$Temp[-ind]))

######################################################
# 1.- Analizar la distribución univariante           #
######################################################
par(mfrow=c(1,2))
plot(ecdf(airquality$Ozone),main="Ozono")
plot(ecdf(airquality$Temp),main="Temperatura")

par(mfrow=c(1,2))
plot(density(airquality$Ozone),main="Ozono")
plot(density(airquality$Temp),main="Temperatura")

shapiro.test(airquality$Ozone)
shapiro.test(airquality$Temp)


length(airquality$Ozone)
length(airquality$Temp)

#--Serie de tiempo
Ozone<-ts(airquality$Ozone,frequency=30)
windows()
plot(Ozone, type="o", ylab="Ozono")

Temperatura<-ts(airquality$Temp,frequency=30)
windows()
plot(Temperatura, type="o", ylab="Temperatura")

summary(airquality$Ozone)
summary(airquality$Temp)


#######################################################
# 2.-Estimador local lineal                           #
#######################################################

K<-function(mas){
return(1/sqrt(2*pi)*exp(-mas^2/2))
}


#S<-function(t,mas,bandwidth,r){
#-- t vector de puntos donde se va a evaluar la función.
#-- mas valores de la variable emasplicativa
#-- y valores de la variable respuesta    
#-- r exponente

#    return(apply(K(outer(mas,t,"-")/bandwidth)%*%outer(mas,t,"-")%*%outer(mas,t,"-"),1,sum))
#}

EstLL<-function(t,mas,y,bandwidth){
#-- Esta función me devuelve la estimación por el estimador noparamétrico de la
#-- regresión del estimador local-lineal
#-- t vector de puntos donde se va a evaluar la función.
#-- mas valores de la variable emasplicativa
#-- y valores de la variable respuesta
#-- bandwidth ventana de suavizado
    S1<-K(outer(mas,t,"-")/bandwidth)%*%outer(mas,t,"-")
    S2<-K(outer(mas,t,"-")/bandwidth)%*%outer(mas,t,"-")%*%outer(mas,t,"-")
    return((K(outer(mas,t,"-")/bandwidth)%*%(S2-(outer(mas,t,"-")%*%S1))/apply(K(outer(mas,t,"-")/bandwidth)%*%(S2-(outer(mas,t,"-")%*%S1)),1,sum))%*%y)
}



####################################################
#    Validación cruzada: Estimador lineal_local    #
####################################################
library(MASS)

CV<- function(bandwidth,mas,y){
        #Lii<-K(0)*S(mas,mas,bandwidth,2)/apply(K(outer(mas,mas,"-")/bandwidth)%*%S(mas,mas,bandwidth,2)-outer(mas,mas,"-")%*%S(mas,t,bandwidth,1),1,sum)
        S1<-K(outer(mas,mas,"-")/bandwidth)%*%outer(mas,mas,"-")
        S2<-K(outer(mas,mas,"-")/bandwidth)%*%outer(mas,mas,"-")%*%outer(mas,mas,"-")
        return(apply((outer(y,EstLL(mas,mas,y,bandwidth),"-")%*%ginv(1-(K(0)*S2/apply(K(outer(mas,mas,"-")/bandwidth)%*%(S2-(outer(mas,mas,"-")%*%S1)),1,sum)),tol = sqrt(.Machine$double.eps)))^2,1,mean))
}


cross_validation_bandwidth<-optimize(CV,interval=c(0,max(dist(mas))),mas=mas,y=y)$minimum


#######################################################
#  Analizar resultados                                #
#######################################################

windows()
plot(airquality$Temp,airquality$Ozone)
Local-lineal<-EstLL(airquality$Temp,airquality$Temp,airquality$Ozone,bandwidth_CV))


########################################################
#     Cálculo de nueva ventana                         #
########################################################
library(KernSmooth)
dpill(airquality$Temp[-ind],airquality$Ozone[-ind])



########################################################
#     Nadaraya-Watson                                  #
########################################################
nw<-function(t,mas,y,bandwidth){
#-- Esta función me devuelve la estimación por el estimador noparamétrico de la
#-- regresión de Nadaraya-Watson
#-- t vector de puntos donde se va a evaluar la función.
#-- mas valores de la variable emasplicativa
#-- y valores de la variable respuesta
#-- bandwidth ventana de suavizado

return((K(outer(t,mas,"-")/bandwidth)/apply(K(outer(t,mas,"-")/bandwidth),1,sum))%*%y)

}



#####################################################
#    Regresión lineal simple                        #
#####################################################
windows()
plot(airquality$Temp[-ind],airquality$Ozone[-ind])
abline(lm(Ozone~Temp))
abline(lm(Ozone[-ind]~bs(Temp[-ind],degree=1),data=airquality))
summary(lm(Ozone[-ind]~bs(Temp[-ind],degree=1),data=airquality))
