rm(list=ls())
library(rpart)

### point (a) ###

myAdaBoost <- function(X, y, B)
{
  # 1: initialization
  N <- NROW(X)
  f <- list()
  w <- rep(1/N, N)
  b  <- matrix(NA, ncol = N, nrow = B)
  err <- alpha <- NULL
  
  # 2: iterations
  for (m in 1:B)
  {
    f[[m]] <- rpart(as.factor(y) ~ ., data = as.data.frame(X), weights = w, control = rpart.control(maxdepth = 1))
    b[m, ] <- 2*(as.numeric(predict(f[[m]], newdata = as.data.frame(X), type = 'class')) - 1.5)
    # note that rpart gives 1 and 2 instead of -1, 1
    # multiply x2 because the estimare - 1.5 gives -0.5 or 0.5
    err[m] <- w %*% (b[m,] != y) / sum(w)
    alpha[m] <- log((1 - err[m]) / err[m])
#    w <- w * exp(alpha[m] * as.numeric(b[m, ] != y))
    w <- w * exp(alpha[m] * as.numeric(b[m, ] != y) - alpha[m]/2) # edited to take into account constant
    w <- w/sum(w)
    print(m)
  }

  # 3: final classifier
  g <- rep(0, N)
  for (m in 1:B) g <- g + alpha[m] * b[m, ]
  list(prediction = ifelse(g > 0, 1, -1), alpha = alpha, b = b, f = f)
}



### point (b) ###

# generate the data
set.seed(1)
X <- matrix(rnorm(20000), ncol = 10)
y <- apply(X^2, 1, sum)
y[y <= qchisq(p = 0.5, df = 10) ] <- -1
y[y > qchisq(p = 0.5, df = 10) ] <- 1
table(y)
B <- 400
N <- length(y)

# fit the algorithm
model <- myAdaBoost(X, y, B)

# wrt boosting iterations
train.error<-NULL
for (i in 1:B)
{
  g <- rep(0, N)
  for (m in 1:i) g <- g + model$alpha[m] * model$b[m, ]
  train.error[i] <- sum((g > 0) != (y > 0)) / N
}
plot(1:B, train.error, type = 'l')

# test error
set.seed(2)
X.test <- matrix(rnorm(100000), ncol = 10)
y.test <- apply(X.test^2, 1, sum)
y.test[y.test <= qchisq(p = 0.5, df = 10) ] <- -1
y.test[y.test > qchisq(p = 0.5, df = 10) ] <- 1
table(y.test)
N.test <- length(y.test)

test.error <- NULL
b.test  <- matrix(NA, ncol = N.test, nrow = B)
g.test <- rep(0, N)
for (m in 1:B)
{
  b.test[m,] <- 2*(as.numeric(predict(model$f[[m]], newdata = as.data.frame(X.test), type = 'class')) - 1.5)
  g.test <- g.test + model$alpha[m] * b.test[m, ]
  test.error[m] <- sum((g.test > 0) != (y.test > 0)) / N.test
}  
plot(1:B, test.error, type = 'l')

stump <- rpart(as.factor(y) ~ ., data = as.data.frame(X), control = rpart.control(maxdepth = 1))
g_stump <- as.numeric(predict(stump, newdata = as.data.frame(X.test), type = 'class')) - 1.5
test.error_stump <- sum((g_stump > 0) != (y.test > 0)) / N.test
abline(h = test.error_stump, lty = 3)
text(200, 0.47, 'single stump')

single.tree <- rpart(as.factor(y) ~ ., data = as.data.frame(X), control = rpart.control(maxdepth = 30))
g_single.tree <- as.numeric(predict(single.tree, newdata = as.data.frame(X.test), type = 'class')) - 1.5
test.error_single.tree <- sum((g_single.tree > 0) != (y.test > 0)) / N.test
abline(h = test.error_single.tree, lty = 3)
text(200, 0.3, 'single tree')



### point(c) ###
B <- 1000
# rerun with B = 1000
which.min(test.error)



### point (d) ###

set.seed(1)
X <- matrix(rnorm(20000), ncol = 10)
y <- apply(X^2, 1, sum)
y[y <= 12 ] <- -1
y[y > 12 ] <- 1
table(y)
N <- length(y)
B <- 1000

model.2 <- myAdaBoost(X, y, B)

# wrt boosting iterations
train.error <- NULL
for (i in 1:B)
{
  g <- rep(0, N)
  for (m in 1:i) g <- g + model.2$alpha[m] * model.2$b[m, ]
  train.error[i] <- sum((g > 0) != (y > 0)) / N
}
plot(1:B, train.error, type = 'l')


# test error
set.seed(2)
X.test <- matrix(rnorm(100000), ncol = 10)
y.test <- apply(X.test^2, 1, sum)
y.test[y.test <= 12 ] <- -1
y.test[y.test > 12 ] <- 1
table(y.test)
N.test <- length(y.test)

test.error <- NULL
b.test  <- matrix(NA, ncol = N.test, nrow = B)
g.test <- rep(0, N.test)
for (m in 1:B)
{
  b.test[m,] <- 2*(as.numeric(predict(model.2$f[[m]], newdata = as.data.frame(X.test), type = 'class')) - 1.5)
  g.test <- g.test + model.2$alpha[m] * b.test[m, ]
  test.error[m] <- sum((g.test > 0) != (y.test > 0)) / N.test
}  
plot(1:B, test.error, type = 'l')

stump <- rpart(as.factor(y) ~ ., data = as.data.frame(X), control = rpart.control(maxdepth = 1))
g_stump <- as.numeric(predict(stump, newdata = as.data.frame(X.test), type = 'class')) - 1.5
test.error_stump <- sum((g_stump > 0) != (y.test > 0)) / N.test
abline(h = test.error_stump, lty = 3)
text(500, 0.27, 'single stump')

single.tree <- rpart(as.factor(y) ~ ., data = as.data.frame(X), control = rpart.control(maxdepth = 30))
g_single.tree <- as.numeric(predict(single.tree, newdata = as.data.frame(X.test), type = 'class')) - 1.5
test.error_single.tree <- sum((g_single.tree > 0) != (y.test > 0)) / N.test
abline(h = test.error_single.tree, lty = 3)
text(500, 0.22, 'single tree')

which.min(test.error)