library(pdist)
n<-200
X<-matrix(10*runif(n),ncol=1)
y<-sin(X[,1])#+rnorm(n,sd=.2)
#plot(X[,1],y)
#xnew<-3

#span<-1

my_loess<-function(xnew,X,y,span=.75){
  #xnew is a vector with length=ncol(X)
  #nn=number of nearest neighbors to consider
  nn=ceiling(length(y)*span)
  d<-sqrt(colSums((xnew-t(X))^2))
  rk<-rank(d,ties.method="random")
  ii<- rk<=nn
  d[ii]<-d[ii]/max(d[ii])
  w<-0*d
  w[ii]<-(1-abs(d[ii])^3)^3
  X<-cbind(1,X)
  c(1,xnew)%*%solve(crossprod(X,w*X),crossprod(X,w*y))
}

my_loess_vec<-function(Xnew,X,y,span=.75){
  apply(Xnew,1,my_loess,X,y,span)
}


Xnew<-matrix(seq(from=-10,to=20,length.out=100),ncol=1)
ynew<-my_loess_vec(Xnew,X,y,span=1)
plot(X,y,xlim=c(-5,15))
lines(Xnew,ynew)

fit<-loess(y~X,span=1,degree=1)
ydefault<-predict(fit,Xnew)
lines(Xnew,ydefault,col="red")

my_gp<-function(Xnew,X,y){
  Xnew<-cbind(1,Xnew)
  X<-cbind(1,X)
  D<-as.matrix(dist(X))
  M<-max(D)
  K<-(1-(D/M)^3)^3
  Dnew<-as.matrix(pdist(Xnew,X))
  Knew<-(1-(Dnew/M)^3)^3
  Knew%*%solve(K,y)
}

ynew2<-my_gp(Xnew,X,y)
lines(Xnew,ynew2,col="blue")
legend("bottomleft",c("loess_default","loess_custom","GP"),lty=rep(1,3),col=c("black","red","blue"))

#possible extension of tricube kernel to real domain
tricube<-function(d){
  u<-2/(1+exp(-d))-1 #map reals to (-1,1)
  70/81*(1-abs(u)^3)^3
}

curve(tricube,from=-10,to=10)