支持向量机原理及R语言实现

Posted 2021-04-18 Super越

距离上一篇文章竟过去快一个月了。本次主要是想讲一下支持向量机。

支持向量机是一种二分类模型，最基本的模型定义为特征空间中最大间隔的线性分类器，其学习的优化目标便是间隔最大化。什么意思呢？

假设现在的数据是线性可分的，用一条直线将其划分开，再假设样本分布如上图所示，能将样本点分开的直线有L1/L2/L3.....，有很多很多，但是什么样的直线使得分类的鲁棒性最好呢？也就是当有新的样本放入时，直线对这些点的适应性最强，错误分类的可能性最小。例如下图：红色的点为新样本，那么L1这条直线就将其错误分类了,说明L1这条直线的鲁棒性并不是最好的，因此支持向量机就是找到一条鲁棒性最好的直线。

那么如何寻找这条直线呢？

简单的说可以理解为两个支持向量相连的中垂线。

 

假设样本中有一类为正类（>1），另外一类为负类（<-1）

两条虚线对应着安全边际，黑色样本点称为支持向量。如果支持向量发生移动，就会导致边际和决策范围发生改变。分类器之间的距离被称为边际（两条虚线距离）。

分离边缘最大化等价于使权重向量范数||w||最小化。即求1/||w||的最大值相当于求||w||^2的最小值。

目标函数：min_w,b (1/2*||w||²)

但是存在约束条件：即蓝色圆点的样本点用w^Tx+b=2,3,4,5....(w^Tx+b>1),蓝色三角形的样本点用w^Tx+b=-2,-3,-4,-5.....(w^Tx+b<-1)表示；那么所有不是支持向量的点，有y*(w^Tx+b)>1；什么意思呢？y:表示真实值，圆点为正类为正数，正正得正，而三角形为负类为负数，负负也得正。

因此目标函数受到约束：

min_w,b (1/2*||w||²)

    s.t. y⁽ⁱ⁾(w^Tx⁽ⁱ⁾+b)≥1；i=1,2,3,...,n

约束条件转化一下：

g(w)=-y⁽ⁱ⁾(w^Tx⁽ⁱ⁾+b)+1≤0

题外话：最优化问题

1、min f(x)

属于无条件约束，对变量求导并令其等于0，求极值即可。

2、s.t.hi(x)=0 i=1,2,3....n

属于等式约束，需构建拉格朗日函数，将有约束问题转化为无约束条件，再按照无约束条件求解。

其中，λ为拉格朗日乘子

3、gi(x)≤0 i=1,2,3,4....n

属于等式约束和不等式约束。对于等式约束引入拉格朗日函数，不等式引入松弛变量。将这种约束条件转化为无约束条件，再进行求解。

其中，λ为拉格朗日乘子；μ为KKT乘子（注：KKT乘子是大于0的）

回归正题：

因此可以将目标函数拉格朗日化：

对于拉格朗日乘子发的对偶问题

min_w,b max_a L(w,b,a)  等价于 max_a min_w,b L(w,b,a)

第一步：求L(w,b,a)最小化。

因此，分别对w和b求偏导，并令其等于0；

将w带入目标函数拉格朗日化中：

第二步：对a 求极大值

求极大值不好求，可以返过来求极小值，即将等式添加一个负号，等式变为

由于篇幅有限，不再往下推到了。求导求出a值后即可计算出最优的w值

那么b值是什么呢？

通过带入公式wx+b≥1和wx+b≤1可以计算最优偏置:b=1-wx;那么线性可分问题得到的最优分类判别函数为：

当样本中存在一些异常值时，即样本不满足y(wx+b)≥1条件的约束，而出现分类误差。

像红色的点就不满足约束条件吧，其真实值是负数，预测值是正数。两者乘积不满足≥1的条件。

因此必须放宽约束条件:y(wx+b)≥1-§，其中§为松弛变量，并且§≥0。

松弛变量：度量一个样本点对线性可分约束条件的偏离程度。当§=0时，样本点就是一个支持向量了；0≤§≤1，时样本点是落在分离区域内的，且在分类超平面的正确一边。§≥1时，样本点进入分类超平面错误一边。

目标函数：

剩下的推到过程就和上面是一样的步骤，这里就不再展示了。

核函数

当样本数据线性不可分时，可以使用核函数进行非线性支持向量机，其是将输入向量映射到一个高维特征向量空间。例如：

样本分布如上图所示如何划分？

可以通过二维的上升到三维中。例如：高斯分布，将中间红色的点向上提，蓝色的在底部（如下图）这样子可以在中间用一个平面就可以划分开了。

【图片来源：网络】

支持向量机数学原理部分暂时讲到这里，如果想了解整个数学推导过程的可以去看一些机器学习的书籍，例如很多人推崇的西瓜书，李航的书籍讲得比较完整，也可以来问笔者，一起学习。若笔者上述原理中表达有错误的欢迎指出，下面是R语言实践部分。

R语言实现

数据集：MASS包中的糖尿病数据集

该数据包含532条样本即532位观察者的信息，分别存储在两个数据集中。其中数据表有8个特征变量，分别为：

npreg ：表示观察者怀孕次数

glu ：表示观察者血糖浓度，由口服葡萄糖耐量测试给出

bp ：表示观察者舒张压（单位为mm Hg）

skin ：表示观察者三头肌皮褶厚度（单位为mm）

bmi ：表示观察者身体质量指数

ped ：表示观察者糖尿病家族影响因素

age ：表示观察者年龄

type ：表示观察者是否患有糖尿病

从数据结构从可以看出变量之间的值域是不等的，因此需进行标准化后再进行后续分析。

pima.scale <- scale(pima[, 1:7]) %>% as_tibble() %>% bind_cols(type = pima$type)pima.scale %>% melt(id.var = "type") %>% ggplot(aes(type, value)) + geom_boxplot(outlier.colour = "red") + facet_wrap(~variable,ncol =2)

貌似除了血糖浓度(glu）外其他的变量看不出存在明显的差异，间接说明血糖浓度与观察者是否患有糖尿病存在很大的关系。从图中也可以看出全部的变量都存在疑似异常值的数据（红色点），以下的分析均将其归为正常值。

pima.scale[,-8] %>% cor() %>% corrplot::corrplot.mixed()

从图中可以看到，变量npreg和变量age、 变量skin和变量bmi 存在一定的相关性。可以研究一下是否存在共线性问题（我这里就不再研究了）

检查样本数据分布是否均衡

table(pima.scale$type) No Yes 355 177

从结果来看样本数据分布是均衡的。若存在不均衡的样本数据，可以通过一、欠采样法：该方法主要是对大类进行处理它会减少大类的观测数来使得数据集平衡，，也就是从大类中抽一部分数据进行后续分析，但是会导致不少的重要信息缺失，这种方法一般是数据集整体很大时使用。二、过采样法：针对小类进行处理它会以重复小类的观测的方式来平衡数据，类似bootstrop抽样，缺点就是容易造成过拟合。三、人工合成数据：其也属于过采样法的一种，利用模型或其他抽样方法生成新的数据而不是重复原始观测来解决不平衡性。例如使用K近邻法，将附近的N个点的数据取均值生成一个新的数据，也可以使用加权K近邻考虑相关性在里面。方法还有很多这里就不列举了。

按照三七分的方式划分训练集和测试集

set.seed(2020)index <- sample(1:2,nrow(pima.scale),replace = TRUE,prob = c(0.7,0.3))pima_train <- pima.scale[index==1,]pima_text <- pima.scale[index==2,]

建模

本文涉及到e1071包中的svm和tune.svm函数；其参数设置差不多。

svm(formula, data = NULL, ..., subset,na.action =na.omit, scale = TRUE)svm(x, y = NULL, scale = TRUE, type = NULL,kernel ="radial", degree = 3,gamma = if (is.vector(x)) 1 else 1 / ncol(x),coef0 = 0, cost = 1, nu = 0.5,class.weights = NULL, cachesize = 40,tolerance = 0.001, epsilon = 0.1,shrinking = TRUE,cross = 0, probability = FALSE, fitted = TRUE,..., subset, na.action = na.omit)tune.svm(x, y = NULL, data = NULL, degree = NULL, gamma = NULL, coef0 = NULL,cost = NULL, nu = NULL, class.weights = NULL, epsilon = NULL, ...)best.svm(x, tunecontrol = tune.control(), ...)

参数说明 ：

formula：表示R中公式，例如：y~.。

data：表示训练数据集。

subset：为一个索引向量，指定训练样本数据。

na.action： 逻辑值，是对缺 失值的处理方法，默认是删除缺失值所在的行。

scale：逻辑值，是否对变量进行标准化处理，默认是进行标准化处理。

x：表示特征变量，可以是矩阵、向量，也可以是稀疏矩阵。

y：表示为分类变量。

type：指定建模的类别，支持向量机通常用于分类、回归和异常值检测，默认情况下，svm模型根据因变量y是否为因子，type选择C-classification或eps-regression。用于分类的有：C-classification、nu-classification、one-classification。回归的有：eps-regression、nu-regression。

kernel：指定建模过程中使用的核函数，目的在于解决支持向量机线性不可分问题。有四类核函数可选，即线性核函数（linear）、多项式核函数（polynomial）、径向基（高斯）核函数（radial）和神经网络核函数（sigmoid）。对应的核函数需要指定的参数：例如多项式核函数，可以指定：gamma、coef0和degree。可能有人会有疑问，是不是认为线性核函数就没有可以指定的参数了？不知道还记得前面讲原理部分，当样本存在异常值时，引入松弛变量，有一个惩罚系数，在R中用cost参数表示。

degree：用于多项式核函数的参数，默认为3。

gamma：用于除线性核函数之外的所有核函数参数，默认为1。

coef0：用于多项式核函数和神经网络核函数的参数，默认为0。

nu：用于nu-classification、nu-regression和one-classification回归类型中的参数。

class.weights：指定类权重。

cachesize：默认缓存大小为40M。

cross：可为训练集数据指定k重交叉验证。

probability：逻辑值，指定模型是否生成各类的概率预测，默认不产生概率值。

fitted：逻辑值，是否将分类结果包含在模型中，默认生成拟合值。

线性核函数

#线性核函数svm.linear <- tune.svm(type~.,data = pima_train,kernel="linear",cost=c(0.01,0.05,0.1,0.5,1,5,10))summary(svm.linear)Parameter tuning of ‘svm’:- sampling method: 10-fold cross validation 
- best parameters: cost 0.01
- best performance: 0.2125178 
- Detailed performance results: cost error dispersion1 0.01 0.2125178 0.070637502 0.05 0.2312945 0.072947073 0.10 0.2313656 0.086956724 0.50 0.2207681 0.083157165 1.00 0.2234708 0.089077316 5.00 0.2261024 0.088172717 10.00 0.2261024 0.08817271

从模型使用7个（0.01,0.05,0.1,0.5,1,5,10）不同的惩罚系数训练模型，通过十阶交叉验证得到最优的模型是惩罚系数为0.01的模型，模型的分类误差率为21.2%。将模型在测试集上进行预测和检验：

svm.linear_pred <- predict(svm.linear$best.model,newdata = pima_text)flag <- table(svm.linear_pred,pima_text$type) %>% print(.) svm.linear_pred No Yes No 95 27 Yes 7 26#计算准确率accu.svm.linear <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)[1] 0.7806452

线性核函数的准确率为78%，再看看非线性模型表现怎么样，依然使用交叉验证方式对模型参数进行选择。注：多项式的阶数不要设置太高，容易过拟合。

svm.poly <- tune.svm(type~.,data = pima_train,kernel="polynomial",degree = c(2:5),coef0 = c(seq(0,3,by=0.5)))summary(svm.poly)Parameter tuning of ‘svm’:- sampling method: 10-fold cross validation 
- best parameters: degree coef0 2 0.5
- best performance: 0.2386913 
- Detailed performance results: degree coef0 error dispersion1 2 0.0 0.3128023 0.045069742       3   0.0 0.2492888 0.04119076**************展示部分**************27 4 3.0 0.2996444 0.0748187128 5 3.0 0.3339972 0.07827219

一共有4*7=28个模型，其中最优的模型多项式是2阶的，它的核系数为0.5，将这两个参数形成的模型在测试集上进行预测

svm.poly_pred <- predict(svm.poly$best.model,newdata = pima_text)flag <- table(svm.poly_pred,pima_text$type) %>% print(.) 
svm.poly_pred No Yes No 94 29 Yes 8 24#计算准确率accu.svm.poly <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)[1] 0.7612903

从结果上看多项式核函数的准确率才76%，模型的表现不如线性核函数的模型。

接下来就是径向基核函数，该核函数仅需设定gamma参数，注：若该参数设置过小，模型就不能解决决策边界的复杂性。若参数设置过大。模型容易造成过拟合。

#径向基核函数svm.radial <- tune.svm(type~.,data = pima_train,kernel="radial",gamma = c(seq(0.5,3,by=0.5)))summary(svm.radial)
Parameter tuning of ‘svm’:
- sampling method: 10-fold cross validation 
- best parameters: gamma 0.5
- best performance: 0.283357 
- Detailed performance results: gamma error dispersion1 0.5 0.2833570 0.080148332 1.0 0.3016358 0.109434303 1.5 0.3152205 0.099800374 2.0 0.3285917 0.110078295 2.5 0.3285917 0.104790076 3.0 0.3285917 0.10479007

最优的模型中参数 gamma 值是0.5，在测试集上测试

svm.radial_pred <- predict(svm.radial$best.model,newdata = pima_text)flag <- table(svm.radial_pred,pima_text$type) %>% print(.) 
svm.radial_pred No Yes No 93 30 Yes 9 23> #计算准确率> accu.svm.radial <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)[1] 0.7483871

径向基核函数的表现更差了，准确率仅有74%。最后只能看看S函数核函数的表现了。

svm.sigmoid <- tune.svm(type~.,data = pima_train,kernel="sigmoid",gamma = c(0.1,0.5,1,2,3),coef0 = c(0.1,0.5,1,2,3))summary(svm.sigmoid)
Parameter tuning of ‘svm’:
- sampling method: 10-fold cross validation 
- best parameters: gamma coef0 0.1 2
- best performance: 0.2440256 
- Detailed performance results: gamma coef0 error dispersion1 0.1 0.1 0.2466572 0.039342002    0.5   0.1 0.3020626 0.04304675**************展示部分**************24 2.0 3.0 0.3285917 0.0841480825 3.0 3.0 0.3127312 0.06223051

S函数核函数有5*5=25个模型，其中最优的模型参数为：gamma=0.1 coef0=2。

svm.sigmoid_pred <- predict(svm.sigmoid$best.model,newdata = pima_text)flag <- table(svm.sigmoid_pred,pima_text$type) %>% print(.) 
svm.sigmoid_pred No Yes No 93 24 Yes 9 29 #计算准确率accu.svm.sigmoid <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)[1] 0.7870968

将四种核函数的准确率合在一样。

tibble(linear=accu.svm.linear,polynomial=accu.svm.poly,radial=accu.svm.radial,sigmoid=accu.svm.sigmoid) %>% print()# A tibble: 1 x 4 linear polynomial radial sigmoid <dbl> <dbl> <dbl> <dbl>1 0.781 0.761 0.748 0.787

sigmoid核函数的模型准确率和线性核函数是差不多的，而径向基核函数的准确率是最低的。

拓展一 ：采用sigmoid核函数，进一步提高分类效果，可以考虑使用svm函数中的class.weights参数为不同的类指定权重，规则是：某个类与其他类存在明显差异时，考虑给予较小的权重，某些类之间差异较小时，考虑给予较大的权重。

weights <- c(2.5,2)names(weights) <- c("Yes","No")model <- svm(type~.,data = pima_train,kernel="sigmoid",gamma=0.1,coef0=2,class.weights = weights)pred <- predict(model, newdata = pima_text)flag <- table(pred,pima_text$type) %>% print(.)
pred No Yes No 93 22 Yes 9 31#计算准确率accu <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)[1] 0.8

虽然提高得不是很明显，但是还是有作用的。 由于篇幅有限就不 再对其他的参数进行尝试了，私下可自行去尝试其他的参数。

拓展二 ：多类分类

运用拆分思想：“一对一”、“一对其余”

一对一：假定有N个真实类别，将N个类别进行两两配对（一个正一个反），从而产生（N-1）/2个二分类学习器，最终通过投票方式产生结果。

一对其余：假定有N个真实类别，每一次将其中一个类做正类，其余类别为负类，从而产生N个分类器，则得到N个结果，再通过决策方式确定最终结果。

当然还有其他的多类分类方法，有一些可能理解起来比较困难，暂时就不讲先。

R完整代码

#清除缓存rm(list=ls())gc()#加载数据包library(MASS)library(dplyr)library(reshape2)library(ggplot2)library(e1071)
#描述性分析pima <- rbind(Pima.tr, Pima.te)str(pima)table(pima$type)pima.scale <- scale(pima[, 1:7]) %>% as_tibble() %>% bind_cols(type = pima$type)pima.scale %>% melt(id.var = "type") %>% ggplot(aes(type, value)) + geom_boxplot(outlier.colour = "red") + facet_wrap(~variable,ncol =2)pima.scale[,-8] %>% cor() %>% corrplot::corrplot.mixed()table(pima.scale$type)
#划分训练集和测试集set.seed(2020)index <- sample(1:2,nrow(pima.scale),replace = TRUE,prob = c(0.7,0.3))pima_train <- pima.scale[index==1,]pima_text <- pima.scale[index==2,]dim(pima_train)dim(pima_text)
#线性核函数svm.linear <- tune.svm(type~.,data = pima_train,kernel="linear",cost=c(0.01,0.05,0.1,0.5,1,5,10))summary(svm.linear)svm.linear_pred <- predict(svm.linear$best.model,newdata = pima_text)flag <- table(svm.linear_pred,pima_text$type) %>% print(.) #计算准确率accu.svm.linear <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)
#多项式核函数svm.poly <- tune.svm(type~.,data = pima_train,kernel="polynomial",degree = c(2:5),coef0 = c(seq(0,3,by=0.5)))summary(svm.poly)svm.poly_pred <- predict(svm.poly$best.model,newdata = pima_text)flag <- table(svm.poly_pred,pima_text$type) %>% print(.) #计算准确率accu.svm.poly <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)
#径向基核函数svm.radial <- tune.svm(type~.,data = pima_train,kernel="radial",gamma = c(seq(0.5,3,by=0.5)))summary(svm.radial)svm.radial_pred <- predict(svm.radial$best.model,newdata = pima_text)flag <- table(svm.radial_pred,pima_text$type) %>% print(.) #计算准确率accu.svm.radial <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)
#S函数核函数svm.sigmoid <- tune.svm(type~.,data = pima_train,kernel="sigmoid",gamma = c(0.1,0.5,1,2,3),coef0 = c(0.1,0.5,1,2,3))summary(svm.sigmoid)svm.sigmoid_pred <- predict(svm.sigmoid$best.model,newdata = pima_text)flag <- table(svm.sigmoid_pred,pima_text$type) %>% print(.)#计算准确率accu.svm.sigmoid <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)
#四种核函数准确率整合tibble(linear=accu.svm.linear,polynomial=accu.svm.poly,radial=accu.svm.radial,sigmoid=accu.svm.sigmoid) %>% print()
#加入权重weights <- c(2.5,2)names(weights) <- c("Yes","No")model <- svm(type~.,data = pima_train,kernel="sigmoid",gamma=0.1,coef0=2,class.weights = weights)pred <- predict(model, newdata = pima_text)flag <- table(pred,pima_text$type) %>% print(.)#计算准确率accu <- (sum(diag(flag))/nrow(pima_text)) %>% print(.)

感谢您的支持