##
##  Moving plot showing how a locally weighted regression traces out as the weight
##  distribution slides over a scatterplot.
##
##  jamesHonaker  20/2/07
##

rm( list=ls() )     					# Clear all variables in memory

sim=TRUE

### Simulate Dataset ###############################

form<-function(x){
  yhat<- -0.5*x -1.5*x^2 + 1*x^3
  return(yhat)
}

if (sim){
  n=100
  upper.bound<-2.5
  lower.bound<- -1.3
  x<-runif(n)*(upper.bound - lower.bound) + lower.bound

  y<-form(x)
  y<-y + rnorm(n)*1.4

}else{
  x<-as.vector(read.table(file="x.txt"))
  y<-as.vector(read.table(file="y.txt"))
}

  dataset<-as.data.frame(cbind(x,y))


ncuts<-500
increment<-1
xlb<-seq(min(x)-increment/2,max(x)-increment/2,length=ncuts)
xseq<-seq(from=min(x), to=max(x), length=100)
yhat<-form(xseq)

all.ylim<-c(-5,5)
myTitle<-"LOESS Model"


for(i in 1:ncuts){


plot(x,y,xlim=c(lower.bound,upper.bound),ylim=all.ylim,ylab="",xlab="",main=myTitle)
par(new=TRUE)
plot(xseq,yhat,xlim=c(lower.bound,upper.bound),ylim=all.ylim,type="l",col="blue",ylab="",xlab="",lty=2,main=myTitle,lwd=2)

  abline(v=xlb[i])
  abline(v=xlb[i]+increment)

  flag<-x<xlb[i]+increment & x>xlb[i]

  w<- 1- abs((xlb[i]+increment/2 - x)/(increment/2))
  output<-lm(y[flag] ~ x[flag], weights=w[flag])

  segments(lower.bound, output$coefficients[1] + output$coefficients[2]*lower.bound, upper.bound,output$coefficients[1] + output$coefficients[2]*upper.bound,col="blue",lwd=2)

  segments(xlb[i], output$coefficients[1] + output$coefficients[2]*xlb[i], xlb[i]+increment, output$coefficients[1] + output$coefficients[2]*(xlb[i]+increment),col="red",lwd=2)
  points(xlb[i]+increment/2,all.ylim[1]+1,col="green",pch=24)
  points(xlb[i]+increment/2,all.ylim[2]-1,col="green",pch=25)
  segments(xlb[i]+increment/2,all.ylim[1]+1,xlb[i]+increment/2,all.ylim[2]-1,col="green")

  yhat.now<-output$coefficients[1] + output$coefficients[2]*(xlb[i]+increment/2)
  trace<-matrix(c(xlb[i]+increment/2,yhat.now),1,2)
  if(i==1){
    historyw<-trace
  }else{
    historyw<-rbind(historyw,trace)
  }

  myl<-min(all.ylim)
 polygon(x=c(xlb[i],xlb[i]+increment/2,xlb[i]+increment),y=c(myl-0.3,myl+1-0.3,myl-0.3),col="grey")

  points(historyw,pch=20)
  par(ask=FALSE)
}