背景重述

本文是ESL: 12.3 支持向量機和核中表12.2的重現過程。具體問題如下：
在兩個類別中產生100個觀測值。第一類有4個標準正態獨立特征\(X_1,X_2,X_3,X_4\)。第二類也有四個標準正態獨立特征，但是條件為\(9\le \sum X_j^2\le 16\)。這是個相對簡單的問題。同時考慮第二個更難的問題，用6個標準高斯噪聲特征作為增廣特征。

生成數據

## #####################################
## generate dataset
## 
## `No Noise Features`: num_noise = 0
## `Six Noise Features`: num_noise = 6
## #####################################
genXY <- 
 
function(n = 100, num_noise = 0)
{
  ## class 1
  m1 = matrix(rnorm(n*(4+num_noise)), ncol = 4 + num_noise)
  ## class 2
  m2 = matrix(nrow = n, ncol = 4 + num_noise)
  for (i in 1:n) {
    while (TRUE) {
      m2[i, ] = rnorm(4 + num_noise)
      tmp = sum(m2[i, 1:4]^2)
      if(tmp >= 9 & tmp <= 
 
16)
        break
    }
  }
  X = rbind(m1, m2)
  Y = rep(c(1, 2), each = n)
  return(data.frame(X = X, Y = as.factor(Y)))
}

模型訓練

1. SVM直接調用e1071包中的svm函數
2. BRUTO和MARS都是調用mda包，且由於兩者都是用於回歸，所以轉換為分類時，是比較擬合值與類別標簽的距離，劃分到越靠近的那一類
3. 原書中提到實驗中MARS不限定階數，但實際編程時，設置階數為10

交叉驗證選擇合適的\(C\)

我分兩步進行選擇：

1. 粗選：在較大範圍內尋找最優的\(C\)
2. 細分：在上一步選取的最優值附近進行細分

註意避免最優值取在邊界值。以SVM/poly5為例進行說明，其他類似

## SVM/poly5
set.seed(123)
poly5 = tune.svm(Y~., data = dat, kernel = "polynomial", degree = 5, cost = 2^(-4:8))
summary(poly5)

此時選取的最優\(C\)為32，進一步細化

set.seed(1234)
poly5 = tune.svm(Y~., data = dat, kernel = "polynomial", degree = 5, cost = seq(16, 64, by = 2))
summary(poly5)

所以\(C\)取28。

類似地，得到其它方法的最優\(C\)，比如某次實驗結果如下：

當然，實際中我們並不需要重新設置參數來訓練模型，因為tune.svm()的返回結果就包含了最優模型，直接調用，比如poly5$best.model

計算測試誤差

predict.mars2 <- function(model, newdata)
{
  pred = predict(model, newdata)
  ifelse(pred < 1.5, 1, 2)
}

calcErr <- function(model, n = 1000, nrep = 50, num_noise = 0, method = "SVM")
{
  err = sapply(1:nrep, function(i){
    dat = genXY(n, num_noise = num_noise)
    datX = dat[, -ncol(dat)]
    datY = dat[, ncol(dat)]
    if (method == "SVM")
      pred = predict(model, newdata = datX)
    else if (method == "MARS")
      pred = predict.mars2(model, newdata = datX)
    else if (method == "BRUTO")
      pred = predict.mars2(model, newdata = as.matrix(datX))
    sum(pred != datY)/(2*n) # Attention!! The total number of observations is 2n, not n
  })
  return(list(TestErr = mean(err),
              SE = sd(err)))
}

值得說明的是，對於BRUTO和MARS，因為程序是將其視為回歸模型處理的，需要進一步轉換為類別標簽。因為程序中類別用1和2編號，所以判斷擬合值是否大於1.5，大於則劃為第二類，否則第一類。

結果

將之與表12.2進行比較，可以看出各個方法的誤差率及標準差的相對大小都比較一致。

貝葉斯誤差率

對於類別1，
\[ \sum X_j^2\sim \chi^2(4) \]
對於類別2，
\[ \sum X_j^2\sim \frac{\chi^2(4)I(9\le\chi^2(4)\le 16)}{\int_9^{16} f(t)dt} \]
其中\(f(t)\)是\(\chi^2(4)\)的密度函數。

於是貝葉斯誤差率為

\[ \frac{1}{2}\int_{9}^{16}f(t)dt\approx 0.029 \]

完整代碼可以參見skin-of-the-orange.R

本文永久鏈接：模擬：Tab. 12.2

用一個簡單的例子比較SVM，MARS以及BRUTO（R語言）