[机器学习入门笔记] 3. 监督学习单模型部分

---

前言

进一步巩固机器学习基础，除了巩固现有理论，增加了代码实现，基于Numpy与sklearn，完全适合小白简单实战，本书主要参考了鲁伟.《机器学习》和李航.《统计学习方法》。

1、机器学习预备知识

1.1 关键术语与任务类型

机器学习的定义

系统通过计算手段利用经验来改善自身性能的过程（即通过分析和计算数据来归纳出数据普遍规律的学科）

关键术语

数据集：数据的集合

样本（实例）：数据集中的一条记录

特征维度：特征数量

任务类型

根据机器学习任务，将整个数据集划分为训练集和测试集

训练集：训练机器学习模型

测试集：验证模型在未知数据上的效果

1.2 机器学习三要素

机器学习方法由模型、策略、算法三要素组成，理解为机器学习模型在一定的优化策略下使用相应求解算法来达到最优目标的过程。

模型

要学习的决策函数或者条件概率分布

策略

在假设空间的各种模型中，按照一定标准选择最优模型。对于给定模型，模型输出F(x)和真实输出Y之间的误差用损失函数来度量。回归任务一般用均方误差，分类任务一般使用对数损失函数或者交叉熵损失函数

算法

学习模型的具体优化方法。当机器学习的模型和损失函数确定时，机器学习就具体形式化为最优化问题。常见的优化算法有梯度下降法、牛顿法、拟牛顿法等。

1.3 机器学习的核心

监督机器学习

所有监督学习都可以以通用的损失函数计算公式概括，其中第一项为针对训练集的经验误差项（训练误差），第二项为正则化项（惩罚项），用于堆模型复杂度的约束和惩罚，因此监督机器学习的核心任务是正则化参数的同时最小化经验误差。

过拟合

机器学习模型训练中，模型对训练数据学习过度，将数据中包含噪声和误差也学习了，使模型在训练集上表现的很好，而在测试集上表现的很差的一种现象。（特征工程、扩大训练集数量、算法设计和超参数调优可以防止过拟合）

1.4 机器学习流程

需求分析：

明确机器学习目标、输入是什么、目标输出是什么、是回归任务还是分类任务、关键性能指标，是结构化的机器学习任务还是基于深度学习的图像和文本任务，市面上项目相关产品有什么，对应的SOTA模型有哪些，相关领域的前沿研究和进展如何，项目的利弊。

数据采集：

企业数据源、公开数据或竞赛数据集、采用数据采集技术获取

数据清洗：

在生产环境下，数据比较脏，需要花大量时间清清洗数据

数据分析与可视化：

探索性数据分析和数据可视化，数据特征、目标变量分布、各自变量与目标变量是否需要可视化展示、数据中各变量缺失值如何，如何处理缺失值？

建模调优与特征工程：

训练完一个极限模型后，需要花大量时间进行模型调参和优化，结合业务的精细化特征工程比模型调参更能改善模型表现。

模型结果展示与分析报告：

经过一定的特征工程和模型调优后，一般有一个阶段性的最优模型结果，模型相应的关键性能指标达到最优状态。这是需要呈现模型，并对模型的业务含义解释。

模型部署与上线反馈优化

第 2 章 线性回归

线性回归是线性模型的一种典型方法，比如产品的销量预测、薪资水平预测。某种程度上讲回归分析不再是局限于线性回归这一具体模型和算法，更包含了广泛的由资源变量到因变量的机器学习建模思想。

2.1 线性回归的原理推导

线性回归就是通过训练学习得到一个线性模型最大程度的根据输入x拟合输出y，线性回归的关键问题在于确定参数w和b，使得你和拟合输出f(x)与真实输出y尽可能接近。

在回归任务中通常使用均方误差来度量预测与标签之间的损失，所以回归任务的优化目标就是拟合输出和真实输出之间的均方误差最小化。

为求得w和b的最小化参数w^*和b^*，基于上式分别对w和b求一阶导数并令其为0，推导过程如下：

2.2 线性回归的代码实现

线性回归模型代码编写思路

Numpy：

• 模型主体
  • 回归模型公式
  • 均方损失函数
  • 参数求偏导
• 训练过程
  • 参数初始化
  • 多轮迭代训练过程
  • 梯度下降的参数优化更新
• 模型测试
  • 测试结果
  • 可视化展示

sklearn：

• linear_model.LinearRegression
• 使用范例

2.2.1 基于Numpy的手动实现

定义线性回归模型主体

# 导入numpy模块
import numpy as np
### 定义模型主体部分
### 包括线性回归模型公式、均方损失函数和参数求偏导三部分
def linear_loss(X, y, w, b):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    w：变量参数权重矩阵
    b：偏置
    输出：
    y_hat：线性回归模型预测值
    loss：均方误差
    dw：权重系数一阶偏导
    db：偏置一阶偏导
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = x.shape[1]
    # 线性回归预测值
    y_hat = np.dot(X, w) + b
    # 计算预测值与实际标签之间的均方误差
    loss = np.sum(y_hat-y)**2 / num_train
    # 基于均方误差对权重系数的一阶梯度
    dw = np.dot(X.T, (y_hat-y)) / num_train
    # 基于均方误差对偏置的一阶梯度
    db = np.sum(y_hat-y) / num_train
    return y_hat, loss, dw, db

初始化模型参数

### 初始化模型参数
def initialize_params(dims):
    '''
    输入：
    dims：训练数据的变量维度
    输出：
    w：初始化权重系数
    b：初始化偏置参数
    '''
    # 初始化权重系数为零向量
    w = np.zeros((dims, 1))
    # 初始化偏置参数为零
    b = 0
    return w, b

定义线性回归的模型的训练过程

### 定义线性回归模型训练过程
def linear_train(X, y, learing_rate=0.01, epochs=10000):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    learning_rate：学习率
    epochs：训练迭代次数
    输出：
    loss_his：每次迭代的均方误差
    params：优化后的参数字典
    grads：优化后的参数梯度字典
    '''
    # 记录训练损失的空列表
    loss_his = []
    # 初始化模型参数
    w, b = initialize_params(X.shape[1])
    # 迭代训练
    for i in range(1, epochs):
        # 计算当前迭代的预测值、均方误差和梯度
        y_hat, loss, dw, db = linear_loss(X, y, w, b)
        # 基于梯度下降法的参数更新
        w += -learing_rate * dw
        b += -learing_rate * db
        # 记录当前迭代的损失
        loss_his.append(loss)
        # 每 10000 次迭代打印当前损失信息
        if i % 10000 == 0:
            print('epoch %d loss %f' %(i, loss))
        # 将当前迭代步优化后的参数保存到字典中
        params = {
            'w':w,
            'b':b
        }
        # 将当前迭代步的梯度保存到字典中
        grads = {
            'dw':dw,
            'db':db
        }
    return loss_his, params, grads

基于上述实现，使用sklearn的diabetes数据集测试，从sklearn中导入数据集，将其划分为测试集和训练集。

导入数据集

# 导入load_diabetes模块
from sklearn.datasets import load_diabetes
# 导入打乱数据函数
from sklearn.utils import shuffle
# 获取diabetes数据集
diabetes = load_diabetes()
# 获取输入
data = diabetes.data
# 获取标签
target = diabetes.target 

# 打乱数据
X, y = shuffle(data, target, random_state=13)
X = X.astype(np.float32)

# 训练集与测试集的简单划分
offset = int(X.shape[0] * 0.9)
# 训练集
X_train, y_train = X[:offset], y[:offset]
# 测试集
X_test, y_test = X[offset:], y[offset:]
# 将训练集改为列向量的形式
y_train = y_train.reshape((-1,1))
# 将测试集改为列向量的形式
y_test = y_test.reshape((-1,1))
# 打印训练集和测试集的维度
print('X_train=', X_train.shape)
print('X_test=', X_test.shape)
print('y_train=', y_train.shape)
print('y_test=', y_test.shape)

模型训练过程

# 线性回归模型训练
loss_list, loss, params, grads = linar_train(X_train, y_train, 0.001, 100000)
# 打印训练后得到的模型参数
print(params)

回归模型的预测函数

### 定义线性回归模型的而预测函数
def predict(X, params):
    '''
    输入：
    X：测试集
    params：模型训练参数
    输出：
    y_pred：模型预测结果
    '''
    # 获取模型参数
    w = params['w']
    b = params['b']

    # 预测
    y_pred = np.dot(X, w) + b    
    return y_pred
# 基于测试集的预测
y_pred = predict(X_test, params)
y_pred[:5]

回归模型R²系数

除均方误差外，回归模型的一个重要评估指标是R²系数，用来判断模型拟合水平

### 定义R^2系数函数
def r2_score(y_test, y_pred):
    '''
    输入：
    y_test：测试集标签值
    y_pred：测试集预测值
    输出：
    r2：R^2系数
    '''
    # 测试集标签均值
    y_avg = np.mean(y_test)
    # 总离差平方和
    ss_tot = np.sum((y_test-y_avg)**2)
    # 残差平方和
    ss_res = np.sum((y_test-y_pred)**2)
    ## R^2计算
    r2 = 1 - (ss_res/ss_tot)
    return r2
# 计算测试集的R^2系数
print(r2_score(y_test,y_pred))

利用 matplotlib 对预测结果和真值进行展示

import matplotlib.pyplot as plt
f = X_test.dot(params['w']) + params['b']

plt.scatter(range(X_test.shape[0]), y_test)
plt.plot(f, color = 'darkorange')
plt.xlabel('X')
plt.ylabel('y')
plt.show()

训练过程中损失的下降

plt.plot(loss_list, color = 'blue')
plt.xlabel('epochs')
plt.ylabel('loss')
plt.show()

2.2.2 基于sklearn的模型实现

基于sklearn的LinearRegression类给出对于该数据集的拟合效果，LinearRegression函数位于sklearn的linear_model模块下，定义该类的一个线性回归实例后，直接调用其fit方法拟合训练集

基于sklearn的线性回归模型

# 导入线性回归模块
from sklearn import linear_model
from sklearn.metrics import mean_suqared_error, r2_score
# 定义模型实例
regr = linear_model.LinearRegression()
# 模型拟合训练数据
regr.fit(X.train, y_train)
# 模型预测值
y_pred = regr.predict(X_test)
# 输出模型均方误差
print("Mean squared error:%.2f"% mean_squared_error(y_test, y_pred))
# 计算R^2系数
print('R Square score: %.2f' % r2_score(y_test,y_pred))

2.3 小结

基于均方误差最小化的最小二乘法是线性回归模型求解的基本方法，通过最小化均方误差和R²系数可以评估线性回归的拟合效果

第 3 章 对数几率回归

对数几率回归是一种线性分类模型，应用于业务场景，如信用卡场景下基于客户数据对其进行违约预测，互联网广告场景下预测用户是否点击广告，社交场景下判断垃圾邮件，基于体检数据判断是否患有某种疾病。

3.1 对数几率回归的原理推导

只需要知道一个单调可微函数将分类任务的真实标签y于线性回归模型的预测值进行映射，在对数几率回归中，需要找到一个映射函数将线性回归模型的预测值转换为0/1值。

Sigmoid函数单调可微，取值范围(0,1)。

将线性回归模型带入Sigmoid函数得，


为确定上式模型参数w和b，需要推导出损失函数并使其最小化，得到w和b得估计值，将y视为类后验概率p(y=1|x)，则上式写为
那么有：

简单综合，写成：
写成对数形式就是熟知的交叉熵损失函数了，这也是交叉熵损失的推导由来：

基于上式分别对w和b求偏导，

基于w和b得梯度下降对交叉熵损失函数最小化，相应得参数即为模型最优参数。
[这里具体推导一笔带过，详解机器学习入门笔记2 周志华.《机器学习》推导过程，本篇笔记着重梳理机器学习经典模型推导大致流程]

3.2 对数几率回归的代码实现

对数几率回归模型代码的编写思路

Numpy：

• 模型主体
  • Sigmoid
  • 模型公式
  • 交叉熵损失
  • 参数偏导
• 训练过程
  • 参数初始化
  • 多轮训练迭代过程
  • 梯度下降的参数优化更新
• 预测函数
  • 概率阈值
  • 可视化展示
  • 分类效果评估

sklearn：

• linear_model.LogisticRegression
• 使用范例

3.2.1 基于Numpy的对数几率回归实现

定义辅助函数

# 导入numpy模块
import numpy as np
### 定义Sigmoid函数
def sigmoid(x):
    '''
    输入：
    x：数组
    输出：
    z：经过sigmoid函数计算后的数组
    '''
    z = 1 / (1 + np.exp(-x))    
    return z

### 定义参数初始化数组 
def initialize_params(dims):
    '''
    输入：
    dims：参数维度
    输出：
    z：初始化后的参数向量W和参数值b
    '''
    # 将权重向量初始化为零向量
    W = np.zeros((dims, 1))
    # 将偏置初始化为零
    b = 0
    return W, b

定义对数几率回归模型主体

# 定义对数几率回归模型主体
def logistic(X, y, W, b):
    '''
    输入：
    X：输入特征矩阵
    y：输入标签向量
    W：权重系数
    b：偏置参数
    输出：
    a：对数几率回归模型输出
    cost：损失
    dW：权重梯度
    db：偏置梯度
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = X.shape[1]

    # 对数几率回归模型输出
    a = sigmoid(np.dot(X, W) + b)
    # 交叉熵损失
    cost = -1/num_train * np.sum(y*np.log(a) + (1-y)*np.log(1-a))
   
    # 权重梯度
    dW = np.dot(X.T, (a-y))/num_train
    # 偏置梯度
    db = np.sum(a-y)/num_train
    # 压缩损失数组维度
    cost = np.squeeze(cost) 

    return a, cost, dW, db

定义对数几率回归模型训练过程

### 定义对数几率回归模型训练过程
def logistic_train(X, y, learning_rate, epochs):    
    '''
    输入：
    X：输入特征矩阵
    y：输出标签向量
    learing_rate：学习率
    epochs：训练轮数
    输出：
    cost_list：损失列表
    params：模型参数
    grads：参数梯度
    '''
    # 初始化模型参数
    W, b = initialize_params(X.shape[1])  
    # 初始化损失列表
    cost_list = []  

    # 迭代训练
    for i in range(epochs):       
        # 计算当前次的模型计算结果、损失和参数梯度
        a, cost, dW, db = logistic(X, y, W, b)    
        # 参数更新
        W = W -learning_rate * dW
        b = b -learning_rate * db        

        # 记录损失
        if i % 100 == 0:
            cost_list.append(cost)   
        # 打印训练过程中的损失 
        if i % 100 == 0:
            print('epoch %d cost %f' % (i, cost)) 

    # 保存参数
    params = {            
        'W': W,            
        'b': b
    }        
    # 保存梯度
    grads = {            
        'dW': dW,            
        'db': db
    }           
    return cost_list, params, grads

定义预测函数

### 定义预测函数
def predict(X, params):
    '''
    输入：
    X：输入特征矩阵
    params：训练好的模型参数
    输出：y_pred：转换后的模型预测值
    '''
    
    # 模型预测值
    y_pred = sigmoid(np.dot(X, params['W']) + params['b']) 
    # 基于分类阈值对概率预测值进行类别转换
    for i in range(len(y_pred)):        
        if y_pred[i] > 0.5:
            y_pred[i] = 1
        else:
            y_pred[i] = 0
    return y_pred

生成模拟二分类数据集

# 导入matplotlib绘图库
import matplotlib.pyplot as plt
# 导入生成分类数据函数
from sklearn.datasets.samples_generator import make_classification
# 生成100X2的模拟二分类的数据集
X,labels=make_classification(n_samples=100, n_features=2, n_redundant=0, n_informative=2, random_state=1, n_clusters_per_class=2)
# 设置随机数种子
rng=np.random.RandomState(2)
# 对生成的特征数据添加一组均匀分布噪声
X+=2*rng.uniform(size=X.shape)
# 标签类别数
unique_lables=set(labels)
# 根据标签类别数设置颜色
colors=plt.cm.Spectral(np.linspace(0, 1, len(unique_lables)))
# 绘制模拟数据的散点图
for k, col in zip(unique_lables, colors):
    x_k=X[labels==k]
    plt.plot(x_k[:, 0], x_k[:, 1], 'o', markerfacecolor=col, markeredgecolor="k",
             markersize=14)
plt.title('data by make_classification()')
plt.show()

划分数据集

offset = int(X.shape[0] * 0.9)
X_train, y_train = X[:offset], labels[:offset]
X_test, y_test = X[offset:], labels[offset:]
y_train = y_train.reshape((-1,1))
y_test = y_test.reshape((-1,1))

print('X_train=', X_train.shape)
print('X_test=', X_test.shape)
print('y_train=', y_train.shape)
print('y_test=', y_test.shape)

模型训练和预测

# 执行对数几率回归模型训练
cost_list, params, grads = lr_train(X_train, y_train, 0.01, 1000)
# 打印训练好的模型参数
print(params)
# 基于训练参数对测试集进行预测
y_pred = predict(X_test, params)
print(y_pred)

测试集上的分类准确率评估

基于sklearn的分类评估方法衡量表现：

# 导入classification_report模块
from sklearn.metrics import classification_report
# 打印测试集分类预测评估报告
print(classification_report(y-test,y_pred))

绘制分类决策边界

### 绘制对数几率回归分类决策边界
def plot_logistic(X_train, y_train, params):
    '''
    输入：
    X_train:训练集输入
    y_train：训练集标签
    params：训练好的模型参数
    输出：分类决策边界图
    '''
    # 训练样本量
    n = X_train.shape[0]
    # 初始化类别坐标点列表
    xcord1 = []
    ycord1 = []
    xcord2 = []
    ycord2 = []    
    # 获取两类坐标点并存入列表
    for i in range(n):        
        if y_train[i] == 1:
            xcord1.append(X_train[i][0])
            ycord1.append(X_train[i][1])        
        else:
            xcord2.append(X_train[i][0])
            ycord2.append(X_train[i][1])
    # 创建绘图    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    # 绘制两类散点以不同颜色表示
    ax.scatter(xcord1, ycord1,s=32, c='red')
    ax.scatter(xcord2, ycord2, s=32, c='green')
    # 取值范围
    x = np.arange(-1.5, 3, 0.1)、
    # 分类决策边界公式
    y = (-params['b'] - params['W'][0] * x) / params['W'][1]
    # 绘图
    ax.plot(x, y)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

plot_logistic(X_train, y_train, params)

3.2.2基于sklearn的对数几率回归实现

LogisticRegression函数位于sklearn的linear_model模块下，定义该类的一个对数几率回归实例后，直接调用fit方法拟合训练集即可。

# 导入对数几率回归模块
from sklearn.linear_model import LogisticRegression
# 拟合训练集
clf = LogisticRegression(random_state=0).fit(X.train,y_train)
# 预测测试集
y_pred= clf.predict(X_test)
# 打印预测结果
print(y_pred)

3.3 小结

对数几率回归是用线性回归的结果来拟合真实标签的对数几率，可以将对数几率回归看成由条件概率分布表示的分类模型，对数几率回归是感知机模型、神经网络和支持向量机等模型的基础。

第 4 章 LASSO回归

4.1 LASSO回归原理推导

LASSO全称为 Least absolute shrinkage and selection operator ，译为 最小收缩与选择算子，LASSO回归是在线性回归的损失函数后加一个1-范数项（正则化项）。
$$L(w)=(y-wX{)}^2+\lambda ||w|{|}_1$$
，其中||w||即为矩阵的1-范数，λ为1-范数项的系数。

范数的概念

数学分析中，范数视为一种长度或距离概念的函数。针对向量或矩阵而言，常用范数包括0-范数、1-范数、2-范数和p-范数。

• 矩阵的0-范数表示矩阵中非零元素的个数
• 矩阵的1-范数表示矩阵中所有元素的绝对值之和
• 矩阵的2-范数表示矩阵中个元素的平方和再求均方根的结果

从机器学习角度，上述式子相当于给最初的线性回归损失函数添加了一个L1正则化项。从防止模型过拟合，正则化项相当于对目标参数施加了一个惩罚项，使得模型不过于复杂。在优化过程中，正则化的存在能够使那些不重要的特征系数逐渐为零，从而保留关键特征，使模型简化。

因此上述式子可简化，表示权重系数矩阵所有元素绝对值之和小于指定常数s，s取值越小。特征参数中被压缩到零的特征越多。
$$argmin(y-wX{)}^{2}s.t.\sum |w_{ij}|<s$$

由于L1正则化项的存在使得代价函数连续不可导，无法使用梯度下降法寻优，主要采用坐标下降法。

坐标下降法

坐标下降法使在当前坐标轴上搜索损失函数最小值，无需计算函数梯度， 而是沿着某一个坐标轴方向，通过一次一次的迭代更新权重系数的值，来渐渐逼近最优解。

坐标下降法属于一种非梯度优化的方法，它在每步迭代中沿一个坐标的方向进行线性搜索（线性搜索是不需要求导数的），通过循环使用不同的坐标方法来达到目标函数的局部极小值。

4.2 LASSO回归的代码实现

LASSO回归模型代码编写思路

Numpy：

• 模型主体
  • 回归模型公式
  • L1损失
  • 符号函数
  • 基于L1损失的梯度计算
• 训练过程
  • 参数初始化
  • 多轮训练迭代过程
  • 梯度下降的参数优化更新
• 数据测试
  • 测试结果
  • 可视化展示

sklearn：

• linear_model.Lasso
• 使用范例

4.2.1 基于Numpy的LASSO回归实现

定义符号函数

# 导入numpy模块
import numpy as np
### 定义符号函数
def sign(x):
    '''
    输入：
    X：浮点数值
    输出：
    整型数值
    '''
    if x > 0:
        return 1
    elif x < 0:
        return -1
    else:
        return 0
 
# 利用numpy对符号函数进行向量化
vec_sign = np.vectorize(sign)
vec_sign(np.zeros((3,1))) 

LASSO回归模型主体

### 定义LASSO回归损失函数
def l1_loss(X, y, w, b, alpha):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    w：变量参数权重矩阵
    b：偏置
    alpha：正则化系数
    输出：
    y_hat：线性模型预测输出
    loss：均方损失值
    dw：权重系数一阶偏导
    db：偏置一阶偏导
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = X.shape[1]
    # 回归模型预测输出
    y_hat = np.dot(X, w) + b
    # L1损失函数 
    loss = np.sum((y_hat-y)**2)/num_train + np.sum(alpha*abs(w))
    # 基于向量化符号函数的参数梯度计算
    dw = np.dot(X.T, (y_hat-y)) /num_train + alpha * vec_sign(w)
    db = np.sum((y_hat-y)) /num_train
    return y_hat, loss, dw, db

初始化模型参数

# 初始化模型参数
def initialize(dims):
    '''
    输入：
    dims：训练数据变量维度
    输出：
    w：初始化权重系数
    b：初始化偏置参数
    '''
    # 初始化权重系数为零向量
    w = np.zeros((dims, 1))
    # 初始化偏置系数为0
    b = 0
    return w, b

LASSO回归模型的训练过程

# 定义LASSO回归模型训练过程
def lasso_train(X, y, learning_rate=0.01, epochs=300):
    '''
    输入：
    X：输入变量矩阵
    y：输入标签向量
    learing_rate：学习率
    epochs：训练迭代次数
    输出：
    loss_his：每次迭代的L1损失列表
    params：优化后的参数字典
    grads：优化后的参数梯度字典
    '''
    # 记录训练损失的空列表
    loss_list = []
    # 初始化模型参数
    w, b = initialize(X.shape[1])
    # 迭代训练
    for i in range(1, epochs):
        # 记录当前迭代的预测值、损失和梯度
        y_hat, loss, dw, db = l1_loss(X, y, w, b, 0.1)
        # 基于梯度下降法的参数更新
        w += -learning_rate * dw
        b += -learning_rate * db
        loss_list.append(loss)
        
        # 每50次迭代打印当前损失信息
        if i % 50 == 0:
            print('epoch %d loss %f' % (i, loss))
        # 将当前迭代步优化后的参数保存在字典中 
        params = {
            'w': w,
            'b': b
        }
        # 将当前迭代步的梯度保存到字典中
        grads = {
            'dw': dw,
            'db': db
        }
    return loss, loss_list, params, grads

导入数据集

# 读取示例数据
data = np.genfromtxt('example.dat', delimiter=',')
# 选择特征与标签
x = data[:,0:100]
y = data[:,100].reshape(-1,1)
# 加一列
X = np.column_stack((np.ones((x.shape[0],1)),x))
# 划分训练集与测试集
X_train,y_train = X[:70],y[:70]
X_test,y_test = X[70:],y[70:]
print(X_train.shape,y_train.shape,X_test.shape,y_test.shape)

LASSO回归训练

# 执行训练示例
loss, loss_list, params, grads = lasso_train(X_train, y_train, 0.01, 300)
# 获取训练参数
params

定义预测函数

# 定义预测函数
def predict(X, params):
    w = params['w']
    b = params['b']
    
    y_pred = np.dot(X, w) + b
    return y_pred

y_pred = predict(X_test, params)
y_pred[:5]

4.2.2 基于sklearn的LASSO回归实现

# 导入线性模型模块
from sklearn import linear_model
# 创建lasso模型实例
sk_lasso = linear_model.Lasso(alpha=0.1)
# 对训练集进行拟合
sk_lasso.fit(X_train, y_train)
# 打印模型相关系数
print("sklearn Lasso intercept :", sk_lasso.intercept_)
print("\nsklearn Lasso coefficients :\n", sk_lasso.coef_)
print("\nsklearn Lasso number of iterations :", sk_lasso.n_iter_)

【近端梯度下降使Lasso线性回归问题求得最优解的证明，详见西瓜书p253】

第 5 章 Ridge回归

5.1 Ridge回归的原理推导

Ridge回归（岭回归）是一种使用2-范数作为惩罚项改造线性回归损失函数的模型。
$$L(w)=(y-wX{)}^2+\lambda ||w|{|}_2，其中\lambda ||w|{|}_2=\lambda \sum _{i=1}^n{w}_i^2，即L2正则化项$$
采用2-范数进行正则化的原理是最小化参数矩阵的每个元素，使其无限接近0但又不像L1那样等于0

5.2 Ridge回归的代码实现

5.2.1 基于Numpy的Ridge回归实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

划分数据集

# 划分数据集
data = pd.read_csv('./abalone.csv')
data['Sex'] = data['Sex'].map({'M':0, 'F':1, 'I':2})
X = data.drop(['Rings'], axis=1)
y = data[['Rings']]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)
X_train, X_test, y_train, y_test = X_train.values, X_test.values, y_train.values, y_test.values
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)

模型参数初始化

# 定义参数初始化函数
def initialize(dims):
    '''
    输入：
    dims：训练数据变量维度
    输出：
    w：初始化权重系数值
    b：初始化偏置参数
    '''
    # 初始化权重系数值为零向量
    w = np.zeros((dims, 1))
    # 初始化偏执参数为0
    b = 0
    return w, b

定义L2损失函数

# 定义ridge损失函数
def l2_loss(X, y, w, b, alpha):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    w：变量参数权重矩阵
    b：偏置
    alpha：正则化系数
    输出：
    y_hat：线性模型预测输出
    loss：均方损失
    dw：权重系数一阶偏导
    db：偏置一阶偏导
    '''
    # 训练样本量
    num_train = X.shape[0]
    # 训练特征数
    num_feature = X.shape[1]
    # 回归模型预测输出
    y_hat = np.dot(X, w) + b
    # L2损失函数
    loss = np.sum((y_hat-y)**2)/num_train + alpha*(np.sum(np.square(w)))
    # 参数梯度计算
    dw = np.dot(X.T, (y_hat-y)) /num_train + 2*alpha*w
    db = np.sum((y_hat-y)) /num_train
    return y_hat, loss, dw, db

定义Ridge回归模型得训练过程

# 定义Ridge回归模型得训练过程
def ridge_train(X, y, learning_rate=0.001, epochs=5000):
    '''
    输入：
    X：输入变量矩阵
    y：输出标签向量
    learing_rate：学习率
    epochs：训练迭代次数
    输出：
    loss_his：每次迭代得L1损失列表
    params：优化后的参数字典
    grads：优化后的参数梯度字典
    '''
    # 记录训练损失的空列表
    loss_list = []
    # 初始化模型参数
    w, b = initialize(X.shape[1])
    # 迭代训练
    for i in range(1, epochs):
        y_hat, loss, dw, db = l2_loss(X, y, w, b, 0.1)
        # 基于梯度下降的参数更新
        w += -learning_rate * dw
        b += -learning_rate * db
        # 记录当前迭代的损失
        loss_list.append(loss)
        
        # 每五十次迭代打印当前损失信息
        if i % 50 == 0:
            print('epoch %d loss %f' % (i, loss))
        # 将当前迭代步优化后的参数保存在字典中 
        params = {
            'w': w,
            'b': b
        }
        # 将当前迭代步的梯度保存在字典中
        grads = {
            'dw': dw,
            'db': db
        }
    return loss, loss_list, params, grads

执行Ridge回归模型训练

# 执行训练示例
loss, loss_list, params, grads = ridge_train(X_train, y_train, 0.01, 1000)

5.2.2 基于sklearn的Ridge回归实现

# 导入线性模型模块
from sklearn.linear_model import Ridge
# 创建Ridge模型实例
clf = Ridge(alpha=1.0)
# 对训练集进行拟合
clf.fit(X_train, y_train)
# 打印模型相关系数
print("sklearn Ridge intercept :", clf.intercept_)
print("\nsklearn Ridge coefficients :\n", clf.coef_)

5.3 小结

数学角度看，LASSO回归和Ridge回归都是在X^TX不可逆的情况下，通过给常规的平方损失函数添加L1正则化项和L2正则化项，回归问题有可行解。

业务可解释的角度看，影响一个变量的因素很多，但关键因素永远只会是少数。当影响一个变量的因素很多（特征数可能会大于样本量），LASSO回归和Ridge回归通过对损失函数施加正则化项的方式，，使回归建模过程中大量不重要的特征系数被压缩成0或接近0，从而找出对目标边框有较强影响的关键特征。

第 6 章 LDA线性判别分析

线性判别分析（Linear Discriminant Analysis，LDA）是一种经典的线性分类方法。

6.1 基本思想

LDA的基本思想是将数据投影到低维空间后，使得同一类数据尽可能接近，不同类数据尽可能疏远。所以，LDA是一种有监督的线性分类算法。

6.2 LDA数学推导

6.3 LDA的代码实现

LDA算法代码实现思路

Numpy：

• 数据标准化
• LDA流程
  • 数据分组
  • 均值和协方差计算
  • 计算类间散度
  • 计算均值差
  • 奇异值分解
  • 计算权重矩阵
  • 计算投影后的数据
• 数据测试

sklearn：

• sklearn.discriminant_analysis.LinearDiscriminantAnalysis

6.3.1 基于Numpy的LDA算法实现

LDA 算法实现

# 导入numpy库
import numpy as np

# 定义LDA类
class LDA():
    def __init__(self):
        # 初始化权重矩阵
        self.w = None
        
    # 协方差矩阵计算方法
    def calc_cov(self, X, Y=None):
        m = X.shape[0]
        # 数据标准化
        X = (X - np.mean(X, axis=0))/np.std(X, axis=0)
        Y = X if Y == None else (Y - np.mean(Y, axis=0))/np.std(Y, axis=0)
        return 1 / m * np.matmul(X.T, Y)
    
    # 数据投影方法
    def project(self, X, y):
        # LDA拟合获取模型权重
        self.fit(X, y)
        # 数据投影
        X_projection = X.dot(self.w)
        return X_projection
    
    # LDA拟合方法
    def fit(self, X, y):
        # 1. 按类分组
        X0 = X[y == 0]
        X1 = X[y == 1]

        # 2. 分别计算两类数据自变量的协方差矩阵
        sigma0 = self.calc_cov(X0)
        sigma1 = self.calc_cov(X1)
        # 3. 计算类内散度矩阵
        Sw = sigma0 + sigma1

        # 4. 分别计算两类数据自变量的均值和差
        u0, u1 = np.mean(X0, axis=0), np.mean(X1, axis=0)
        mean_diff = np.atleast_1d(u0 - u1)

        # 5. 对类内散度矩阵进行奇异值分解
        U, S, V = np.linalg.svd(Sw)
        # 6. 计算类内散度矩阵的逆
        Sw_ = np.dot(np.dot(V.T, np.linalg.pinv(np.diag(S))), U.T)
        # 7. 计算w
        self.w = Sw_.dot(mean_diff)

    
    # LDA分类预测
    def predict(self, X):
        # 初始化预测结果为空列表
        y_pred = []
        # 遍历待预测样本
        for sample in X:
            # 模型预测
            h = sample.dot(self.w)
            y = 1 * (h < 0)
            y_pred.append(y)
        return y_pred

LDA算法数据集测试

# 导入库
from sklearn import datasets
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

# 导入iris数据集
data = datasets.load_iris()
# 数据
X = data.data
#标签
y = data.target
# 取标签不为2的数据
X = X[y != 2]
y = y[y != 2]
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=41)
print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)
# 创建LDA模型实例
lda = LDA()
# LDA模型拟合
lda.fit(X_train, y_train)
# LDA模型预测
y_pred = lda.predict(X_test)

# 测试集上的分类准确率
from sklearn.metrics import accuracy_score
accuracy = accuracy_score(y_test, y_pred)
print(accuracy)

6.3.2 基于sklearn的LDA算法实现

# 导入模块
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
# 创建LDA分类器
clf = LinearDiscriminantAnalysis()
# 模型拟合
clf.fit(X_train, y_train)
# 模型预测
y_pred = clf.predict(X_test)
# 测试集上的分类准确率
accuracy = accuracy_score(y_test, y_pred)
print(accuracy)

6.4 小结

• LDA基本思想是将数据投影到低维空间后，使得同一类数据尽可能接近，不同类数据尽可能疏远。
• 数学角度看，使训练样本的类内散度尽可能小，而类间散度尽可能大，从而设计出LDA的优化目标
• 多分类LDA将样本投影到低维空间，降低数据集的原有维度，并在投影过程中使用类别信息，他是经典的监督降维技术

• 第 8 章 决策树

决策树基于特征对数据实例按照条件不断进行划分，最终达到分类或回归的目的，决策树模型的预测过程既可以看作一组if-then条件的集合，也可以视作定义在特征空间与类空间中的条件概率分布。

决策树模型的核心概念包括特征选择方法、决策树构造过程和决策树的剪枝。常见的特征选择方法包括信息增益、信息增益比和基尼指数，对应的三种常见的决策树算法为ID3、C4.5和CART。

8.1 决策树

决策树通过树形结构对数据进行分类。一棵完整的结构树由结点和有向边构成，其中内部节点表示特征，叶子结点表示类别。决策树从根结点开始，选取数据中某一特征，根据特征取值对实例进行分配，通过不断选取特征进行实例分配，决策树可以达到对所有实例进行分类的目的。

给定训练集D={(x₁,y₁), (x₂,y₂), …, (x_N, y_N)}，其中x_i=(x_i⁽¹⁾, x_i⁽²⁾,…,x_i⁽ⁿ⁾)^T表示输入实例，即特征向量，n为特征个数，y_i为类别标记，N为样本容量。决策树的学习目标是构造一个决策树模型，对输入实例进行最大可能的正确分类。

决策树模型没有偏离机器学习模型训练的一般范式，即正则化参数的同时最小化经验误差。假设树T的叶子结点个数为|T|，t为树T的叶子结点，每个叶子结点有N_t个样本，假设k类样本有N_tk个，其中k=1,2,…,K，H_t(T)we为叶子结点上的经验熵，α≥0为正则化参数，那么决策树学习的损失函数可表示为：
$$L_{\alpha }(T)=\sum _{t=1}^{|T|}{N}_t{H}_t(T)+\alpha |T|$$
具体而言，递归选择最优特征，按照该特征取值将训练集划分成不同子集，使得各子集有一个在当前条件下的最优分类，如果这些子集能被正确分类，即可将这些子集归类到叶子结点，否则从这些子集中继续选取最优特征，如此递归下去直到所有子集被正确分类。但该方法可能会造成过拟合，因此需要对决策树进行剪枝。

因此，决策树模型包括特征选择、决策树构建、决策树剪枝。

8.2 特征选择：从信息增益到基尼指数

8.2.1 特征选择

特征对数据集分类能力低下情况：一个特征数据集进行分类的效果与随机选取的分类效果并无差异

特征对数据集分类能力较强情况：一个特征数据集进行分类后的分支结点尽可能属于同一类别，即该结点具有较高的纯度

特征选择：从数据集中选择具备较强分类能力的特征来进行数据集划分

选择方式：信息增益、信息增益比、基尼指数

8.2.2 信息增益

熵：信息论和概率统计中，熵是一种描述随机变量不确定性的度量方式，也可以用来描述样本集合的纯度，信息熵越低，样本的不确定性越小，样本的纯度越高。

假设当前样本数据集D中第k个类别所占比例为p_k(k=1, 2, …, Y)，那么该样本数据集的熵为：
$$E(D)=-\sum _{k=1}^Y{p}_{k}log{p}_k$$
假设离散随机变量(X, Y)的联合概率分布为：
$$P(X=x_i,Y=y_j)=p_{ij}(i=1,2,...,m,j=1,2,...,n)$$
条件熵E(Y|X)表示在已知随机变量X的条件下Y的不确定性的度量，E(Y|X)定义为在给定X的条件下Y的条件概率分布的熵对X的数学期望。条件熵可以表示为：
$$E(Y|X)=\sum _{i=1}^m{p}_{i}E(Y|X=x_i)$$
其中p_i=P(X = X_i)，i=1, 2, …, n。在实际数据计算中，熵和条件熵的概率是基于极大似然估计得到，对应的熵和条件熵称经验熵和经验条件熵。

信息增益

指的是由于得到特征X的信息而是类Y的信息不确定性减少的程度，即信息增益是一种描述目标类别确定性增加的量，特征的信息增益越大，目标的确定性越大。

假设训练集D的经验熵为E(D)，给定特征A的条件下D的经验条件熵为E(D|A)，那么信息增益定义为经验熵E(D)与经验条件熵E(D|A)之差：
$$g(D,A)=E(D)-E(D|A)$$
给定训练集D和特征A，经验熵E(D)表示为对数据集D进行分类的不确定性，经验条件熵E(D|A)表示在给定特征A之后对数据集D进行分类的不确定性，二者的差是两个不确定性之间的差，即信息增益。

具体到数据集D中，每个特征一般有不同的信息增益，信息增益越大，代表对应特征分类能力越强

[ID3算法基于信息增益进行特征选择]

信息熵计算定义

# 导入numpy库
import numpy as np
# 导入对数计算模块
from math import log
# 定义信息熵计算函数
def entropy(ele):
    '''
    输入：
    ele：包含类别取值的列表
    输出；
    信息熵值
    '''
    # 计算列表取值的概率分布
    probs = [ele.count(i)/len(ele) for i in set(ele)]
    # 计算信息熵
    entropy = -sum([prob*log(prob,2) for prob in probs])
    return entropy

信息增益计算（天气特征对于高尔夫数据集的信息增益）

# 导入pandas库
import pandas as pd
# 以数据框结构读取高尔夫数据集
df = pd.read_csv('./golf_data.csv')
# 计算数据集的经验熵
# 以'play'为目标变量，即是否打高尔夫
entropy_D = entropy(df['play'].tolist())
# 计算天气特征的经验条件熵
# 其中subset1~subset3为根据天气特征三个取值划分后的子集
entropy_DA = len(subset1)/len(df)*entropy(subset1['play'].tolist()) +
			 len(subset2)/len(df)*entropy(subset2['play'].tolist()) +
             len(subset3)/len(df)*entropy(subset3['play'].tolist()) 
# 计算天气特征的信息增益
info_gain = entropy_D - entropy_DA
print('天气特征对于数据集分类的信息增益为：', info_gain)

8.2.3 信息增益比

信息增益存在问题：当某个特征分类取值较多时，该特征的信息增益计算结果就比较大，因此基于信息增益选择特征时，会偏向于取值较大的特征。

使用信息增益比对上述问题进行校正。特征A对数据集D的信息增益比定义为其信息增益g(D, A)与数据集D关于特征A取值的熵E_A(D)的比值：
$$g_R(D,A)=\frac{g(D,A)}{E_A(D)}，其中E_A(D)=-\sum _{i=1}^n\frac{|D_i|}{|D|}lo{g}_2\frac{D_i}{D}，n表示特征A的取值个数。$$
跟信息增益一样，在基于信息增益比进行特征选择时，选择信息增益比最大的特征作为决策树分裂结点。

C4.5算法基于信息增益比进行特征选择

8.2.4 基尼指数

假设样本有K类，样本属于第k类的概率为p_k，则样本类别概率分布的基尼指数定义为：
$$Gini(p)=\sum _{k=1}^K{p}_k(1-p_k)=1-\sum _{k=1}^K{p}_k^2$$
如果给定训练集D，C_k是属于第k类样本的集合，则该训练集的基尼指数定义为：
$$Gini(D)=1-\sum _{k=1}^K(\frac{|C_k|}{|D|}{)}^2$$
如果训练集D根据特征A某一取值a划分为D₁和D₂两个部分，那么在特征A这个条件下，训练集D的基尼指数定义为：
$$Gini(D,A)=\frac{D_1}{D}Gini(D_1)+\frac{D_2}{D}Gini(D_2)$$
训练集D的基尼指数Gini(D)表示该集合的不确定性，Gini(D, A)表示训练集D经过A = a划分后的不确定性。，Gini(D, A)越小，训练集的不确定性越小，对应的特征对训练样本的分类能力越强。

CART算法基于基尼指数进行特征选择

基尼指数计算

# 导入numpy库
import numpy as np
# 定义基尼指数计算函数
def gini(nums):
    '''
    输入：
    nums：包含类别取值的列表
    输出：基尼指数值
    '''
    # 获取列表类别的概率分布
    probs = [nums.count(i)/len(nums) for i in set(nums)]
    # 计算基尼指数
    gini = sum([p*(1-p) for p in probs])
    return gini

天气特征条件下的基尼指数

# 计算天气特征的基尼指数
# 导入pandas库
import pandas as pd
# 以数据框结构读取高尔夫数据集
df = pd.read_csv('./golf_data.csv')
# 其中subset1和subset2为根据天气特征取值为晴或者非晴的划分的两个子集
gini_DA = len(subset1)/len(df) * gini(subset1['play'].tolist()) +
          len(subset2)/len(df) * gini(subset2['play'].tolist())
print('天气特征取值为晴的基尼指数为', gini_DA)    

8.3 决策树模型

8.3.1 ID3

ID3算法的全称为Iterative Dichotomiser3，即三代迭代二叉树，其核心是基于信息增益递归地选择最优特征构造决策树。

具体方法

• 预设决策树根结点
• 对所有特征计算信息增益，选择信息增益最大的作为最优特征，根据该特征的不同取值建立子结点
• 对每个结点递归调用上述方法，直到信息增益很小或者没有特征可选时，即可构建最终的ID3决策树

给定训练集D、特征集合A以及信息增益阈值δ，ID3算法流程

• step1 如果D中所有实例属于同一类别C_k，那么所构建的决策树T为单结点树，并且类C_k即为该结点的类的标记
• step2 如果T不是单结点树，则计算特征集合A中各特征对D的信息增益，选择信息增益最大的特征A_g
• step3 如果A_g的信息增益小于阈值δ，则将T视为单结点树，并将D中所属数量最多的类C_K作为该结点的类的标记并返回T
• step4 否则，可对A_g的每一特征取值a_i，按照A_g=a_i将D划分为若干非空子集D_i，以D_i所属数量最多的类作为标记并构建子结点，由结点和子结点构成树T并返回
• step5 对第i个子结点，以D_i为训练集，以A-A_g为特征集，递归调用step1 – step4，即可得到决策子树T_i并返回

数据集划分函数

# 根据数据集和指定特征定义数据集划分函数
def df_split(df, col):
    '''
    输入：
    df：待划分的训练数据
    col：划分数据的依据特征
    输出：
    res_dict：根据特征取值划分后的不同数据集字典
    '''
    # 获取依据特征的不同取值
    unique_col_val = df[col].unique()
    # 创建划分结果的数据框字典
    res_dict = {elem : pd.DataFrame for elem in unique_col_val}
    # 依据特征取值进行划分
    for key in res_dict.keys():
        res_dict[key] = df[;][df[col] == key]
    return res_dict    

选择最优特征

# 根据训练集和标签选择信息增益最大的特征作为最优特征
def choose_best_feature(df, label):    
    '''
    输入：
    df：待划分的训练数据
    label：训练标签
    输出：
    max_value：最大信息增益
    best_feature：最优特征
    max_splited：根据最优特征划分后的数据字典
    '''
    # 计算训练标签的信息熵
    entropy_D = entropy(df[label].tolist())    
    # 特征集
    cols = [col for col in df.columns if col not in [label]]    
    # 初始化最大信息增益、最优特征和划分后的数据集
    max_value, best_col = -999, None
    max_splited = None
    # 遍历特征并根据特征取值进行划分
    for col in cols:
        # 根据当前特征划分后的数据集
        splited_set = split_dataframe(df, col)
        # 初始化经验条件熵
        entropy_DA = 0
        # 对划分后的数据集遍历计算
        for subset_col, subset in splited_set.items():            
            # 计算划分后的数据子集的标签信息熵
            entropy_Di = entropy(subset[label].tolist())            
            # 计算当前特征的经验条件熵
            entropy_DA += len(subset)/len(df) * entropy_Di        
        # 计算当前特征的信息增益
        info_gain = entropy_D - entropy_DA        
        # 获取最大信息增益
        if info_gain > max_value:
            max_value, best_col = info_gain, col
            max_splited = splited_set    
        return max_value, best_col, max_splited

构建ID3决策树

# ID3算法类
class ID3Tree:    
    # 定义决策树结点类
    class TrNode:
        # 定义树结点
        def __init__(self, name):
            self.name = name
            self.connections = {}    
         # 定义树连接   
        def connect(self, label, node):
            self.connections[label] = node    
    # 定义全局变量，包括数据集、特征集、标签和根结点
    def __init__(self, data, label):
        self.columns = data.columns
        self.data = data
        self.label = label
        self.root = self.Node("Root")    
    
    # 构建树的调用
    def construct_tree(self):
        self.construct(self.root, "", self.data, self.columns)    
    
    # 决策树构建方法
    def construct(self, parent_node, parent_label, sub_df, columns):
        # 选择最优特征
        max_value, best_feature, max_splited = choose_best_feature(sub_df[columns], self.label)        
        # 如果选不到最优特征，则构造单结点树
        if not best_feature:
            node = self.Node(input_data[self.label].iloc[0])
            parent_node.connect(parent_connection_label, node)            
        return

        # 根据最优特征以及子结点构建树
        node = self.Node(best_feature)
        parent_node.connect(parent_label, node)

        # 以A-Ag为新的特征集
        new_columns = [col for col in columns if col != best_feature]        
        # 递归的构建决策树
        for splited_value, splited_data in max_splited.items():
            self.construct(node, splited_value, splited_data, new_columns)

基于高尔夫数据集的ID3决策树

# 读取高尔夫数据集
df = pd.read_csv('./example_data.csv')
# 创建ID3决策树实例
id3_tree = ID3Tree(df,'play')
# 构造ID3决策树
id3_tree.construct_tree()

基于sklearn的ID3决策树构建

from sklearn.datasets import load_iris
from sklearn import tree
import graphviz

iris = load_iris()
# criterion选择entropy，这里表示选择ID3算法
clf = tree.DecisionTreeClassifier(criterion='entropy', splitter='best')
clf = clf.fit(iris.data, iris.target)

dot_data = tree.export_graphviz(clf, out_file=None,
                               feature_names=iris.feature_names,
                               class_names=iris.target_names,
                               filled=True, 
                               rounded=True,
                               special_characters=True)
graph = graphviz.Source(dot_data)
graph

8.3.2 C4.5

C4.5算法整体与ID3类似，不同之处在于C4.5在构造决策树时使用信息增益比作为特征选择方法。

给定训练集D、特征集合A以及信息增益阈值δ，C4.5算法流程

• step1 如果D中所有实例属于同一类别C_k，那么所构建的决策树T为单结点树，并且类C_k即为该结点的类的标记
• step2 如果T不是单结点树，则计算特征集合A中各特征对D的信息增益比，选择信息增益比最大的特征A_g
• step3 如果A_g的信息增益比小于阈值δ，则将T视为单结点树 ，并将D中所属数量最多的类C_k作为该结点的类的标记并返回T
• step4 否则，可对A_g的每一特征值a_i，按照A_g=a_i将D划分为若干非空子集D_i，以D_i中所属数量最多的类作为标记并构建子结点，由结点和子结点构成树T并返回
• step5 对第i个子结点，以D_i为训练集，以A-A_g为特征集，递归调用step1 – step2，即可得到决策树子树T_i并返回

8.3.3 CART分类树

CART算法全称分类与回归树（classification and regression tree），既可以用于分类，也可以用于回归，CART算法特征选择方法基于基尼指数，CART算法不仅仅包括决策树生成算法，还包括决策树剪枝算法。

CART 分类树生成算法

输入：训练数据集D，停止计算的条件

输出：CART决策树

根据训练数据集，从根结点开始，递归地对每个结点进行以下操作，构建二叉决策树：

• step1 设结点的训练数据集为D，计算现有特征对该数据集的基尼指数，此时，对每个特征A，对其可能取的每个值a，根据样本A=a的测试为“是”或“否”，将D分割成D₁和D₂两部分，基于下式计算A=a基尼指数。

$$Gini(D,A)=\frac{D_1}{D}Gini(D_1)+\frac{D_2}{D}Gini(D_2)$$

• step2 取基尼指数最小的特征及其对应的划分点作为最优特征和最优划分点，据此将当前结点划分为两个子结点，将训练集根据特征分配到两个子结点中。
• step3 对两个子结点递归调用step1 – step2，直至满足条件
• 最后即可生成CART分类决策树

8.3.4 CART回归树

假设训练输入X和输出Y，给定训练集D={(x₁, y₁), (x₂, y₂), …, (x_N, y_N)}，CART回归树建立思路如下。

回归树对应特征空间的一个划分以及在该划分单元上的输出值。假设特征空间有M个划分单元R₁，R₂，…，R_M，且每个划分单元都有一个输出权重c_m，那么回归树模型可以表示为：
$$f(x)=\sum _{m=1}^M{c}_{m}I(x\in R_m)(1)$$
根线性回归模型一样，回归树模型训练的目的同样是最小化均方误差，以求得最优输出权重$ \hat{c} $m。用平方误差最小方法求解每个单元上的最优权重，最优输出权重$ \hat{c}_m$可以通过每个单元上所有输入实例x_i对应的输出值y_i来确定，即：
$${\hat{c}}_m=average(y_i|x_i\in R_m)(2)$$

假设随机选取第j个特征x^(j)及其对应的某个取值s，将其作为划分特征和划分点，同时定义两个区域：
R 1 ( j , s ) = { x ∣ x ( j ) ≤ s } ; R 2 ( j , s ) = { x ∣ x ( j ) > s } ( 3 ) R_1(j,s)=\{{x|x^{(j)}}≤s\};R_2(j,s)=\{x|x^{(j)}>s\}\quad(3) R1​(j,s)={x∣x(j)≤s};R2​(j,s)={x∣x(j)>s}(3)

然后求解：
$$mi{n}_{js}[mi{n}_{c1}\sum _{x_i\in R_1(j,s)}(y_i-c_1{)}^2+mi{n}_{c2}\sum _{x_i\in R_2(j,s)}(y_i-c_2{)}^2](4)$$
求解式(4)即可得到输入特征j和最优划分点s。按照上述平方误差最小准则球的全局最优特征和取值，并据此将特征空间划分成两个子区域，对每个子区域重复前述划分过程，直至满足停止条件，即可生成一颗回归树。

CART回归树生成算法描述

• step1 根据式(4)求解最优特征j和最优划分点s，遍历训练集所有特征，对固定划分特征扫描划分点s，求得式(4)最小值。
• step2 通过式(2)和式(3)确定最优的(j, s)来划分特征空间区域并决定相应的输出权重
• step3 对划分的两个子区域递归调用step1 – step2，直至满足条件
• step4 将特征空间划分为M个单元R₁，R₂，…，R_M，生成回归树

$$f(x)=\sum _{m=1}^M{\hat{c}}_{m}I(x\in R_m)(5)$$

8.3.5 CART算法实现

8.3.5.1 CART算法实现思路

无论是分类树还是回归树，二者对于树结点和基础二叉树的实现方式是一致的，差异在于特征选择方法和叶子结点取值预测方法

分类树

• 树结点
  • 分裂特征
  • 特征分裂阈值
  • 叶子结点值
  • 左子树
  • 右子树
• 二叉树
  • 树基本属性
  • 树构建流程
  • 预测函数
• 分类树
  • 基尼不纯度
  • 多数投票
  • 拟合方法

回归树

• 树结点
• 二叉树
• 回归树
  • 方差减少量
  • 结点平均值
  • 拟合方法

辅助函数

• 特征分裂
• 基尼指数计算
• 其他

8.3.5.2 二叉决策树

导入模块

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, mean_squared_error
from utils import feature_split, calculate_gini

utils

import numpy as np

### 定义二叉特征分裂函数
def feature_split(X, feature_i, threshold):
    split_func = None
    if isinstance(threshold, int) or isinstance(threshold, float):
        split_func = lambda sample: sample[feature_i] >= threshold
    else:
        split_func = lambda sample: sample[feature_i] == threshold

    X_left = np.array([sample for sample in X if split_func(sample)])
    X_right = np.array([sample for sample in X if not split_func(sample)])

    return np.array([X_left, X_right])


### 计算基尼指数
def calculate_gini(y):
    # 将数组转化为列表
    y = y.tolist()
    probs = [y.count(i)/len(y) for i in np.unique(y)]
    gini = sum([p*(1-p) for p in probs])
    return gini

定义树结点

### 定义树结点
class TreeNode():
    def __init__(self, feature_i=None, threshold=None,
                 leaf_value=None, left_branch=None, right_branch=None):
        # 特征索引
        self.feature_i = feature_i          
        # 特征划分阈值
        self.threshold = threshold 
        # 叶子节点取值
        self.leaf_value = leaf_value   
        # 左子树
        self.left_branch = left_branch     
        # 右子树
        self.right_branch = right_branch 

定义基础的二叉决策树

### 定义二叉决策树
class BinaryDecisionTree(object):
    ### 决策树初始参数
    def __init__(self, min_samples_split=2, min_gini_impurity=999,
                 max_depth=float("inf"), loss=None):
        # 根结点
        self.root = None  
        # 节点最小分裂样本数
        self.min_samples_split = min_samples_split
        # 节点初始化基尼不纯度
        self.mini_gini_impurity = min_gini_impurity
        # 树最大深度
        self.max_depth = max_depth
        # 基尼不纯度计算函数
        self.gini_impurity_calculation = None
        # 叶子节点值预测函数
        self._leaf_value_calculation = None
        # 损失函数
        self.loss = loss

    ### 决策树拟合函数
    def fit(self, X, y, loss=None):
        # 递归构建决策树
        self.root = self._build_tree(X, y)
        self.loss=None

    ### 决策树构建函数
    def _build_tree(self, X, y, current_depth=0):
        # 初始化最小基尼不纯度
        init_gini_impurity = 999
        # 初始化最优特征索引和阈值
        best_criteria = None    
        # 初始化数据子集
        best_sets = None        

        # 合并输入和标签
        Xy = np.concatenate((X, y), axis=1)
        # 获取样本数和特征数
        n_samples, n_features = X.shape
        # 设定决策树构建条件
        # 训练样本数量大于节点最小分裂样本数且当前树深度小于最大深度
        if n_samples >= self.min_samples_split and current_depth <= self.max_depth:
            # 遍历计算每个特征的基尼不纯度
            for feature_i in range(n_features):
                # 获取第i特征的所有取值
                feature_values = np.expand_dims(X[:, feature_i], axis=1)
                # 获取第i个特征的唯一取值
                unique_values = np.unique(feature_values)

                # 遍历取值并寻找最佳特征分裂阈值
                for threshold in unique_values:
                    # 特征节点二叉分裂
                    Xy1, Xy2 = feature_split(Xy, feature_i, threshold)
                    # 如果分裂后的子集大小都不为0
                    if len(Xy1) > 0 and len(Xy2) > 0:
                        # 获取两个子集的标签值
                        y1 = Xy1[:, n_features:]
                        y2 = Xy2[:, n_features:]

                        # 计算基尼不纯度
                        impurity = self.impurity_calculation(y, y1, y2)

                        # 获取最小基尼不纯度
                        # 最佳特征索引和分裂阈值
                        if impurity < init_gini_impurity:
                            init_gini_impurity = impurity
                            best_criteria = {"feature_i": feature_i, "threshold": threshold}
                            best_sets = {
                                "leftX": Xy1[:, :n_features],   
                                "lefty": Xy1[:, n_features:],   
                                "rightX": Xy2[:, :n_features],  
                                "righty": Xy2[:, n_features:]   
                                }
        
        # 如果计算的最小不纯度小于设定的最小不纯度
        if init_gini_impurity < self.mini_gini_impurity:
            # 分别构建左右子树
            left_branch = self._build_tree(best_sets["leftX"], best_sets["lefty"], current_depth + 1)
            right_branch = self._build_tree(best_sets["rightX"], best_sets["righty"], current_depth + 1)
            return TreeNode(feature_i=best_criteria["feature_i"], threshold=best_criteria[
                                "threshold"], left_branch=left_branch, right_branch=right_branch)

        # 计算叶子计算取值
        leaf_value = self._leaf_value_calculation(y)
        return TreeNode(leaf_value=leaf_value)

    ### 定义二叉树值预测函数
    def predict_value(self, x, tree=None):
        if tree is None:
            tree = self.root

        # 如果叶子节点已有值，则直接返回已有值
        if tree.leaf_value is not None:
            return tree.leaf_value

        # 选择特征并获取特征值
        feature_value = x[tree.feature_i]

        # 判断落入左子树还是右子树
        branch = tree.right_branch
        if isinstance(feature_value, int) or isinstance(feature_value, float):
            if feature_value >= tree.threshold:
                branch = tree.left_branch
        elif feature_value == tree.threshold:
            branch = tree.left_branch

        # 测试子集
        return self.predict_value(x, branch)

    ### 数据集预测函数
    def predict(self, X):
        y_pred = [self.predict_value(sample) for sample in X]
        return y_pred

8.3.5.3 分类树

分类树实现

### CART分类树
class ClassificationTree(BinaryDecisionTree):
    ### 定义基尼不纯度计算过程
    def _calculate_gini_impurity(self, y, y1, y2):
        p = len(y1) / len(y)
        gini = calculate_gini(y)
        gini_impurity = p * calculate_gini(y1) + (1-p) * calculate_gini(y2)
        return gini_impurity
    
    ### 多数投票
    def _majority_vote(self, y):
        most_common = None
        max_count = 0
        for label in np.unique(y):
            # 统计多数
            count = len(y[y == label])
            if count > max_count:
                most_common = label
                max_count = count
        return most_common
    
    # 分类树拟合
    def fit(self, X, y):
        self.impurity_calculation = self._calculate_gini_impurity
        self._leaf_value_calculation = self._majority_vote
        super(ClassificationTree, self).fit(X, y)

分类树测试

# 导入数据集
from sklearn import datasets
# 导入数据划分模块
from sklearn.model_selection import train_test_split
# 导入准确率评估函数
from sklearn.metrics import accuracy_score
# 导入iris数据集
data = datasets.load_iris()
# 获取输入和标签
X, y = data.data, data.target
# 注意！是否要对y进行reshape取决于numpy版本
y = y.reshape(-1,1)
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
# 创建分类树模型实例
clf = ClassificationTree()
# 分类树训练
clf.fit(X_train, y_train)
# 分类树预测
y_pred = clf.predict(X_test)
# 打印模型分类准确率
print(accuracy_score(y_test, y_pred))

sklearn分类树测试

# 导入分类树模块
from sklearn.tree import DecisionTreeClassifier
# 创建分类树实例
clf = DecisionTreeClassifier()
# 分类树训练
clf.fit(X_train, y_train)
# 分类树预测
y_pred = clf.predict(X_test)

print(accuracy_score(y_test, y_pred))

8.3.5.4 回归树

回归树实现

### CART回归树
class RegressionTree(BinaryDecisionTree):
    # 计算方差减少量
    def _calculate_variance_reduction(self, y, y1, y2):
        var_tot = np.var(y, axis=0)
        var_y1 = np.var(y1, axis=0)
        var_y2 = np.var(y2, axis=0)
        frac_1 = len(y1) / len(y)
        frac_2 = len(y2) / len(y)
        # 计算方差减少量
        variance_reduction = var_tot - (frac_1 * var_y1 + frac_2 * var_y2)
        
        return sum(variance_reduction)

    # 结点值取平均
    def _mean_of_y(self, y):
        value = np.mean(y, axis=0)
        return value if len(value) > 1 else value[0]

    # 回归树拟合
    def fit(self, X, y):
        self.impurity_calculation = self._calculate_variance_reduction
        self._leaf_value_calculation = self._mean_of_y
        super(RegressionTree, self).fit(X, y)

回归树测试

# 导入波士顿房价数据集模块
from sklearn.datasets import load_boston
# 导入均方误差评估函数
from sklearn.metrics import mean_suqared_error
# 获取输入和标签
X, y = load_boston(return_X_y=True)
y = y.reshape(-1,1)
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
# 创建回归树模型实例
model = RegressionTree()
# 模型训练
model.fit(X_train, y_train)
# 模型预测
y_pred = model.predict(X_test)
# 评估均方误差
mse = mean_squared_error(y_test, y_pred)

print("Mean Squared Error:", mse)

sklearn回归树测试

# 导入回归树模块
from sklearn.tree import DecisionTreeRegressor
# 创建回归树模型实例
reg = DecisionTreeRegressor()
reg.fit(X_train, y_train)
y_pred = reg.predict(X_test)
mse = mean_squared_error(y_test, y_pred)

print("Mean Squared Error:", mse)

8.4 决策树剪枝

剪枝原因：决策树生成算法递归产生决策树，生成的决策树大而全，容易导致过拟合

决策树剪枝一般包括预剪枝和后剪枝。

预剪枝

在决策树生成过程中提前停止决策树的增长的一种剪枝方法。

主要思路：在决策树结点分裂之前，计算当前结点划分能否提升模型泛化能力，如果不能，则决策树在该结点停止生长。

缺点：预剪枝提前停止树生长的方法，一定程度上存在欠拟合的风险，导致决策树生长不够完全。

后剪枝

后剪枝主要通过极小化决策树整体损失函数来实现。决策树学习的目标是最小化损失函数$L_{\alpha }(T)={\sum }_{t=1}^{|T|}{N}_t{H}_t(T)+\alpha |T|（设k类样本有N_{tk}个，其中k=1,2,...,K，H_t(T)为叶子结点上的经验熵，\alpha \ge 0为正则化参数）$

上损失函数中的经验熵可表示为：
$$H_t(T)=-\sum _k\frac{N_{tk}}{N_t}log\frac{N_{tk}}{N_t}$$
令损失函数第一项为：
$$L(T)=\sum _{t=1}^{|T|}{N}_t{H}_t(T)=-\sum _{t=1}^{|T|}\sum _{k=1}^{K}log\frac{N_{tk}}{N_t}$$
因此，损失函数改写为：
$$L_{\alpha }(T)=L(T)+\alpha |T|$$
其中L(T)为模型的经验误差项，|T|表示决策树复杂度，α≥0即为正则化参数，用于控制经验误差项和正则化项之间的影响。

决策树后剪枝，是在复杂度α确定下，选择损失函数L_α(T)最小的决策树模型。给定生成算法得到的决策树T和正则化参数α。

剪枝算法描述

• step1 计算每个树结点的经验熵H_t(T)
• step2 递归地自底向上回缩，假设一组叶子节点回缩到父结点前后的树分别为T_before与T_after，其对应的损失函数分别为L_α(T_before)和Lα(T_after)，如果L_α(T_after)≤L_α(T_before)，则进行剪枝，将父结点变成新的叶子结点
• step3 重复step2，直到得到损失函数最小的子树T_α。

CART算法剪枝属于后剪枝法。CART后剪枝先通过计算子树的损失函数来实现剪枝并得到一个子树序列，然后通过交叉验证的方法从子树序列中选取最优子树。

8.5 小结

• 决策树一般包括特征选择条件、决策树生成算法、决策树剪枝算法
• 常用的基础决策树算法包括ID3算法[信息增益最大]、C4.5算法[信息增益比最大]、CART算法[基尼指数最小]
• 由于生成的决策树可能过度增长，会造成过拟合，因此需要进行剪枝（包括预剪枝和后剪枝）

第 9 章 神经网络

神经网络（neutral network）可以溯源到原先的感知机（perceptron），单层感知机逐渐发展到多层感知机，加入隐藏层使得感知机发展成能拟合一切的神经网络，反向传播算法是整个神经网络训练的核心。

9.1 感知机

9.1.1 感知机推导

感知机是一个线性模型，旨在建立一个线性分隔超平面对线性可分的数据集进行分类。

模型接收$x_1、x_2、...、x_n$多个输入，将输入与权重系数$w$进行加权求和并经过Sigmoid函数进行激活，将激活结果$y$作为输出。执行完向前计算得到输出后，模型需要根据输出和实际输出按照损失函数计算当前损失，计算损失函数关于权重和偏置使得损失最小。

给定输入实例$x\in X$，输出 y ∈ Y = { + 1 , − 1 } y∈Y=\{+1,-1\} y∈Y={+1,−1}，由输入到输出的感知机模型表示为：
$$y=sign(w\cdot x+b)，其中w为权重系数，b偏置系数，sign为符号函数$$
感知机学习目标是建立一个线性分隔超平面，以将训练数据正例和负例完全分开，可以提供最小化损失函数来确定模型参数$w$和$b$。

假设输入空间中任意一点$x_0$到线性分隔超平面的距离：
$$\frac{1}{||w||}w\cdot x_0+b，其中||w||为w的2-范数。$$
对于任意一误分类点$(x_i,y_i)$，当$w\cdot x_i+b>0$时，$y_i=-1$;当$w\cdot x_i+b<0$时，$y_i=+1$，因而都有$-y_i(w\cdot x_i+b)$成立。所以误分类点到线性分隔超平面的距离$S$为：
$$-\frac{1}{||w||}{y}_i(w\cdot x_0+b)$$
假设总共有M哥误分类点，所有误分类点到线性分隔超平面的总距离为：
$$-\frac{1}{||w||}\sum _{x_i\in M}{y}_i(w\cdot x_0+b)$$
在忽略2-范数$\frac{1}{||w||}$的情况下，感知机的损失函数可表示为：
$$L(w,b)=-\sum _{x_i\in M}{y}_i(w\cdot x_i+b)$$
其中$M$是该分类点的集合。针对上述损失函数可使用随机梯度下降进行优化求解。分别计算损失函数$L(w,b)$关于$w$和$b$的梯度：
$$\frac{\partial L(w,b)}{\partial w}=-\sum _{x_i\in M}{y}_i{x}_i$$

$$\frac{\partial L(w,b)}{\partial b}=-\sum _{x_i\in M}{y}_i$$

根据上式更新权重系数：
$$w=w+\lambda y_i{x}_i$$

$$b=b+\lambda y_i$$

其中$\lambda$为学习步长（神经网络训练调参中的学习率）。

关于感知机模型，当一个实例被误分类时，即实例位于线性分隔超平面的错误一侧时，需要调整参数$w$和$b$的值，使得线性分隔超平面向该误分类点的一侧移动，以缩短该误分类点与分隔超平面的距离，直到线性分隔超平面越过该误分类点使其能够被正确分类。

9.1.2 基于Numpy的感知机实现

定义辅助函数

### 导入numpy模块
import numpy as np
# 定义sign符号函数
def sign(x, w, b):
    '''
    输入：
    x：输入实例
    w：权重系数
    b：偏置参数
    输出：
    符号函数值
    '''
    return np.dot(x,w)+b

# 定义参数初始化函数
def initialize_parameters(dim):
    '''
    输入：
    dim：输入数据维度
    输出：
    w：初始化后的权重系数
    b：初始化后的偏置参数
    '''
    w = np.zeros(dim, dtype=np.float32)
    b = 0.0
    return w, b

定义感知机训练过程

### 定义感知机训练函数
def train(X_train, y_train, learning_rate):
    '''
    输入：
    X_train：训练输入
    y_train：训练标签
    learing_rate：学习率
    输出：
    params：训练得到的参数
    '''
    # 参数初始化
    w, b = initialize_parameters(X_train.shape[1])
    # 初始化误分类状态
    is_wrong = False
    # 当存在误分类点时
    while not is_wrong:
        # 初始化误分类点计数
        wrong_count = 0
        # 遍历训练数据
        for i in range(len(X_train)):
            X = X_train[i]
            y = y_train[i]
            # 如果存在误分类点
            if y * sign(X, w, b) <= 0:
                # 更新参数
                w = w + learning_rate*np.dot(y, X)
                b = b + learning_rate*y
                # 误分类点+1
                wrong_count += 1
        # 直到没有误分类点        
        if wrong_count == 0:
            is_wrong = True
            print('There is no missclassification!')
        
        # 保存更新后的参数
        params = {
            'w': w,
            'b': b
        }
    return params

测试数据准备

# 导入pandas模块
import pandas as pd
# 导入iris数据集
from sklearn.datasets import load_iris
iris = load_iris()
# 转化为pandas数据框
df = pd.DataFrame(iris.data, columns=iris.feature_names)
# 数据标签
df['label'] = iris.target
# 变量重命名
df.columns = ['sepal length', 'sepal width', 'petal length', 'petal width', 'label']
# 取前100行数据
data = np.array(df.iloc[:100,[0, 1, -1]])
# 定义训练输入和输出
X,y = data[:,:-1],data[:,-1]
y = np.array([1 if i == 1 else -1 for i in y])
# 输出训练集大小
print(X.shape, y.shape)

感知机训练

params = perceptron_train(X, y, 0.01)
print(params)

绘制感知机的线性分隔超平面

# 导入matplotlib绘图库
import matplotlib as plt
# 输入实例取值
x_points = np.linspace(4, 7, 10)
# 线性分隔超平面
y_hat = -(params['w'][0]*x_points + params['b'])/params['w'][1]
# 绘制线性分隔超平面
plt.plot(x_points, y_hat)
# 绘制二分类散点图
plt.plot(data[:50, 0], data[:50, 1], color='red', label='0')
plt.plot(data[50:100, 0], data[50:100, 1], color='green', label='1')
plt.xlabel('sepal length')
plt.ylabel('sepal width')
plt.legend()
plt.show()

9.2 从单层到多层

9.2.1 神经网络与反向传播

单层感知机仅包含两层神经元，即输入神经元和输出神经元，容易实现逻辑与、逻辑或、逻辑非等线性可分情形，但异或问题单层感知机难以处理。

线性不可分

对于输入训练数据，不存在一个线性分隔超平面将其线性分类

对于线性不可分的情况，在感知机基础上一般两个处理方向SVM和神经网络模型（多层感知机）。

反向传播

反向传播（back propagation，BP）算法又称误差逆传播

假设输入层为$x$，有$m$个训练样本，输入层与隐藏层之间的权重和偏置分别为$w_1$和$b_1$，线性加权结果为$Z_1=w_{1}x+b_1$,采用Sigmoid激活函数，激活输出为$a_1=\sigma (Z_1)$，隐藏层到输出层权重和偏置为$w_2$和$b_2$，线性加权的计算结果为$Z_2=w_{2}x+b_2$，激活输出为$a_2=\sigma (Z_2)$。所以两层神经网络前向计算过程为$x\to Z_1\to a_1\to Z_2\to a_2$ 。

反向传播是将前向计算过程反过来，但必须是梯度计算的方向反过来，假设这里采用如下交叉熵损失函数：
$$L(y,a)=-(yloga+(1-y)log(1-a))$$
反向传播是基于梯度下降策略的，主要是从目标参数的负梯度方向更新参数，因此基于损失函数对前向计算过程中各个变量进行梯度计算是关键。

将前向计算过程反过来，基于损失函数的梯度计算顺序为：$d{a}_2\to d{Z}_2\to d{w}_2\to d{b}_2\to d{a}_1\to d{Z}_1\to d{w}_1\to d{b}_1$

9.2.2 基于Numpy的神经网络搭建

神经网络代码编写思路

Numpy：

• 基本思路
  • 网络结构
  • 参数初始化
  • 前向传播
  • 损失计算
  • 反向传播
  • 权重更新
  • 模块整合
• 数据测试
  • 数据准备
  • 训练与预测
  • 可视化

sklearn：

• sklearn.neural_network.MLPClassifier

定义网络结构

## 定义网络结构
def layer_sizes(X, Y):
    '''
    输入：
    X：训练输入
    Y：训练输出
    输出：
    n_x：输入层大小
    n_h：隐藏层大小
    n_y：输出层大小
    '''
    # 输入层大小
    n_x = X.shape[0]
     # 隐藏层大小
    n_h = 4
    # 输出层大小
    n_y = Y.shape[0] 
    return (n_x, n_h, n_y)

定义模型参数初始化函数

### 定义模型参数初始化函数
def initialize_parameters(n_x, n_h, n_y):
    '''
    输入：
    n_x：输入层神经元数
    n_h：隐藏层神经元数
    n_y：输出层神经元数
    输出：
    parameters：初始化后的模型参数
    '''
    # 权重系数随机初始化
    # 偏置参数以零为初始化值
    W1 = np.random.randn(n_h, n_x)*0.01
    b1 = np.zeros((n_h, 1))
    W2 = np.random.randn(n_y, n_h)*0.01
    b2 = np.zeros((n_y, 1)) 

    # 封装为字典
    parameters = {"W1": W1, 
                  "b1": b1,                 
                  "W2": W2,                  
                  "b2": b2}   
                   
    return parameters

定义前向传播过程

以$tanh$函数为隐藏层激活函数，$Sigmoid$函数为输出层激活函数，前向传播计算过程由以下四个公式定义：
$$z^{[1](i)}=W^{[1]}{x}^{(i)}+b^{[1][i]}$$

$$a^{[1][i]}=tanh(z^{[1](i)})$$

$$z^{[2](i)}=W^{[2]}{x}^{(i)}+b^{[2][i]}$$

$${\hat{y}}^{(i)}=a^{[2](i)}=\sigma z^{[2](i)}$$

def sigmoid(x):
    s = 1/(1+np.exp(-x))
    return s
# 定义前向传播过程
def forward_propagation(X, parameters):
    '''
    输入：
    X：训练输入
    parameters：初始化的模型参数
    输出：
    A2：模型输出
    cache：前向传播过程计算的中间值缓存
    '''
    # 获取各参数初始值
    W1 = parameters['W1']
    b1 = parameters['b1']
    W2 = parameters['W2']
    b2 = parameters['b2']    
    # 执行前向计算
    Z1 = np.dot(W1, X) + b1
    A1 = np.tanh(Z1)
    Z2 = np.dot(W2, A1) + b2
    A2 = sigmoid(Z2)    
    assert(A2.shape == (1, X.shape[1]))

    # 将中间结果封装为字典
    cache = {"Z1": Z1,                   
             "A1": A1,                   
             "Z2": Z2,                  
             "A2": A2}    

    return A2, cache

计算当前损失

前向计算输出结果后，将其与真实标签作比较，基于损失函数给出当前迭代的损失。基于交叉熵的损失定义：
$$L=-\frac{1}{m}\sum _{i=0}^m(y^{(i)}log{a}^{[2](i)}+(1-y^{(i)}log(1-a^{[2](i)}))$$

# 定义损失函数
def compute_cost(A2, Y):
    '''
    输入：
    A2：前向计算输出
    Y：训练标签
    输出：
    cost：当前损失
    '''
    # 训练样本量
    m = Y.shape[1] 
    # 计算交叉熵损失
    logprobs = np.multiply(np.log(A2),Y) + np.multiply(np.log(1-A2), 1-Y)
    cost = -1/m * np.sum(logprobs)
    # 维度压缩
    cost = np.squeeze(cost)     
    return cost

执行反向传播

前向传播和损失计算完后，神经网络最关键、最核心的部分是执行反向传播。损失函数关于各参数的梯度下降计算公式如下：
$$d{z}^{[2]}=a^{[2]}-y$$

$$d{W}^{[2]}=d{z}^{[2]}{a}^{[1{]}^T}$$

$$d{b}^{[2]}=d{z}^{[2]}$$

$$d{z}^{[1]}=W^{[2]T}d{z}^{[2]}\cdot g^{[1]}(z^{[1]})$$

$$d{W}^{[1]}=d{z}^{[1]}{x}^T$$

$$d{b}^{[1]}=d{z}^{[1]}$$

### 定义反向传播过程
def backward_propagation(parameters, cache, X, Y):
    '''
    输入：
    parameters：神经网络参数字典
    cache：神经网络前向计算中间缓存字典
    X：训练输入
    Y：训练输出
    输出：
    grads：权重梯度字典
    '''
    # 样本量
    m = X.shape[1]    
    # 获取W1和W2
    W1 = parameters['W1']
    W2 = parameters['W2']    
    # 获取A1和A2
    A1 = cache['A1']
    A2 = cache['A2']    
    # 执行反向传播
    dZ2 = A2-Y
    dW2 = 1/m * np.dot(dZ2, A1.T)
    db2 = 1/m * np.sum(dZ2, axis=1, keepdims=True)
    dZ1 = np.dot(W2.T, dZ2)*(1-np.power(A1, 2))
    dW1 = 1/m * np.dot(dZ1, X.T)
    db1 = 1/m * np.sum(dZ1, axis=1, keepdims=True)

    # 将权重梯度封装为字典
    grads = {"dW1": dW1,
             "db1": db1,                      
             "dW2": dW2,             
             "db2": db2}   
    return grads

更新权重

有了梯度计算结果后，根据权重更新公式更新权重和偏置参数，具体计算公式如下，其中$\eta$为学习率。
$$w=w-\eta dw$$
反向传播完成后，基于权重梯度更新权重，按照计算公式，对权重按照负梯度方向不断迭代（梯度下降法），即可一步一步达到最优值。

# 定义权重更新过程
def update_parameters(parameters, grads, learning_rate=1.2):
    '''
    输入：
    parameters：神经网络参数字典
    grads：权重梯度字典
    learning_rate：学习率
    输出：
    parameters：更新后的权重字典
    '''
    # 获取参数
    W1 = parameters['W1']
    b1 = parameters['b1']
    W2 = parameters['W2']
    b2 = parameters['b2']    
    # 获取梯度
    dW1 = grads['dW1']
    db1 = grads['db1']
    dW2 = grads['dW2']
    db2 = grads['db2']    
    # 参数更新
    W1 -= dW1 * learning_rate
    b1 -= db1 * learning_rate
    W2 -= dW2 * learning_rate
    b2 -= db2 * learning_rate

   # 将更新后的权重封装为字典
    parameters = {"W1": W1, 
                  "b1": b1,            
                  "W2": W2,   
                  "b2": b2}    
    return parameters

模块整合

### 神经网络模型封装
def nn_model(X, Y, n_h, num_iterations=10000, print_cost=False):
    '''
    输入：
    X：训练输入
    Y：训练输出
    n_h：隐藏层结点数
    num_iterations：迭代次数
    print_cost：训练过程是否打印损失
    输出：
    parameters：神经网络训练优化后的权重系数
    '''
    # 设置随机数种子
    np.random.seed(3)
    # 输入和输出结点数
    n_x = layer_sizes(X, Y)[0]
    n_y = layer_sizes(X, Y)[2]    
    # 初始化模型参数
    parameters = initialize_parameters(n_x, n_h, n_y)
    W1 = parameters['W1']
    b1 = parameters['b1']
    W2 = parameters['W2']
    b2 = parameters['b2']    
    # 梯度下降和参数更新循环
    for i in range(0, num_iterations):        
    # 前向传播计算
        A2, cache = forward_propagation(X, parameters)        
        # 计算当前损失
        cost = compute_cost(A2, Y, parameters)        
        # 反向传播
        grads = backward_propagation(parameters, cache, X, Y)        
        # 参数更新
        parameters = update_parameters(parameters, grads, learning_rate=1.2)        
        # 打印损失
        if print_cost and i % 1000 == 0:            
            print ("Cost after iteration %i: %f" %(i, cost))    
            
    return parameters

[以上是完整的两层全连接神经网络代码示例]

生成模拟数据集

# 生成非线性可分数据集
def create_dataset():
    '''
    输入：
    无
    输出：
    X：模拟数据集输入
    Y：模拟数据集输出
    '''
    # 设置随机数种子
    np.random.seed(1)
    # 数据量
    m = 400
    # 每个标签的实例数
    N = int(m/2)
    # 数据维度
    D = 2 
    # 数据矩阵
    X = np.zeros((m,D)) 
    # 标签维度
    Y = np.zeros((m,1), dtype='uint8') 
    a = 4 
    # 遍历生成数据
    for j in range(2):
        ix = range(N*j,N*(j+1))
        t = np.linspace(j*3.12,(j+1)*3.12,N) + np.random.randn(N)*0.2 # theta
        r = a*np.sin(4*t) + np.random.randn(N)*0.2 # radius
        X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
        Y[ix] = j
        
    X = X.T
    Y = Y.T

    return X, Y

模型训练

### 模型训练
parameters = nn_model(X, Y, n_h = 4, num_iterations=10000, print_cost=True)

9.2.3 基于sklearn神经网络搭建

# 导入sklearn神经网络模块
from sklearn.neural_network import MLPClassifier
# 创建神经网络分类器
clf = MLPClassifier(solver='sgd', alpha=1e-5, hidden_layer_sizes=(4), random_state=1)
clf.fit(X.T, Y.T)

9.3 小结

神经网络的核心是反向传播算法，一个典型的神经网络一般有前向计算到反向传播的算法流程，链式求导法则是反向传播的核心操作。

第 10 章 支持向量机

支持向量机是一种二分类模型，通过不同的间隔最大化策略，支持向量机模型可以分为线性可分支持向量机、近似线性可分支持向量机和线性不可分支持向量机。

10.1 由感知机到支持向量机

感知机是一种通过寻找一个线性分隔超平面将正负实例分隔开来的分类模型。

支持向量机是从感知机的无穷多个解中选取一个到两边实例最大间隔的线性分隔超平面。

• 当训练数据线性可分，SVM通过求硬间隔最大化来求最优线性分隔超平面
• 当训练数据近似线性可分，SVM通过求软间隔最大化来求最优线性分隔超平面
• 当训练数据线性不可分，SVM通过使用核函数和软间隔最大化，将线性不可分问题转换成线性可分问题，从而实现分类

10.2 线性可分支持向量机

10.2.1 线性可分支持向量机的原理推导

当训练数据线性线性可分时，能够通过硬间隔（hard margin）最大化求解对应的凸二次规划问题得到ui有线性分隔超平面$w^{\ast }\cdot x+b^{\ast }=0$，以及相应的分类决策函数$f(x)=sign(w^{\ast }\cdot x+b^{\ast })$，该情况称为线性可分支持向量机。

对于支持向量机而言，一个实例点到线性分隔超平面的距离可以表示为分类预测的可靠度，当分类的线性分隔超平面$w\cdot x+b=0$确定时，$|w\cdot x+b|$可以表示点$x$与该超平面的距离。同时用$w\cdot x+b$的符号与分类标记$y$符号的一致性来判定分类是否正确。因此可以用量$y(w\cdot x+b)$来表示分类的正确性以及确信度，这就是函数间隔的概念。

对于给定训练样本和线性分隔超平面$w\cdot x+b=0$，线性分隔超平面关于任意样本点$(x_i,y_i)$的函数间隔表示为：
$$\widehat{d_i}=y_i(w\cdot x_i+b)$$
那么该训练集与线性分隔超平面的间隔可以由该超平面与所有样本点的最小函数间隔决定，即
$$\hat{d}=\underset{i=1,...,N}{\min }\widehat{d_i}$$
为了使间隔不受线性分隔超平面参数$w$和$b$的变化影响，对$w$加一个规范化约束$||w||$使间隔确定，这时函数间隔变为几何间隔，即线性分隔超平面关于任意样本点$(x_i,y_i)$的几何间隔表示为：
$$d_i=y_i(\frac{w}{||w||}\cdot x_i+\frac{b}{||w||})$$
硬间隔最大化理解为以足够高的可靠度对训练数据进行分类，以此求得的线性分隔超平面不仅将正负实例点分开，而且对于最难分的实例点也能以最够高的可靠度将其分类。

将硬间隔最大化表示为条件约束最优化问题：
$$\begin{gathered} \underset{w,b}{\max }d \\ s.t.y_i(\frac{w}{||w||}\cdot x_i+\frac{b}{||w||})\ge d,i=1,2,...,N \end{gathered}$$

即我们希望最大化超平面$(w,b)$关于训练数据集的集合间隔$d$，约束条件表示超平面$(w,b)$关于每个训练样本点的几何间隔至少是$d$。

考虑函数间隔与几何间隔的关系$d=\frac{\hat{d}}{||w||}$，将问题改写为：
$$\begin{gathered} \underset{w,b}{\max }\frac{\hat{d}}{||w||} \\ s.t.y_i(w\cdot x_i+b)\ge \hat{d},i=1,2,...,N \end{gathered}$$

函数间隔的取值并不影响最优化问题的解，即使将$w$和$b$按比例改为$\lambda w$和$\lambda b$，函数间隔成为$\lambda \hat{d}$，函数间隔的改变对不等式约束没有影响，假设取函数间隔为1，将$\hat{d}=1$代入最优化问题，由于最大化$\frac{1}{||w||}$等价于最小化$\frac{1}{2}||w|{|}^2$，于是线性可分支持向量机学习的最优化问题为：
$$\begin{gathered} \underset{w,b}{\min }\frac{1}{2}||w|{|}^2 \\ s.t{y}_i(w\cdot x_i+b)-1\ge 0,i=1,2,...,N \end{gathered}$$

至此，硬间隔最大化问题转换为典型的凸二次规划问题。

为了求解线性可分支持向量机的最优化问题，将上最优化问题应用拉格朗日对偶性，通过求解对偶问题得到问题的最优解，这便是线性可分支持向量机的对偶算法。

首先构建拉格朗日函数，如下：
$$L(w,b,a)=\frac{1}{2}||w|{|}^2-\sum _{i=1}^N{a}_i{y}_i(w\cdot x_i+b)+\sum _{i=1}^N{a}_i，其中a=(a_1,a_2,...,a_N{)}^T为拉格朗日乘子向量。$$
根据拉格朗日对偶性，原始问题的对偶问题是极大极小问题：
$$\underset{a}{\max }\underset{w,b}{\min }L(w,b,a)$$
第一步，求$\underset{w,b}{\min }L(w,b,a)$

将拉格朗日函数$L(w,b,a)$分别对$w$和$b$求偏导并令其等于$0$
$$\begin{gathered} \frac{\partial L}{\partial w}=w-\sum _{i=1}^N{a}_i{y}_i{x}_i=0 \\ \frac{\partial L}{\partial b}=\sum _{i=1}^N{a}_i{y}_i=0 \end{gathered}$$

得，
$$\begin{gathered} w=\sum _{i=1}^N{a}_i{y}_i{x}_i \\ \sum _{i=1}^N{a}_i{y}_i=0 \end{gathered}$$

将$w={\sum }_{i=1}^N{a}_i{y}_i{x}_i$代入拉格朗日函数，有：
$$\begin{gathered} \underset{w,b}{\min }L(w,b,a)=\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{a}_i{a}_j{y}_i{y}_j(x_i\cdot x_j)-\sum _{i=1}^N{a}_i{y}_i((\sum _{j=1}^N{a}_j{y}_j{x}_j)\cdot x_i+b)+\sum _{i=1}^N{a}_i \\ =-\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{a}_i{a}_j{y}_i{y}_j(x_i\cdot x_j)+\sum _{i=1}^N{a}_i \end{gathered}$$
第二步，求$\underset{a}{\max }\underset{w,b}{\min }L(w,b,a)$，即求$\underset{w,b}{\min }L(w,b,a)$对$a$的极大

假设$a^{\ast }=(a_1^{\ast },a_2^{\ast },...,a_i^{\ast }{)}^T$是对偶最优化问题的解，根据拉格朗日对偶性相关理论，拉格朗日函数满足KKT条件：

可以解得，

相应的线性可分支持向量机的线性分隔超平面表达为：
$$\sum _{i=1}^N{a}_i^{\ast }{y}_i(x\cdot x_i)+b^{\ast }=0$$

10.2.2 线性可分支持向量机的算法实现

cvxopt 用法

经典的二次规划问题可表示为如下形式：

假设要求解如下二次规划问题：

将目标函数和约束条件写成矩阵形式：

from cvxopt import matrix, solvers

# 定义二次规划参数
P = matrix([[1.0,0.0],[0.0,0.0]])
q = matrix([3.0,4.0])
G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])
h = matrix([0.0,0.0,-15.0,100.0,80.0])

# 构建求解
sol = solvers.qp(P, q, G, h)

生成模拟二分类数据集

# 导入相关库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split

# 导入sklearn模拟二分类数据生成模块
from sklearn.datasets.samples_generator import make_blobs
# 生成模拟二分类数据集
X, y =  make_blobs(n_samples=150, n_features=2, centers=2, cluster_std=1.2, random_state=40)
# 将标签转换为1/-1
y_ = y.copy()
y_[y_==0] = -1
y_ = y_.astype(float)
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y_, test_size=0.3, random_state=43)
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)
# 设置颜色参数
colors = {0:'r', 1:'g'}
# 绘制二分类数据集的散点图
plt.scatter(X[:,0], X[:,1], marker='o', c=pd.Series(y).map(colors))
plt.show();

定义线性可分支持向量机类

### 实现线性可分支持向量机
### 硬间隔最大化策略
class Hard_Margin_SVM:
    ### 线性可分支持向量机拟合方法
    def fit(self, X, y):
        # 训练样本数和特征数
        m, n = X.shape

        # 初始化二次规划相关变量：P/q/G/h
        self.P = matrix(np.identity(n + 1, dtype=np.float))
        self.q = matrix(np.zeros((n + 1,), dtype=np.float))
        self.G = matrix(np.zeros((m, n + 1), dtype=np.float))
        self.h = -matrix(np.ones((m,), dtype=np.float))

        # 将数据转为变量
        self.P[0, 0] = 0
        for i in range(m):
            self.G[i, 0] = -y[i]
            self.G[i, 1:] = -X[i, :] * y[i]
        
        # 构建二次规划求解
        sol = solvers.qp(self.P, self.q, self.G, self.h)

        # 对权重和偏置寻优
        self.w = np.zeros(n,) 
        self.b = sol['x'][0] 
        for i in range(1, n + 1):
            self.w[i - 1] = sol['x'][i]
        return self.w, self.b

    ### 定义模型预测函数
    def predict(self, X):
        return np.sign(np.dot(self.w, X.T) + self.b)

基于cvxopt的线性可分支持向量机训练数据

# 创建线性可分支持向量机模型实例
hard_margin_svm = Hard_Margin_SVM()
# 执行训练
hard_margin_svm.fit(X_train, y_train)
# 模型预测
y_pred = hard_margin_svm.predict(X_test)
from sklearn.metrics import accuracy_score
# 计算测试集准确率
print(accuracy_score(y_test, y_pred))

10.2.3 基于sklearn的线性可分支持向量机的算法实现

# 导入sklearn线性SVM分类模块
from sklearn.svm import LinearSVC
# 创建模型实例
clf = LinearSVC(random_state=0, tol=1e-5)
# 训练
clf.fit(X_train, y_train)
# 预测
y_pred = clf.predict(X_test)
# 计算测试集准确率
print(accuracy_score(y_test, y_pred))

10.3 近似线性可分支持向量机

10.3.1 近似线性可分支持向量机的原理推导

近似线性可分是指训练集中大部分实例点是线性可分的，只是一些特殊实例点的存在使得这种数据集不适用于直接使用线性可分支持向量机进行处理，但也没到完全线性不可分的程度。

近似线性可分问题需要采用软间隔最大化的策略处理。线性不可分意味着少数样本点不能满足函数间隔大于等于1的约束条件，近似线性可分支持向量机的解决方案是对每个这样的特殊实例引入一个松弛变量${\xi }_i\ge 0$，使得函数间隔加上松弛变量大于等于1，因此约束条件变为$y_i(w\cdot x_i+b)\ge 1-{\xi }_i$，每个松弛变量${\xi }_i$，支付一个代价${\xi }_i$。目标函数由原来的$\frac{1}{2}||w|{|}^2$变成$\frac{1}{2}||w|{|}^2+C{\sum }_{i=1}^N{\xi }_i$。

这里$C$为惩罚系数，C值大对误分类的惩罚增大，C值小对误分类的惩罚减小。最小化目标函数指

• 使$\frac{1}{2}||w|{|}^2$尽可能小，即间隔尽量大
• 误分类点的个数尽量小

根线性可分支持向量机一样，它可以变成如下的凸二次规划问题：

将其转换成对偶问题求最优解，上式对偶问题是：
$$\begin{gathered} \underset{a}{\min }\frac{1}{2}\sum _{i=1}^N\sum _{j=1}^N{a}_i{a}_j{y}_i{y}_j-\sum _{i=1}^N{a}_i \\ s.t.\sum _{i=1}^N{a}_i{y}_i=0 \\ 0\le a_i\le C,,i=1,2,...,N \end{gathered}$$
因此原始最优化问题的拉格朗日函数为：

，对偶问题为拉格朗日函数的极大极小问题。基于该拉格朗日函数对$w$、$b$和$\xi {}_i$求偏导：

将上述三个解代入拉格朗日函数中，

再对$\underset{w,b,\xi }{\min }L(w,b,\xi ,\alpha ,\mu )$求$\alpha$的极大，得到线性支持向量机的对偶问题：
$$\begin{gathered} \underset{\alpha }{\max }-\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.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ C-{\alpha }_i-{\mu }_i=0 \\ {\alpha }_i\ge 0 \\ {\mu }_i\ge 0,i=1,2,...,N \end{gathered}$$
消去变量${\mu }_i$后简化约束条件为：
$$0\le {\alpha }_i\le C$$
通过求解对偶问题得到原始问题的解，进而确定分离超平面和决策函数。假设${\alpha }^{\ast }=({\alpha }_1^{\ast },{\alpha }_2^{\ast },...,{\alpha }_N^{\ast }{)}^T$是对偶最优化问题的解，根据拉格朗日对偶性理论，原始问题是凸二次规划问题，解满足KKT条件，有：

可解得，

因此，分离超平面写成，
$$\sum _{i=1}^N{\alpha }_i^{\ast }{y}_i(x\cdot x_i)+b^{\ast }=0$$
分类决策函数写成，
$$f(x)=sign(\sum _{i=1}^N{\alpha }_i^{\ast }{y}_i(x\cdot x_i)+b^{\ast })$$

10.3.2 近似线性可分支持向量机的算法实现

与线性可分的数据集不同：该数据集有部分数据重叠，使分类任务近似线性可分。

生成近似线性可分模拟数据集

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 给定二维正态分布均值矩阵
mean1, mean2 = np.array([0, 2]), np.array([2, 0])
# 给定二维正态分布协方差矩阵
covar = np.array([[1.5, 1.0], [1.0, 1.5]])
# 生成二维正态分布样本X1
X1 = np.random.multivariate_normal(mean1, covar, 100)
# 生成X1的标签 1
y1 = np.ones(X1.shape[0])
# 生成二维正态分布样本X2
X2 = np.random.multivariate_normal(mean2, covar, 100)
# 生成X2的标签 -1
y2 = -1 * np.ones(X2.shape[0])
# 设置训练集和测试集
X_train = np.vstack((X1[:80], X2[:80]))
y_train = np.hstack((y1[:80], y2[:80]))
X_test = np.vstack((X1[80:], X2[80:]))
y_test = np.hstack((y1[80:], y2[80:]))
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)
(160, 2) (160,) (40, 2) (40,)

# 设置颜色参数
colors = {1:'r', -1:'g'}
# 绘制二分类数据集的散点图
plt.scatter(X_train[:,0], X_train[:,1], marker='o', c=pd.Series(y_train).map(colors))
plt.show();

定义近似线性可分支持向量机类

from cvxopt import matrix, solvers

### 定义一个线性核函数
def linear_kernel(x1, x2):
    '''
    输入:
    x1: 向量1
    x2: 向量2
    输出:
    np.dot(x1, x2): 两个向量的点乘
    '''
    return np.dot(x1, x2)
from cvxopt import matrix, solvers

### 定义近似线性可分支持向量机
### 软间隔最大化策略
class Soft_Margin_SVM:
    ### 定义基本参数
    def __init__(self, kernel=linear_kernel, C=None):
        # 软间隔svm核函数，默认为线性核函数
        self.kernel = kernel
        # 惩罚参数
        self.C = C
        if self.C is not None: 
            self.C = float(self.C)
    
    ### 定义线性支持向量机拟合方法
    def fit(self, X, y):
        # 训练样本数和特征数
        m, n = X.shape
        
        # 基于线性核计算Gram矩阵
        K = self._gram_matrix(X)
                
        # 初始化二次规划相关变量：P/q/G/h
        P = matrix(np.outer(y,y) * K)
        q = matrix(np.ones(m) * -1)
        A = matrix(y, (1, m))
        b = matrix(0.0)
        
        # 未设置惩罚参数时的G和h矩阵
        if self.C is None:
            G = matrix(np.diag(np.ones(m) * -1))
            h = matrix(np.zeros(m))
        # 设置惩罚参数时的G和h矩阵
        else:
            tmp1 = np.diag(np.ones(m) * -1)
            tmp2 = np.identity(m)
            G = matrix(np.vstack((tmp1, tmp2)))
            tmp1 = np.zeros(m)
            tmp2 = np.ones(m) * self.C
            h = matrix(np.hstack((tmp1, tmp2)))

        # 构建二次规划求解
        sol = solvers.qp(P, q, G, h, A, b)
        # 拉格朗日乘子
        a = np.ravel(sol['x'])

        # 寻找支持向量
        spv = a > 1e-5
        ix = np.arange(len(a))[spv]
        self.a = a[spv]
        self.spv = X[spv]
        self.spv_y = y[spv]
        print('{0} support vectors out of {1} points'.format(len(self.a), m))

        # 截距向量
        self.b = 0
        for i in range(len(self.a)):
            self.b += self.spv_y[i]
            self.b -= np.sum(self.a * self.spv_y * K[ix[i], spv])
        self.b /= len(self.a)

        # 权重向量
        self.w = np.zeros(n,)
        for i in range(len(self.a)):
            self.w += self.a[i] * self.spv_y[i] * self.spv[i]

    ### 定义Gram矩阵计算函数
    def _gram_matrix(self, X):
        m, n = X.shape
        K = np.zeros((m, m))
        # 遍历计算Gram矩阵
        for i in range(m):
            for j in range(m):
                K[i,j] = self.kernel(X[i], X[j])
        return K
    
    ### 定义映射函数
    def project(self, X):
        if self.w is not None:
            return np.dot(X, self.w) + self.b
    
    ### 定义模型预测函数
    def predict(self, X):
        return np.sign(np.dot(self.w, X.T) + self.b)

基于cvxopt的近似线性可分支持向量机训练数据

# 导入准确率评估模块
from sklearn.metrics import accuracy_score
# 构建线性可分支持向量机实例，设置惩罚参数为 0.1
soft_margin_svm = Soft_Margin_SVM(C=0.1)
# 模型拟合
soft_margin_svm.fit(X_train, y_train)
# 模型预测
y_pred = soft_margin_svm.predict(X_test)
# 计算测试集准确率
print('Accuracy of soft margin svm based on cvxopt: ', 
      accuracy_score(y_test, y_pred))

10.3.3 基于sklearn的近似线性可分支持向量机的实现

from sklearn import svm
# 创建svm模型实例
clf = svm.SVC(kernel='linear')
# 模型拟合
clf.fit(X_train, y_train)
# 模型预测
y_pred = clf.predict(X_test)
# 计算测试集准确率
print('Accuracy of soft margin svm based on sklearn: ', 
      accuracy_score(y_test, y_pred))

10.4 线性不可分支持向量机

10.4.1 线性不可分与核技巧

非线性可分问题，就是对于给定数据集，如果能用一个超曲面将正负实例正确分开，则这个问题为非线性可分问题。非线性问题的一个关键在于将原始数据空间转换到一个新的数据空间，在原始空间中的非线性可分问题到新空间就是是线性可分问题。

核技巧

线性可分方法来解决非线性可分问题可分为两步

• 首先用一个变换将原始空间的数据映射到新空间
• 再在新空间中用线性分类学习方法训练分类模型。

核函数

假设输入空间为$\chi$，特征空间为KaTeX parse error: Can't use function '\H' in math mode at position 1: \̲H̲，若存在一个从$\chi$到KaTeX parse error: Can't use function '\H' in math mode at position 1: \̲H̲的映射$\varphi (x)$，使得所有的$x，z\in \chi$，函数$K(x,z)=\varphi (x)\cdot \varphi (z)$，则$K(x,z)$是核函数，$\varphi ()x$和$\varphi (z)$是映射函数。

在线性支持向量机的对偶问题中，无论目标函数还是决策函数（分离超平面）只涉及输入实例与实例之间的内积，因此在对偶问题的目标函数$\underset{\alpha }{\min }\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$中的内积$x_i\cdot x_j$用核函数$K(x,z)=\varphi (x)\cdot \varphi (z)$代替。此时对偶问题的目标函数成为
$$\underset{\alpha }{\min }\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$$

同样分类函数中的内积也可以用核函数代替，因此分类决策函数式成为：
$$\begin{gathered} f(x)=sign(\sum _{i=1}^N{\alpha }_i^{\ast }{y}_i\varphi (x_i)\cdot \varphi (x)+b^{\ast }) \\ =sign(\sum _{i=1}^N{\alpha }_i^{\ast }{y}_{i}K(x_i,x)+b^{\ast }) \end{gathered}$$
这等价于经过映射函数$\varphi$将原来的输入空间变换到新的特征空间，将输入空间中的内积$x_1\cdot x_2$变换为特征空间中的内积$\varphi (x_i)\cdot \varphi (x_j)$，在新的特征空间里从训练样本中学习线性支持向量机。当映射函数是非线性函数时，学习到的含有核函数的支持向量机是非线性模型。也就是说在核函数$K(x,z)$给定的条件下，可以利用解线性分类问题的方法求解非线性分类问题的支持向量机。

线性不可分的支持向量机学习算法

给定训练数据 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x N , y N ) } \{(x_1,y_1),(x_2,y_2),...,(x_N,y_N)\} {(x1​,y1​),(x2​,y2​),...,(xN​,yN​)}，其中$x_i\in R^n$， y i ∈ Y = { + 1 , − 1 } y_i∈Y=\{+1,-1\} yi​∈Y={+1,−1}，$i=1,2,...,N$

• 选取合适的核函数$K(x,z)$和参数$C$，构造并求解最优化问题：

$$\begin{gathered} \underset{\alpha }{\min }\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.\sum _{i=1}^N{\alpha }_i{y}_i=0 \\ 0\le {\alpha }_i\le C,i=1,2,...,N \end{gathered}$$

可求得最优解${\alpha }^{\ast }=({\alpha }_1^{\ast },{\alpha }_2^{\ast },...,{\alpha }_N^{\ast }{)}^T$。

• 选择${\alpha }^{\ast }$的一个正分量$0<{\alpha }_j^{\ast }<C$，计算：

$$b^{\ast }=y_i-\sum _{i=1}^N{\alpha }_i^{\ast }{y}_{i}K(x_i,x_j)$$

• 构造决策函数

$$f(x)=sign(\sum _{i=1}^N{\alpha }_i^{\ast }{y}_{i}K(x_i,x)+b^{\ast })$$

当$K(x,z)$是正定核时，该最优化问题是凸二次规划问题，一般使用SMO（序列最小最优化）算法来解。

[序列最小最优化算法 参考 李航. 《统计学习方法》第 7 章 支持向量机]

10.4.2 线性不可分支持向量机算法实现

生成完全线性不可分的模拟数据集

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# 给定二维正态分布均值矩阵
mean1, mean2 = np.array([-1, 2]), np.array([1, -1])
mean3, mean4 = np.array([4, -4]), np.array([-4, 4])
# 给定二维正态分布协方差矩阵
covar = np.array([[1.0, 0.8], [0.8, 1.0]])
# 生成二维正态分布样本X1
X1 = np.random.multivariate_normal(mean1, covar, 50)
# 合并两个二维正态分布并令其为新的X1
X1 = np.vstack((X1, np.random.multivariate_normal(mean3, covar, 50)))
# 生成X1的标签 1
y1 = np.ones(X1.shape[0])
# 生成二维正态分布样本X2
X2 = np.random.multivariate_normal(mean2, covar, 50)
# 合并两个二维正态分布并令其为新的X2
X2 = np.vstack((X2, np.random.multivariate_normal(mean4, covar, 50)))
# 生成X2的标签 -1
y2 = -1 * np.ones(X2.shape[0])
# 划分训练集和测试集
X_train = np.vstack((X1[:80], X2[:80]))
y_train = np.hstack((y1[:80], y2[:80]))
X_test = np.vstack((X1[80:], X2[80:]))
y_test = np.hstack((y1[80:], y2[80:]))
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)

# 设置颜色参数
colors = {1:'r', -1:'g'}
# 绘制二分类数据集的散点图
plt.scatter(X_train[:,0], X_train[:,1], marker='o', c=pd.Series(y_train).map(colors))
plt.show();

定义线性不可分支持向量机

from cvxopt import matrix, solvers

### 定义高斯核函数
def gaussian_kernel(x1, x2, sigma=5.0):
    '''
    输入:
    x1: 向量1
    x2: 向量2
    输出:
    两个向量的高斯核
    '''
    return np.exp(-1 * np.linalg.norm(x1-x2)**2 / (2 * (sigma ** 2)))

### 定义线性不可分支持向量机
### 借助于高斯核函数转化为线性可分的情形
class Non_Linear_SVM:
    ### 定义基本参数
    def __init__(self, kernel=gaussian_kernel):
        # 非线性可分svm核函数，默认为高斯核函数
        self.kernel = kernel
    
    ### 定义非线性可分支持向量机拟合方法
    def fit(self, X, y):
        # 训练样本数和特征数
        m, n = X.shape
        
        # 基于线性核计算Gram矩阵
        K = self._gram_matrix(X)
                
        # 初始化二次规划相关变量：P/q/A/b/G/h
        P = matrix(np.outer(y,y) * K)
        q = matrix(np.ones(m) * -1)
        A = matrix(y, (1, m))
        b = matrix(0.0)
        G = matrix(np.diag(np.ones(m) * -1))
        h = matrix(np.zeros(m))

        # 构建二次规划求解
        sol = solvers.qp(P, q, G, h, A, b)
        # 拉格朗日乘子
        a = np.ravel(sol['x'])

        # 寻找支持向量
        spv = a > 1e-5
        ix = np.arange(len(a))[spv]
        self.a = a[spv]
        self.spv = X[spv]
        self.spv_y = y[spv]
        print('{0} support vectors out of {1} points'.format(len(self.a), m))

        # 截距向量
        self.b = 0
        for i in range(len(self.a)):
            self.b += self.spv_y[i]
            self.b -= np.sum(self.a * self.spv_y * K[ix[i], spv])
        self.b /= len(self.a)

        # 权重向量
        self.w = None

    ### 定义Gram矩阵计算函数
    def _gram_matrix(self, X):
        m, n = X.shape
        K = np.zeros((m, m))
        # 遍历计算Gram矩阵
        for i in range(m):
            for j in range(m):
                K[i,j] = self.kernel(X[i], X[j])
        return K
    
    ### 定义映射函数
    def project(self, X):
        y_pred = np.zeros(len(X))
        for i in range(X.shape[0]):
            s = 0
            for a, spv_y, spv in zip(self.a, self.spv_y, self.spv):
                s += a * spv_y * self.kernel(X[i], spv)
            y_pred[i] = s
        return y_pred + self.b
    
    ### 定义模型预测函数
    def predict(self, X):
        return np.sign(self.project(X))

基于cvxopt的线性不可分支持向量机训练数据

# 导入sklearn准确率评估函数
from sklearn.metrics import accuracy_score
# 创建非线性可分支持向量机模型实例
non_linear_svm = Non_Linear_SVM()
# 模型拟合
non_linear_svm.fit(X_train, y_train)
# 模型预测
y_pred = non_linear_svm.predict(X_test)
# 计算测试集准确率
print('Accuracy of soft margin svm based on cvxopt: ', 
      accuracy_score(y_test, y_pred))

绘制训练的线性不可分支持向量机的分隔超平面

### 绘制非线性可分支持向量机
def plot_classifier(X1_train, X2_train, clf):
    plt.plot(X1_train[:,0], X1_train[:,1], "ro")
    plt.plot(X2_train[:,0], X2_train[:,1], "go")
    plt.scatter(non_linear_svm.spv[:,0], non_linear_svm.spv[:,1], 
                s=100, c="", edgecolors="b", label="support vector")

    X1, X2 = np.meshgrid(np.linspace(-4,4,50), np.linspace(-4,4,50))
    X = np.array([[x1, x2] for x1, x2 in zip(np.ravel(X1), np.ravel(X2))])
    Z = non_linear_svm.project(X).reshape(X1.shape)
    plt.contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower')
    plt.contour(X1, X2, Z + 1, [0.0], colors='grey', linewidths=1, origin='lower')
    plt.contour(X1, X2, Z - 1, [0.0], colors='grey', linewidths=1, origin='lower')
    plt.legend()
    plt.show()
    
plot_classifier(X_train[y_train==1], X_train[y_train==-1], non_linear_svm)

10.4.3 基于sklearn的线性不可分支持向量机实现

from sklearn import svm
# 创建svm模型实例
clf = svm.SVC(kernel='rbf')
# 模型拟合
clf.fit(X_train, y_train)
# 模型预测
y_pred = clf.predict(X_test)
# 计算测试集准确率
print('Accuracy of soft margin svm based on sklearn: ', 
      accuracy_score(y_test, y_pred))

10.5 小结

• 支持向量机是感知机实现非线性的另一种方法

• 支持向量机通过定义非线性核函数的方式来实现非线性

• 线性可分支持向量机相应模型优化目标转换为硬间隔最大化问题

• 近似线性可分支持向量机相应模型的优化目标转换为软间隔最大化问题

• 线性不可分支持向量机相应模型的优化目标转换为应用核函数后的间隔最大化问题

• 本篇主要给出基本数学推导，介绍模型优化的原始问题和对偶问题以及相应的求解方式，在此基础上，基于cvxopt二次规划求解库，对三种模型进行代码实现