【机器学习实战】-- Titanic 数据集（4）-- 支持向量机

1. 写在前面: 

本篇属于实战部分，更注重于算法在实际项目中的应用。如需对感知机算法本身有进一步的了解，可参考以下链接，在本人学习的过程中，起到了很大的帮助：

统计学习方法 李航

支持向量机原理 https://www.cnblogs.com/pinard/p/6097604.html 等共5篇

空间中任意一点到超平面距离的公式推导 https://www.cnblogs.com/yanganling/p/8007050.html

拉格朗日对偶 https://www.cnblogs.com/ooon/p/5723725.html

 

2. 数据集：

数据集地址：https://www.kaggle.com/c/titanic

Titanic数据集是Kaggle上参与人数最多的项目之一。数据本身简单小巧，适合初学者上手，深入了解比较各个机器学习算法。

数据集包含11个变量：PassengerID、Pclass、Name、Sex、Age、SibSp、Parch、Ticket、Fare、Cabin、Embarked，通过这些数据来预测乘客在Titanic事故中是否幸存下来。

 

3. 算法简介：

支持向量机原本是一种分类模型，和感知机一样，定义了特征空间中的一个超平面，但间隔最大化使其区别于感知机（通常是优于）。此外，和感知机一样，利用Gram矩阵，方便地使用核函数，使其可以成为一个非线性分类器。

支持向量机的学习策略就是间隔最大化，并且等价与合页损失函数（hinge loss）的最小化问题。

 

3.1 线性可分支持向量机模型：

给定一个数据集：$T=\left \{ \left ( x_{1}, y_{1} \right ), \left ( x_{2}, y_{2} \right ), ..., \left ( x_{N}, y_{N} \right ) \right \}$，其中$x_{i}\in \chi\subseteq \bf{R^{n}}$，$y_{i} \in \Upsilon= \left \{+1, -1 \right \}$，$i = 1,2,...,N$。这代表数据集共有 N 对 实例，每个实例 $x_{i}$都是n维的。

假设通过间隔最大化求解得到的分离超平面为：

$$ w^{\ast } \cdot x + b^{\ast} = 0 $$

则如下分类函数即线性可分支持向量机模型：

$$ f(x) = \text{sign}(w^{\ast} \cdot x + b^{\ast}) $$

 

 

3.2 线性可分支持向量机损失函数：

首先，在定义损失函数之前，需要了解函数间隔以及几何间隔：

超平面$(w,b)$关于样本点$(x_{i}, y_{i})$的函数间隔（使用于感知机损失函数）：

$$ \hat{\gamma_{i}} = y_{i} (w \cdot x_{i} + b) $$

超平面$(w,b)$关于训练数据集$T$的函数间隔：

$$ \hat{\gamma} = \min_{i=1,...,N} \hat{\gamma_{i}} $$

 

超平面$(w,b)$关于样本点$(x_{i}, y_{i})$的几何间隔：

$$ \gamma_{i} = \frac{y_{i} (w \cdot x_{i})}{\left \| w \right \|} $$

超平面$(w,b)$关于训练数据集$T$的几何间隔：

$$ \gamma = \min_{i=1,...,N}\gamma_{i} $$

 

回想一下，支持向量机的学习策略就是间隔最大化，即最大化超平面$(w,b)$关于训练数据集的几何间隔$\gamma$，可以将其表示为如下的最优化问题：

$$ \max_{w,b} \gamma \quad s.t. \frac{y_{i} (w \cdot x_{i})}{\left \| w \right \|} \geq \gamma , \quad i=1,2,...,N $$

通过函数间隔与几何间隔的关系，将其改写成如下形式：

$$ \max_{w,b} \frac{\hat{\gamma}}{\left \| w \right \|} \quad s.t. y_{i}(w \cdot x_{i} + b)  \geq \hat{\gamma}, \quad i=1,2,...,N$$

显而易见，函数间隔$\hat{\gamma}$的取值不会影响上述优化问题的解，我们可以令其为1。同时将最大化$\frac{1}{||w||}$问题转化为最小化$||w||^{2}$问题，同时为了便于后续的计算，加上$\frac{1}{2}$的系数，于是优化问题转变为：

$$ \min_{w,b} \frac{1}{2}||w||^{2} \quad s.t. y_{i}(w \cdot x_{i} + b) \geq 1, \quad i=1,2,...,N$$

为了求解上述不等式优化问题，可以应用拉格朗日乘子法，得到拉格朗日函数$L$，L即线性可分支持向量机的损失函数：

$$ L(w, b, \alpha) = \frac{1}{2}\left \| w \right \|^2 + \sum_{i=1}^{N}\alpha_{i}\left [ 1 - y_{i}(w \cdot x + b) \right ] $$

通过拉格朗日函数可以得到：

$$ \min_{w,b}\frac{1}{2}\left \| w \right \|^2 = \min_{w,b} \max_{\alpha_{i} \geq 0} L(w,b,\alpha) $$

然后，将其转化成对偶问题：

$$ \min_{w,b}\frac{1}{2}\left \| w \right \|^2 = \max_{\alpha_{i} \geq 0}\min_{w,b} L(w,b,\alpha)  $$

 

3.3 线性可分支持向量机的学习过程：

线性可分支持向量机的学习过程即求解对偶问题的过程：首先求$L(w,b,\alpha)$对$w,b$的极小，再求其对$\alpha$的极大。

第一步：对求$L(w,b,\alpha)$对$w,b$的极小：

$$ \nabla_{w}L(w,b,\alpha) =  w - \sum_{i=1}^{N} \alpha_{i}y_{i}x_{i} = 0 $$

$$ \nabla_{b}L(w,b,\alpha) = \sum_{i=1}^{N} \alpha_{i}y_{i} = 0 $$

可以得到如下等式：

$$ w = \sum_{i=1}^{N} \alpha_{i}y_{i}x_{i} $$

$$ \sum_{i=1}^{N} \alpha_{i}y_{i} = 0 $$

其中，第一个式子将带回原式$L(w,b,\alpha)$，第二个式子将作为第二步求解对$\alpha$的极大时的约束条件。

将第一个式子带回原式$L(w,b,\alpha)$可得：

$$

\begin{align*}
L(w,b,\alpha) &= \frac{1}{2}\left \| w \right \|^2 + \sum_{i=1}^{N}\alpha_{i}(1 - y_{i}w \cdot x_{i}) - y_{i}b \\
&= \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i} - \sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + 0 \\
&= - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i}
\end{align*}

$$

第二步：求$\min_{w,b} L(w,b,\alpha)$对$\alpha$的极大，配合约束条件$\sum_{i=1}^{N} \alpha_{i}y_{i} = 0$ 可得：

$$ \max_{\alpha} - \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) + \sum_{i=1}^{N} \alpha_{i} \quad s.t. \quad \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \quad \alpha_{i} \geq0, \quad i=1,2,...,N  $$

$$ \Downarrow $$

$$ \min_{\alpha}  \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\left ( x_{i} \cdot x_{j} \right ) - \sum_{i=1}^{N} \alpha_{i} \quad s.t. \quad \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \quad \alpha_{i} \geq0, \quad i=1,2,...,N  $$

以上问题可以通过SMO算法(这里不详细解释)求得$\alpha$的解$\alpha^{\ast}$

最终求得$w^{\ast}$ 和 $ b^{\ast} $ 为：

$$ w^{\ast} = \sum_{i=1}^{N} \alpha_{i}^{\ast}y_{i}x_{i} $$

$$ b_{s}^{\ast} = y_{s} - \sum_{i=1}^{N} \alpha_{i}^{\ast}y_{i}x_{i} \cdot x_{s}， \quad \text{ for $\alpha_{s} > 0$， $(x_{s}, y_{s})$是支持向量 } $$

 

3.4 线性支持向量机：

3.3中介绍的线性可分支持向量机仅适用于线性可分训练数据集，当训练数据集线性不可分时，约束条件 $y_{i}(w \cdot x_{i} + b) \geq \hat{\gamma}, \quad i=1,2,...,N$ 并不对所有数据点都成立。。我们将原来的约束条件称之为硬间隔最大化，为了将支持向量机扩展到线性不可分数据集，我们对其进行一定的修改，修改后的约束条件称之为软间隔最大化。

既然有些样本点$(x_{i}, y_{i})$不能满足函数间隔大于等于1的约束条件，我们对每个样本点$(x_{i}, y_{i})$引入一个松弛变量$\xi_{i} \geq 0$，使函数间隔加上松弛变量大于等于1。同时，由于加入了松弛变量，需要付出代价，在目标函数中加对每个松弛变量加入代价。

$$ \begin{align*}
& \min_{w,b,\xi} \frac{1}{2}\left \| w \right \|^{2} + C\sum_{i=1}^{N}\xi_{i} \\
\text{s.t.} \quad & y_{i}(w \cdot x_{i} + b) \geq 1 - \xi_{i}, \\
& \xi_{i} \geq 0, \quad i=1,2,...,N
\end{align*} $$

依照3.3中相似的步骤，可以得到拉格朗日函数$L$：

$$ L(w,b,\xi,\alpha,\mu) = \frac{1}{2}\left \| w \right \|^2 + C\sum_{i=1}^{N}\xi_{i} + \sum_{i=1}^{N}\alpha_{i}(1-\xi_{i}-y_{i}w \cdot x_{i} - y_{i}b) - \sum_{i=1}^{N} \mu_{i}\xi_{i} \quad s.t. \quad \alpha_{i} \geq 0, \quad \mu_{i} \geq 0  $$

接着求解拉格朗日对偶问题 $ \max_{\alpha, \mu} \min_{w,b,\xi} L(w,b,\xi, \alpha, \mu) $ 可得：

$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}(x_{i}\cdot x_{j}) - \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$

显而易见，其形式和线性可分支持向量机一致，只是约束条件变了， 同样可以利用SMO算法求解$\alpha^{\ast}$，进而解得$w^{\ast}$和$b_{s}^{\ast}$。

3.5 线性支持向量机与合页损失函数：

如3.4节所述：线性支持向量机，其模型为分离超平面$w^{\ast} \cdot x + b^{\ast} = 0$ 以及决策函数 $f(x) = \text{sign}(w^{\ast} \cdot x + b^{\ast})$，其学习策略为软间隔最大化。但其还有另一种解释，就是带有正则化项的合页损失函数最小化：

$$ \min_{w,b} \quad \sum_{i=1}^{N}\left [ 1-y_{i}(w \cdot x_{i} +b) \right ]_{+} +  \lambda \left \| w \right \|^2$$

上述式子中第一项$ [1-y(w \cdot x + b)]_{+} $为合页损失函数(hinge loss)，第二项为参数$w$的L2正则化项。

令 $1-y_{i}(w \cdot x_{i} + b) = \xi_{i}， \quad \xi_{i} \geq 0$，且 $ \lambda = \frac{1}{2C}$，则：

$$ \min \quad \sum_{i=1}^{N} [1-y_{i}(w \cdot x_{i} +b)]_{+} + \lambda||w||^2 $$

$$ \Downarrow $$

$$ \min \quad \sum_{i=1}^{N} \xi_{i} + \lambda||w||^2 $$

$$ \Downarrow $$

$$ \min \quad \frac{1}{C}(\frac{1}{2}||w||^2 + C\sum_{i=1}^{N}\xi_{i} ) $$

可以看到其最小化的函数和3.4节中的一致。

3.6 非线性支持向量机与核函数：

非线性分类问题，是指通过利用非线性模型才能很好地进行分类的问题。求解非线性问题的一个方法是进行非线性变换，将非线性问题变换为线性问题，通过求解变换后的线性问题的方法求解原来的非线性问题。（类似于将二元的多项式回归转换成五元的线性回归）。总结一下，使用线性分类方法求解非线性分类问题分为两步：

第一步：使用一个变换将原空间的数据映射到新空间

第二步：在新空间里用线性分类学习方法从训练数据中学习分类模型

按照上述理论，我们定义一个从输入空间 $\mathcal{X}$ 到特征空间 $\mathcal{H}$ 的映射：$\phi(x): \mathcal{X} \rightarrow \mathcal{H}$，将这一映射应用到线性支持向量机的优化函数中，得到：

$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}\phi(x_{i})\cdot \phi(x_{j}) - \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$

但这样做还存在一个问题，那就是实例的内积 $\phi(x_{i}) \cdot \phi(x_{j})$ 是在映射后的空间$\mathcal{H}$中进行的，这是一个高维空间，计算内积并不容易且费时，于是，我们使用了核技巧，定义了核函数：

$$ K(x,z) = \phi(x) \cdot \phi(z) $$

核技巧的想法是，在学习与预测中只定义核函数$K(x,z)$，而不显示地定义映射函数$\phi$。通常，计算$K(x,z)$比较容易，因为这是在低维空间上进行的。于是，应用了核函数之后的支持向量机的优化函数为：

$$ \begin{align*}
&\min_{\alpha} \quad \frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_{i}\alpha_{j}y_{i}y_{j}K(x_{i},x_{j})- \sum_{i=1}^{N}\alpha_{i} \\
s.t. \quad & \sum_{i=1}^{N}\alpha_{i}y_{i}=0, \\
& 0 \leq \alpha_{i} \leq C, \quad i=1,2,...,N
\end{align*} $$

定义一个核函数，需要各种满足复杂的数学条件，一般我们使用的常见核函数包括：线性核函数(线性支持向量机中所使用的)、高斯核函数(也称作径向基核函数)，其他还包括多项式核函数、Sigmoid核函数：

线性核函数：$ K(x,z) = x \cdot z $

高斯核函数：$ K(x,z) = \exp( -\gamma \| x-z \|^2 ) $

多项式核函数：$ K(x,z) = (\gamma x \cdot z + r)^d $

Sigmoid核函数：$ K(x,z) = \tanh( \gamma x \cdot z + r ) $

以上核函数中的参数 $\gamma$, $r$, $d$ 都需要根据实际情况，调参选取最优值。

 

4. 实战：

 sklearn中，其实有两个模块可以应用支持向量机，分别是 sklean.SVM，其中包括了线性 LinearSVC 以及 带核函数的 SVC；另一个是 sklearn.linear_model.SGDClassifier，只需要将loss设置为 ‘hinge’ 就等同于线性支持向量机。

以下代码中，我们分别使用了 LinearSVC()，SGDClassifier，以及 SVC(kernel='rbf') 进行模型训练，供参考：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler, StandardScaler, OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.model_selection import StratifiedKFold, GridSearchCV
from sklearn.pipeline import Pipeline, FeatureUnion
from sklearn.linear_model import SGDClassifier
from sklearn.svm import LinearSVC, SVC
from sklearn.metrics import accuracy_score, precision_score, recall_score
from sklearn.base import BaseEstimator, TransformerMixin


class DataFrameSelector(BaseEstimator, TransformerMixin):
    def __init__(self, attribute_name):
        self.attribute_name = attribute_name

    def fit(self, x, y=None):
        return self

    def transform(self, x):
        return x[self.attribute_name].values


# Load data
data_train = pd.read_csv('train.csv')

train_x = data_train.drop('Survived', axis=1)
train_y = data_train['Survived']

# Data cleaning
cat_attribs = ['Pclass', 'Sex', 'Embarked']
dis_attribs = ['SibSp', 'Parch']
con_attribs = ['Age', 'Fare']

# encoder: OneHotEncoder()、OrdinalEncoder()
cat_pipeline = Pipeline([
    ('selector', DataFrameSelector(cat_attribs)),
    ('imputer', SimpleImputer(strategy='most_frequent')),
    ('encoder', OneHotEncoder()),
])

dis_pipeline = Pipeline([
    ('selector', DataFrameSelector(dis_attribs)),
    ('scaler', StandardScaler()),
    ('imputer', SimpleImputer(strategy='most_frequent')),
])

con_pipeline = Pipeline([
    ('selector', DataFrameSelector(con_attribs)),
    ('scaler', StandardScaler()),
    ('imputer', SimpleImputer(strategy='mean')),
])

full_pipeline = FeatureUnion(
    transformer_list=[
        ('con_pipeline', con_pipeline),
        ('dis_pipeline', dis_pipeline),
        ('cat_pipeline', cat_pipeline),
    ]
)

train_x_cleaned = full_pipeline.fit_transform(train_x)

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=2)

# test1：using Linear SVC
clf1 = LinearSVC(tol=1e-4, max_iter=10000)
param_grid = [{
    'penalty': ['l1'],
    'loss':['squared_hinge'],
    'dual':[False],
    'C':[1e-4, 1e-3, 1e-2, 1e-1, 1, 10],
    'class_weight':['balanced', None]
                },
                {
    'penalty': ['l2'],
    'loss':['hinge', 'squared_hinge'],
    'dual':[True],
    'C':[1e-4, 1e-3, 1e-2, 1e-1, 1, 10],
    'class_weight':['balanced', None]
                }]

grid_search = GridSearchCV(clf1, param_grid=param_grid, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True)

grid_search.fit(train_x_cleaned, train_y)
predicted_y = grid_search.predict(train_x_cleaned)

df_cv_results = pd.DataFrame(grid_search.cv_results_)
print('-------Result of LinearSVC-------')
print(grid_search.best_params_)
print(accuracy_score(train_y, predicted_y))
print(precision_score(train_y, predicted_y))
print(recall_score(train_y, predicted_y))

# test2: Using SGDClassifier with hinge loss / squared hinge loss
clf2 = SGDClassifier(tol=1e-4, max_iter=10000)
param_grid2 = [{
    'loss': ['hinge', 'squared_hinge'],
    'penalty': ['l2', 'l1'],
    'alpha': [1e-4, 1e-3, 1e-2, 1e-1, 1, 10, 100, 1000],
    'class_weight':[None, 'balanced']
                }]

grid_search2 = GridSearchCV(clf2, param_grid=param_grid2, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True)

grid_search2.fit(train_x_cleaned, train_y)
predicted_y2 = grid_search2.predict(train_x_cleaned)

df_cv_results2 = pd.DataFrame(grid_search2.cv_results_)
print('-------Result of SGD-------')
print(grid_search2.best_params_)
print(accuracy_score(train_y, predicted_y2))
print(precision_score(train_y, predicted_y2))
print(recall_score(train_y, predicted_y2))

# test3: Using SVC
clf3 = SVC(kernel='rbf', tol=1e-4)
param_grid3 = [{
    'C': [0.1, 0.5, 0.75, 1, 2.5, 5, 10],
    'gamma': ['scale', 'auto', 2e-4, 2e-3, 2e-2]
                }]

grid_search3 = GridSearchCV(clf3, param_grid=param_grid3, cv=cv, scoring='accuracy', n_jobs=-1, return_train_score=True)

grid_search3.fit(train_x_cleaned, train_y)
predicted_y3 = grid_search3.predict(train_x_cleaned)

df_cv_results3 = pd.DataFrame(grid_search3.cv_results_)
print('-------Result of SVC with rbf-------')
print(grid_search3.best_params_)
print(accuracy_score(train_y, predicted_y3))
print(precision_score(train_y, predicted_y3))
print(recall_score(train_y, predicted_y3))

# 导入预测数据，预测结果，并生成csv文件
data_test = pd.read_csv('test.csv')
submission = pd.DataFrame(columns=['PassengerId', 'Survived'])
submission['PassengerId'] = data_test['PassengerId']

test_x_cleaned = full_pipeline.fit_transform(data_test)

submission_LinearSVC = pd.DataFrame(submission, copy=True)
submission_LinearSVC['Survived'] = pd.Series(grid_search.predict(test_x_cleaned))

submission_SVCrbf = pd.DataFrame(submission, copy=True)
submission_SVCrbf['Survived'] = pd.Series(grid_search3.predict(test_x_cleaned))

submission_LinearSVC.to_csv('submission_LinearSVC.csv', index=False)
submission_SVCrbf.to_csv('submission_SVCrbf.csv', index=False)

 

4.1 结果分析：

和前几篇一样，分别将 LinearSVC() 和 SVC(kernel='rbf')的最优参数用于预测集，并将结果上传kaggle，结果如下：

可以看到线性SVM和其他线性方法的结果差不多；而使用了核函数后，结果确实有了一定的提升。

	训练集 accuracy	训练集 precision	训练集 recall	预测集 accuracy（需上传kaggle获取结果）
朴素贝叶斯最优解	0.790	0.731	0.716	0.756
感知机	0.771	0.694	0.722	0.722
逻辑回归	0.807	0.781	0.690	0.768
线性SVM	0.801	0.772	0.684	0.773
rbf核SVM	0.834	0.817	0.731	0.785