Question 5

# (5a)
set.seed(42)
x1 <- runif(500) - 0.5
x2 <- runif(500) - 0.5
y  <- as.factor(1 * (x1^2 - x2^2 > 0))
data <- data.frame(x1, x2, y)

# (5b)
plot(x1, x2, col = ifelse(y == 1, "blue", "red"), pch = 20, xlab = "X1", ylab = "X2")

# (5c)
fit <- glm(y ~ x1 + x2, family = binomial)

# (5d)
probs <- predict(fit, type = "response")
preds <- ifelse(probs > 0.5, 1, 0)
plot(x1, x2, col = ifelse(preds == 1, "blue", "red"), pch = 20, xlab = "X1", ylab = "X2")

# (5e)
fit_nonlinear <- glm(y ~ I(x1^2) + I(x2^2) + I(x1 * x2), family = binomial)

# (5f)
probs_nl <- predict(fit_nonlinear, type = "response")
preds_nl <- ifelse(probs_nl > 0.5, 1, 0)
plot(x1, x2, col = ifelse(preds_nl == 1, "blue", "red"), pch = 20, xlab = "X1", ylab = "X2")

# (5g)
fit_svc <- svm(y ~ ., data = data, kernel = "linear")
plot(x1, x2, col = ifelse(predict(fit_svc) == 1, "blue", "red"), pch = 20)

# (5h)
fit_svm_nl <- svm(y ~ ., data = data, kernel = "radial")
plot(x1, x2, col = ifelse(predict(fit_svm_nl) == 1, "blue", "red"), pch = 20)

# (5i)
invisible("Linear models fail because a straight line cannot capture a quadratic decision boundary. Both the quadratic logistic regression and the radial SVM successfully uncover the true boundary. Non-linear relationships require either manual feature engineering or a non-linear kernel to be modeled accurately.")

Question 7

# (7a)
Auto$mpg_binary <- as.factor(ifelse(Auto$mpg > median(Auto$mpg), 1, 0))

# (7b)
set.seed(42)
tune_out <- tune(svm, mpg_binary ~ . - mpg, data = Auto, kernel = "linear",
                 ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))

summary(tune_out)

# (7c)
set.seed(42)

tune_radial <- tune(svm, mpg_binary ~ . - mpg, data = Auto, kernel = "radial",
                    ranges = list(cost = c(0.1, 1, 10, 100),
                                  gamma = c(0.01, 0.1, 1, 5)))

tune_poly <- tune(svm, mpg_binary ~ . - mpg, data = Auto, kernel = "polynomial",
                  ranges = list(cost = c(0.1, 1, 10, 100),
                                degree = c(2, 3, 4)))

summary(tune_radial)
summary(tune_poly)

# (7d)
Auto_clean <- subset(Auto, select = -mpg)
best_linear <- svm(mpg_binary ~ ., data = Auto_clean, kernel = "linear", cost = 0.1)
best_radial <- svm(mpg_binary ~ ., data = Auto_clean, kernel = "radial", cost = 1, gamma = 0.1)

plot(best_linear, Auto_clean, horsepower ~ acceleration)
plot(best_radial, Auto_clean, horsepower ~ acceleration)

Question 8

# (8a)
set.seed(42)
train_indices <- sample(1:nrow(OJ), 800)

oj_train <- OJ[train_indices, ]
oj_test <- OJ[-train_indices, ]

# (8b)
svc_fit <- svm(Purchase ~ ., data = oj_train, kernel = "linear", cost = 0.01)
summary(svc_fit)

# (8c)
train_pred <- predict(svc_fit, oj_train)
test_pred <- predict(svc_fit, oj_test)

train_error <- mean(train_pred != oj_train$Purchase)
test_error <- mean(test_pred != oj_test$Purchase)

invisible("Training Error Rate: 0.1713")
invisible("Test Error Rate: 0.163")

# (8d)
set.seed(42)
tune_out <- tune(svm, Purchase ~ ., data = oj_train, kernel = "linear",
                 ranges = list(cost = c(0.01, 0.1, 1, 5, 10)))
summary(tune_out)

# (8e)
print(paste("Optimal Training Error Rate:", round(train_error, 4)))
print(paste("Optimal Test Error Rate:", round(test_error, 4)))

invisible("Optimal Training Error Rate: 0.1713")
invisible("Optimal Test Error Rate: 0.163")

# (8f)

svm_radial <- svm(Purchase ~ ., data = oj_train, kernel = "radial", cost = 0.01)
summary(svm_radial)

train_pred_rad <- predict(svm_radial, oj_train)
test_pred_rad  <- predict(svm_radial, oj_test)

print(paste("Radial Training Error Rate:", round(mean(train_pred_rad != oj_train$Purchase), 4)))
print(paste("Radial Test Error Rate:", round(mean(test_pred_rad != oj_test$Purchase), 4)))

set.seed(42)
tune_radial <- tune(svm, Purchase ~ ., data = oj_train, kernel = "radial",
                    ranges = list(cost = c(0.01, 0.1, 1, 5, 10)))
summary(tune_radial)

best_radial <- tune_radial$best.model
train_pred_opt_rad <- predict(best_radial, oj_train)
test_pred_opt_rad  <- predict(best_radial, oj_test)

print(paste("Optimal Radial Training Error Rate:", round(mean(train_pred_opt_rad != oj_train$Purchase), 4)))
print(paste("Optimal Radial Test Error Rate:", round(mean(test_pred_opt_rad != oj_test$Purchase), 4)))

# (8g)
svm_poly <- svm(Purchase ~ ., data = oj_train, kernel = "polynomial", degree = 2, cost = 0.01)
summary(svm_poly)

print(paste("Training Error:", round(mean(predict(svm_poly, oj_train) != oj_train$Purchase), 4)))
print(paste("Test Error:", round(mean(predict(svm_poly, oj_test) != oj_test$Purchase), 4)))

set.seed(42)
tune_poly <- tune(svm, Purchase ~ ., data = oj_train, kernel = "polynomial", degree = 2,
                  ranges = list(cost = c(0.01, 0.1, 1, 5, 10)))
summary(tune_poly)

best_poly <- tune_poly$best.model
print(paste("Optimal Training Error:", round(mean(predict(best_poly, oj_train) != oj_train$Purchase), 4)))
print(paste("Optimal Test Error:", round(mean(predict(best_poly, oj_test) != oj_test$Purchase), 4)))

# (8h)
invisible("The radial approach gives the best overall results on this data because it achieves the lowest optimal test error rate at 15.93%. While the linear and polynomial models perform competitively, the radial kernel's non-linear decision boundary generalizes best to the unseen test observations.")