# code based on files "Code regarding section 4.2", "Code regarding section 4.3-4.4",
# "Code regarding section 4.5" and "Code regarding section 4.6"
# (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa)

library(sm)

auto <- read.table("http://azzalini.stat.unipd.it/Book-DM/auto.dat", header = TRUE)
attach(auto)

y <- city.distance
x <- engine.size

# figure 4.1
add.window <- function(x, y, h, x.eval) {
  # function from "script rc_alter" connect to the book
  #    Bowman & Azzalini (1997) "Applied smoothing techniques 
  #    for data analysis", Oxford University Press.
  polygon(rep(c(x.eval - 2 * h, x.eval + 2 * h), rep(2,2)),
          c(range(y), rev(range(y))), col = 5, border = FALSE)
  lines(rep(x.eval, 2), range(y), lty = 2)
  points(x, y, pch = 16, cex = 1.5)
  xseq <- seq(x.eval - 3 * h, x.eval + 3 * h, length = 20)
  kernel <- dnorm(xseq, x.eval, h)
  kernel <- min(y) + 5 * kernel / max(kernel)
  lines(xseq, kernel, lty = 2)
}
par(mfrow = c(1, 1))
x.eval <- 3
hvec <- c(0.5, 0.5, 0.15, 1)
for (i in 1:4) {
  plot(x, y, pch = 16, cex = 1.5,
       xlab = "Engine size (L)", ylab = "City distance (km/L)")
  if(i>1)  {
    add.window(x, y, hvec[i], x.eval)
    sm.regression(x, y,   h = hvec[i], add = TRUE, ngrid = 250, col = 2, lwd = 2)
    text(4.5, 18,  paste("lambda =", as.character(hvec[i])))
  }
}


par(mfrow = c(1,1))
# figure 4.4a
y <- city.distance
x1 <- engine.size
x2 <- curb.weight
xlab <- "Engine size (L)"
ylab <- "Curb weight (kg)"
zlab <- "City distance (km/L)"

h <- c(0.5, 150)
sm.options(ngrid = 50)
a <- sm.regression(cbind(x1, x2), y, h = h, display = "none",
                   options = list(xlab = xlab, ylab = ylab, zlab = zlab))
persp(a$eval.points[ , 1], a$eval.points[ , 2], a$estimate,
      xlab = xlab, ylab = ylab, zlab = zlab, cex = 1, 
      theta = 120, phi = 20, ticktype = "detailed")

# figure 4.4b
contour(a$eval.points[ , 1], a$eval.points[ , 2], a$estimate,
        xlab = xlab, ylab = ylab)
points(x1, x2, col = 2,  cex = 2, pch = 16)
z <- chull(x1, x2)
z <- c(z, z[1])
lines(x1[z], x2[z], lty = 3, col = 3)

detach(auto)