# code based on files " Code regarding the first part of section 5.6",
# "Code regarding section 4.8", and "Code regarding section 5.7"
# (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa)

# figure 4.16
source("f_true.R")
x <- x.250
y <- f.true250
# par(mfrow = c(2,2))
plot(x, y, type = "l", ylab = "f(x)", lwd = 2)
#
plot(x, y, type = "l", lty = 2, ylab = "f(x)", lwd = 2)
abline(v = 0.8, lty = 3, col = 4, lwd = 2)
low <- (x < 0.8)
m1 <- mean(y[low])
m2 <- mean(y[!low])
lines(x[low], rep(m1, length = sum(low)), lwd = 2, col = 2)
lines(x[!low], rep(m2, length = sum(!low)), lwd = 2, col = 2)
#
plot(x, y, type = "l", lty = 2, ylab = "f(x)", lwd = 2)
abline(v = 0.8, lty = 3, col = 4, lwd = 2)
abline(v = 0.6, lty = 3, col = 4, lwd = 2)
low1 <- (x < 0.6)
low2 <- (x < 0.8) & (!low1)
m11 <- mean(y[low1])
m12 <- mean(y[low2])
lines(x[low1], rep(m11, length = sum(low1)), lwd = 2, col = 2)
lines(x[low2], rep(m12, length = sum(low2)), lwd = 2, col = 2)
lines(x[!low], rep(m2, length = sum(!low)), lwd = 2, col = 2)
#
plot(x, y, type = "l", lty = 2, ylab = "f(x)", lwd = 2)
abline(v = 0.8, lty = 3, col = 4, lwd = 2)
abline(v = 0.6, lty = 3, col = 4, lwd = 2)
abline(v = 1.8, lty = 3, col = 4, lwd = 2)
hi <- !low
hi1 <- hi & (x < 1.8)
hi2 <- hi & (x > 1.8)
m21 <- mean(y[hi1])
m22 <- mean(y[hi2])
lines(x[low1], rep(m11, length = sum(low1)), lwd = 2, col = 2)
lines(x[low2], rep(m12, length = sum(low2)), lwd = 2, col = 2)
lines(x[hi1], rep(m21, length = sum(hi1)), lwd = 2, col = 2)
lines(x[hi2], rep(m22, length = sum(hi2)), lwd = 2, col = 2)
# par(mfrow=c(1,1))

# figure 4.17a
load("f2_true.save") # load matrix z
np <- 61
x <- y <- seq(0, 3, length = np)
ind <- seq(1, np, by = 2)
x0 <- x[ind]
y0 <- y[ind]
z0 <- z[ind, ind]
persp(x0, y0, z0, theta = 150, phi = 20, r = 2.5, xlab = "x1", ylab = "x2", zlab = "y")

library(tree)
# figure 4.17b
nx <- length(x)
ny <- length(y)
xoy <- cbind(rep(x, ny), as.vector(matrix(y, nx, ny, byrow = TRUE)))
X0 <- matrix(xoy, nx * ny, 2, byrow = FALSE)
x1<- X0[ , 1]
x2 <- X0[ , 2]
zz <- matrix(z, nx * ny, 1, byrow = FALSE)
t1 <- tree(zz ~ x1 + x2)
#
zp <- matrix(predict(t1), nrow = nx, ncol = ny)
persp(x, y, zp, theta = 150, phi = 20, r = 2.5, xlab = "x1", ylab = "x2", zlab = "y")

# figure 4.18a
plot(t1, type = "uniform")
text(t1, digits = 2, pretty = 3)

# figure 4.18b
partition.tree(t1)

#------------------------------------------
dataset <- read.table("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", header = TRUE)
x <- rep(dataset$x,2)
y <- c(dataset$y.yesterday, dataset$y.tomorrow)

# figure 4.19
plot(x, y, type = "n")
points(dataset$x, dataset$y.yesterday, pch = 16, col = 2)
points(dataset$x, dataset$y.tomorrow, pch = 17, col = 3)
legend('bottomright', legend = c("Yesterday", "Tomorrow"), pch = 16:17, col = 2:3)

x0 <- dataset$x
y0 <- dataset$y.yesterday
t3 <- tree(y0 ~ x0, control = tree.control(nobs = 30, mindev = 0.001, minsize = 2))
# figure 4.21a
plot(t3)
par(mar = c(3.5, 3.5, 2, 1) + 0.1)
p3 <- prune.tree(t3, newdata = data.frame(x0 = x0, y0 = dataset$y.tomorrow))
# figure 4.21b
plot(p3)

p4 <- prune.tree(t3, best = 4)
# figure 4.21c
plot(p4)
text(p4)
# figure 4.21d
x0 <- seq(0.5, 3, length = 5000)
plot(x, y, type = "n")
points(dataset$x, dataset$y.yesterday, pch = 16, col = 2)
points(dataset$x, dataset$y.tomorrow, pch = 17, col = 3)
legend(2, 0.45, legend = c("Yesterday", "Tomorrow"), pch = 16:17, col = 2:3)
lines(x0, predict(p4, newdata = data.frame(x = x0)), col = 4)