## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width = 6,
  fig.height = 4
)

## ----data---------------------------------------------------------------------
library(psvr)
library(ggplot2)

# Target: highway fuel economy (all > 0)
y_all <- mpg$hwy

# Predictors: engine displacement, model year, cylinder count.
#
# `cty` is deliberately EXCLUDED. City and highway fuel economy are two
# measurements of the same property of the same vehicle (they correlate at
# 0.96), so predicting one from the other is leakage rather than modelling.
X_raw <- as.matrix(mpg[, c("displ", "year", "cyl")])

stopifnot(all(y_all > 0))
cat("N =", nrow(X_raw), "  p =", ncol(X_raw),
    "  y range: [", min(y_all), ",", max(y_all), "]\n")

## ----split--------------------------------------------------------------------
set.seed(42)
n      <- nrow(X_raw)
tr_idx <- sample(n, floor(0.7 * n))

X_raw_tr <- X_raw[tr_idx, ];  y_tr <- y_all[tr_idx]
X_raw_te <- X_raw[-tr_idx, ]; y_te <- y_all[-tr_idx]

# Standardise: centre and scale by training mean/sd
col_mean <- colMeans(X_raw_tr)
col_sd   <- apply(X_raw_tr, 2, sd)
X_tr <- scale(X_raw_tr, center = col_mean, scale = col_sd)
X_te <- scale(X_raw_te, center = col_mean, scale = col_sd)

## ----metrics------------------------------------------------------------------
mape  <- function(y, yhat) mean(abs(y - yhat) / y) * 100
rmspe <- function(y, yhat) sqrt(mean(((y - yhat) / y)^2)) * 100
r2    <- function(y, yhat) 1 - sum((y - yhat)^2) / sum((y - mean(y))^2)

## ----baseline-----------------------------------------------------------------
lm_df_tr <- as.data.frame(X_tr)
lm_df_te <- as.data.frame(X_te)
lm_fit   <- lm(y_tr ~ ., data = lm_df_tr)
lm_pred  <- predict(lm_fit, newdata = lm_df_te)

cat(sprintf("Linear regression — MAPE: %.2f%%  RMSPE: %.2f%%  R²: %.4f\n",
            mape(y_te, lm_pred), rmspe(y_te, lm_pred), r2(y_te, lm_pred)))

## ----mape-svr-----------------------------------------------------------------
# make_kernel() returns a closure K(xi, xj) = exp(-||xi - xj||^2 / (2 sigma^2)).
# sigma is a LENGTH in the units of the preprocessed feature space, so it has
# to be set on that scale -- sigma_heuristic() reads it off the data instead of
# guessing. See "Hyperparameter search ranges" below.
K <- make_kernel("rbf", sigma = sigma_heuristic(X_tr))

# C = 10: per-sample box bound |beta_k| <= 100*C/y_k; eps = 1: tube width (% of y_k)
fit_ep <- psvr_mape(X_tr, y_tr, kernel = K, C = 10, eps = 1)
pred_ep <- predict(fit_ep, X_te)

cat(sprintf("ε-SVR MAPE    — MAPE: %.2f%%  RMSPE: %.2f%%  R²: %.4f\n",
            mape(y_te, pred_ep), rmspe(y_te, pred_ep), r2(y_te, pred_ep)))
cat(sprintf("Support vectors: %d / %d\n", length(fit_ep$beta), fit_ep$n_train))
print(fit_ep)

## ----mape-plot, echo = FALSE--------------------------------------------------
lim <- range(c(y_te, pred_ep))
data.frame(actual = y_te, predicted = pred_ep) |>
  ggplot(aes(actual, predicted)) +
  geom_point(colour = "#d7191c", alpha = 0.65, size = 1.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") +
  coord_equal(xlim = lim, ylim = lim) +
  labs(x = "Actual hwy (mpg)", y = "Predicted hwy (mpg)",
       title = "Model 1: \u03b5-SVR with MAPE")

## ----mape-coef----------------------------------------------------------------
cf_ep <- coef(fit_ep)
# alpha, alpha_star: length-N dual variables (paired); the pre-pruning
#               solution.  Useful as a warm-start across CV folds; for
#               prediction use `beta` instead.
# beta:         beta_k = alpha_k - alpha_k* for each SUPPORT VECTOR only
#               (non-zero only for training points outside the
#               percentage-error ε-tube — sparse)
# b:            bias / intercept term
# support_data: training rows corresponding to support vectors only
cat(sprintf("b = %.4f  |  beta range: [%.4f, %.4f]\n",
            cf_ep$b, min(cf_ep$beta), max(cf_ep$beta)))

## ----rmspe-lssvr--------------------------------------------------------------
# gamma = 5000: regularisation; larger gamma -> smaller Y_Gamma diagonal -> tighter fit.
# This is roughly var(y_tr) * N, the scale cost_psvr_ls_data() computes -- and
# already five times the ceiling of the registered `cost` default. See below.
fit_ls <- psvr_rmspe(X_tr, y_tr, kernel = K, gamma = 5000)
pred_ls <- predict(fit_ls, X_te)

cat(sprintf("LS-SVR RMSPE  — MAPE: %.2f%%  RMSPE: %.2f%%  R²: %.4f\n",
            mape(y_te, pred_ls), rmspe(y_te, pred_ls), r2(y_te, pred_ls)))
print(fit_ls)

## ----rmspe-plot, echo = FALSE-------------------------------------------------
lim <- range(c(y_te, pred_ls))
data.frame(actual = y_te, predicted = pred_ls) |>
  ggplot(aes(actual, predicted)) +
  geom_point(colour = "#2c7bb6", alpha = 0.65, size = 1.8) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") +
  coord_equal(xlim = lim, ylim = lim) +
  labs(x = "Actual hwy (mpg)", y = "Predicted hwy (mpg)",
       title = "Model 3: LS-SVR with RMSPE")

## ----rmspe-coef---------------------------------------------------------------
cf_ls <- coef(fit_ls)
# alpha:        N dual variables; weight each training point's kernel
#               contribution in f(x) = sum_k alpha_k K(x_k, x) + b
#               (all N points, no sparsity)
# b:            bias / intercept term
# support_data: all N training inputs stored for prediction
cat(sprintf("b = %.4f  |  alpha range: [%.4f, %.4f]\n",
            cf_ls$b, min(cf_ls$alpha), max(cf_ls$alpha)))

## ----comparison---------------------------------------------------------------
results <- data.frame(
  Model = c("Linear regression", "\u03b5-SVR MAPE (Model 1)",
            "LS-SVR RMSPE (Model 3)"),
  MAPE  = c(mape(y_te,  lm_pred),
            mape(y_te,  pred_ep),
            mape(y_te,  pred_ls)),
  RMSPE = c(rmspe(y_te, lm_pred),
            rmspe(y_te, pred_ep),
            rmspe(y_te, pred_ls)),
  R2    = c(r2(y_te,    lm_pred),
            r2(y_te,    pred_ep),
            r2(y_te,    pred_ls))
)
results[, 2:4] <- round(results[, 2:4], 2)
knitr::kable(results, col.names = c("Model", "MAPE (%)", "RMSPE (%)", "R²"),
             align = "lrrr",
             caption = paste("Test-set performance on ggplot2::mpg",
                             "(70/30 split, RBF kernel, single run,",
                             "untuned hyperparameters)."))

## ----ranges-defaults, message = FALSE-----------------------------------------
library(parsnip)
library(tune)

show_ranges <- function(ps) {
  for (i in seq_len(nrow(ps))) {
    ob <- ps$object[[i]]
    cat(sprintf("  %-10s %s\n", ps$id[i],
      if (inherits(ob, "quant_param"))
        sprintf("[%s] on the %s scale",
                paste(signif(unlist(ob$range), 4), collapse = ", "),
                if (is.null(ob$trans)) "identity" else ob$trans$name)
      else sprintf("{%s}", paste(ob$values, collapse = ", "))))
  }
  invisible(ps)
}

spec_mape <- psvr_mape_rbf(cost = tune(), margin = tune(),
                           rbf_sigma = tune(), sym_type = tune()) |>
  set_engine("psvr")

extract_parameter_set_dials(spec_mape) |> show_ranges()

## ----ranges-mape--------------------------------------------------------------
extract_parameter_set_dials(spec_mape) |>
  update(
    cost      = cost_psvr(),                      # [-2, 10] log2 — fine for C
    margin    = margin_percentage(),              # 1-20% of each target
    rbf_sigma = rbf_sigma_psvr_data(X_tr),        # data-driven; see below
    sym_type  = sym_type_param(c("even", "odd"))  # drops "none"; see below
  ) |>
  show_ranges()

## ----ranges-rmspe-------------------------------------------------------------
spec_ls <- psvr_rmspe_rbf(cost = tune(), rbf_sigma = tune(),
                          sym_type = tune()) |>
  set_engine("psvr")

extract_parameter_set_dials(spec_ls) |>
  update(
    cost      = cost_psvr_ls_data(y_tr),
    rbf_sigma = rbf_sigma_psvr_data(X_tr),
    sym_type  = sym_type_param(c("even", "odd"))
  ) |>
  show_ranges()

## ----ranges-sigma-------------------------------------------------------------
sigma_heuristic(X_tr)   # the geometric centre of the range printed above,
                        # and the value the fits at the top of this page used