## Source:   Chapter11.r (with modifications) in 
## http://www.stat.colostate.edu/computationalstatistics/CompStatR.zip


############################################################################
### EXAMPLE 11.7 EASY DATA (LOESS)
############################################################################
# easy = observed data
# x    = observed predictor data
# y    = observed response data
# k    = smoothing parameter
############################################################################
## NOTES
# The following uses the lowess function. The loess function may also be
# used, however it does has a different implementation and defaults.
############################################################################

## INITIAL VALUES
easy = read.table(file.choose(),header=T)
x    = easy$X
y    = easy$Y
n    = length(y)
k    = 30

## MAIN
s.hat = lowess(x=x,y=y,f=k/n)

## OUTPUT
s.hat      # LOESS OUTPUT

## OUTPUT PLOTS, Figure 11.9
s = function(x){(x^3) * sin((x+3.4)/2)}
x.plot = seq(min(x),max(x),length.out=1000)
y.plot = s(x.plot)
plot(x,y,xlab="Predictor",ylab="Response")
lines(x.plot,y.plot,lty=2)
lines(s.hat)

## without data points
plot(x,y,xlab="Predictor",ylab="Response", main="Loess Smoother", type="n")
lines(x.plot,y.plot,lty=2)
lines(s.hat)
legend("bottomright",c("True relation","k=30"),
       lty=c(2,1))


# running-line
TMA = function(k,y){
  n = length(y)
  S = matrix(0,n,n)
  b = (k-1)/2
  if(k>1){
    for(i in 1:b){
      S[i,1:(b+i)] = 1/(k-b+i-1)
      S[n-i+1,(n-b-i+1):n] = 1/(k-b+i-1)
    }
    for(i in (b+1):(n-b)){
      S[i,(i-b):(i+b)] = 1/k
    }}
  if(k==1){S = diag(1,n)}
  out = S%*%y
  return(out)
}
## output
s23  = TMA(23,y)
plot(x,y,xlab="Predictor",ylab="Response", main="Running-line Smoother", type="n")
lines(x.plot,y.plot,lty=2)
lines(x,s23,lty=1)
legend("bottomright",c("True relation","k=23"),
       lty=c(2,1))


## Example 11.7
x1    = c(x, 1,1,1)
y1    = c(y,-8,-8,-8)
n1    = length(y1)

## OUTPUT PLOTS, 
# LOESS
## MAIN
s1.hat = lowess(x=x1,y=y1,f=k/n1)

x.plot = seq(min(x),max(x),length.out=1000)
y.plot = s(x.plot)
plot(x,y,xlab="Predictor",ylab="Response")
lines(s1.hat,lty=1,lwd=2)
lines(s.hat,lty=2,lwd=2)

# running-line
temp=order(x1)
x2=x1[temp]
y2=y1[temp]
## output
s23  = TMA(23,y)
s23a  = TMA(23,y2)
plot(x,y,xlab="Predictor",ylab="Response")
lines(x,s23,lty=1,lwd=2)
lines(x2,s23a,lty=2,lwd=2)

## replot Figure 11.10
par(mfrow=c(1,2))
plot(x,y,xlab="Predictor",ylab="Response",type="n", main="Running-line smooth")
lines(x,s23,lty=3,lwd=2)
lines(x2,s23a,lty=1,lwd=2)
plot(x,y,xlab="Predictor",ylab="Response",type="n", main="Loess smooth")
lines(s.hat,lty=3,lwd=2)
lines(s1.hat,lty=1,lwd=2)

## replot Figure 11.10 with legend
par(mfrow=c(1,2))
plot(x,y,xlab="Predictor",ylab="Response",type="n", main="Running-line smooth")
lines(x,s23,lty=3,lwd=2)
lines(x2,s23a,lty=1,lwd=2)
legend("bottomright",c("Without outliers","With outliers"),
       lty=c(3,1), cex=0.7)
plot(x,y,xlab="Predictor",ylab="Response",type="n", main="Loess smooth")
lines(s.hat,lty=3,lwd=2)
lines(s1.hat,lty=1,lwd=2)
legend("bottomright",c("Without outliers","With outliers"),
       lty=c(3,1), cex=0.7)



############################################################################
### EXAMPLE 11.8 DIFFICULT DATA (SUPERSMOOTHER)
############################################################################
# tough       = observed data
# x           = observed predictor data
# y           = observed response data
# h           = smoothing parameters
# RLSMOOTH    = running-line smoother (for h odd > 1)
############################################################################
## NOTES
# This example uses the supersmoother implemented in R via the
# "supsmu" function. See the ?supsmu help page for more information.
############################################################################

## INITIAL VALUES
tough = read.table("toughsmooth.dat",header=T)
x     = tough$X
y     = tough$Y
n     = length(y)
h     = c(0.05*n, 0.2*n, 0.5*n)

## OUTPUT PLOTS
par(mfrow=c(1,1))
plot(x,y,xlab="Predictor",ylab="Response",main="Running-line Smooths")

## FUNCTION (running lines smoother).  The supersmoother is called
## using supsmu()

RLSMOOTH = function(k,y,x){
      if(k>1){
      n = length(y)
      s.hat = rep(0,n)
      b = (k-1)/2
      x.bar = mean(x[1:(b+1)])
      y.bar = mean(y[1:(b+1)])
      xy.bar = mean(x[1:(b+1)]*y[1:(b+1)])
      x2.bar = mean(x[1:(b+1)]^2)
      beta1 = (xy.bar - x.bar*y.bar)/(x2.bar - x.bar^2)
      s.hat[1] = y.bar + beta1*(x[1]-x.bar)
      for(i in 2:(b+1)){
            x.bar  = (x.bar*(i+b-1) + x[i+b])/(i+b)
            y.bar  = (y.bar*(i+b-1) + y[i+b])/(i+b)
            xy.bar = (xy.bar*(i+b-1) + x[i+b]*y[i+b])/(i+b)
            x2.bar = (x2.bar*(i+b-1) + (x[i+b]^2))/(i+b)
            beta1  = (xy.bar - x.bar*y.bar)/(x2.bar - x.bar^2)
            s.hat[i] = y.bar + beta1*(x[i]-x.bar)
      }
      for(i in (b+2):(n-b)){
            x.bar = x.bar + (x[i+b] - x[i-b-1])/k
            y.bar = y.bar + (y[i+b] - y[i-b-1])/k
            xy.bar = xy.bar + (x[i+b]*y[i+b] - x[i-b-1]*y[i-b-1])/k
            x2.bar = x2.bar + ((x[i+b]^2)-(x[i-b-1]^2))/k
            beta1 = (xy.bar - x.bar*y.bar)/(x2.bar - x.bar^2)
            s.hat[i] = y.bar + beta1*(x[i]-x.bar)
      }
      for(i in (n-b+1):n){
            x.bar  = (x.bar*(k-i+n-b+1) - x[i-b-1])/(k-i+n-b)
            y.bar  = (y.bar*(k-i+n-b+1) - y[i-b-1])/(k-i+n-b)
            xy.bar = (xy.bar*(k-i+n-b+1) - x[i-b-1]*y[i-b-1])/(k-i+n-b)
            x2.bar = (x2.bar*(k-i+n-b+1) - (x[i-b-1]^2))/(k-i+n-b)
            beta1  = (xy.bar - x.bar*y.bar)/(x2.bar - x.bar^2)
            s.hat[i] = y.bar + beta1*(x[i]-x.bar)
      }
      }
      if(k==1){s.hat = y}
      return(s.hat)
}

## MAIN
ssmooth.val  = supsmu(x,y)
rlsmooth.val = mapply(RLSMOOTH,h+1,MoreArgs = list(y,x))

## OUTPUT PLOTS
par(mfrow=c(1,2))
plot(x,y,xlab="Predictor",ylab="Response",main="Running-line Smooths")
lines(x,rlsmooth.val[,1],type="l")
lines(x,rlsmooth.val[,2],type="l",lty=2)
lines(x,rlsmooth.val[,3],type="l",lty=3)
legend("bottomright",c("k=0.05n","k=0.2n","k=0.5n"),lty=1:3, cex=0.5)

plot(x,y,xlab="Predictor",ylab="Response",main="Supersmooth")
lines(ssmooth.val)

## without data points
par(mfrow=c(1,2))
plot(x,y,xlab="Predictor",ylab="Response",main="Running-line Smooths", type="n")
lines(x,rlsmooth.val[,1],type="l")
lines(x,rlsmooth.val[,2],type="l",lty=2)
lines(x,rlsmooth.val[,3],type="l",lty=3)
legend("bottomright",c("k=0.05n","k=0.2n","k=0.5n"),lty=1:3, cex=0.5)

plot(x,y,xlab="Predictor",ylab="Response",main="Supersmooth", type="n")
lines(ssmooth.val)