用一个简单的例子比较SVM，MARS以及BRUTO（R语言）

Posted 2020-10-26 szcf715

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了用一个简单的例子比较SVM，MARS以及BRUTO（R语言）相关的知识，希望对你有一定的参考价值。

背景重述

本文是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)))
}

模型训练

SVM直接调用e1071包中的svm函数
BRUTO和MARS都是调用mda包，且由于两者都是用于回归，所以转换为分类时，是比较拟合值与类别标签的距离，划分到越靠近的那一类
原书中提到实验中MARS不限定阶数，但实际编程时，设置阶数为10

交叉验证选择合适的\$C\$

我分两步进行选择：

粗选：在较大范围内寻找最优的\$C\$
细分：在上一步选取的最优值附近进行细分

注意避免最优值取在边界值。以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\$，比如某次实验结果如下：

Method	best cost
SV Classifier	2.6
SVM/poly 2	1
SVM/poly 5	28
SVM/poly 10	0.5

当然，实际中我们并不需要重新设置参数来训练模型，因为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，大于则划为第二类，否则第一类。