k<-function(x){     #kernel
return((1/sqrt(2*pi))*exp(-x^2/2))
}

nw<-function(t,x,y,h){#Nadaraya-Watson
    txh<-outer(t,x,"-")/h     #Da todas las restas t-x
    ktxh<-k(txh)              #Hace k de la matriz
    L<-ktxh/apply(ktxh,1,sum) #da una matriz donde las filas suman uno
                        #divido cada fila por cada fila 
    return(L%*%y)}
    
m<-function(x){
return(dnorm(x,2,1)+dnorm(x,4,0.5)+dnorm(x,6,0.1)+dnorm(x,8,0.05))
    }


##función de regresión
x=(1:100)/10
y=m(x)
plot(x,y)
##media del estimador de NW
mbarra=nw(x,x,m(x),0.27)
plot(x,y,type="l")
points(x,mbarra,type="l",col="red")

##muestra
y=m(x)+rnorm(100,mean=0,sd=0.2)
mn=nw(x,x,y,0.27)
points(x,mn,type="l",col="blue")


##muestra aleatoria
x_muestra<-seq(1,100,by=0.05)
y_muestra=m(x_muestra)+rnorm(length(x_muestra),mean=0,sd=0.2)
matrizNW<-function(x,h){
    dnorm(outer(x,x,"-")/h)/apply(dnorm(outer(x,x,"-")/h),1,sum)}
y_muestra<-y_muestra[order(x_muestra)]
y1<-y_muestra[-1]
y2<-y_muestra[-length(y_muestra)]
sigma<-(1/(2*(length(y_muestra)-1)))*sum((y1-y2)^2)

##para la media
alpha=0.05
banda_inf=nw(x_muestra,x_muestra,y_muestra,0.27)-qnorm(1-(alpha/2))*sqrt(sum(matrizNW(x_muestra,0.27)^2)*sigma)
banda_sup=nw(x_muestra,x_muestra,y_muestra,0.27)+qnorm(1-(alpha/2))*sqrt(sum(matrizNW(x_muestra,0.27)^2)*sigma)
lines(x_muestra,banda_inf,col="brown")
lines(x_muestra,banda_sup,col="brown")
