python学智能算法（十六）|机器学习支持向量机简单示例-EW帮帮网

【1】引言

前序学习了逻辑回归等算法，相关文章链接包括且不限于：

python学智能算法（十）|机器学习逻辑回归（Logistic回归）_逻辑回归算法python-CSDN博客

python学智能算法（十一）|机器学习逻辑回归深入（Logistic回归）_np.random.logistic()-CSDN博客

今天在此基础上更进一步，学习支持向量机，为实现较好地理解，先解读一个简单算例。

【2】代码解读

【2.1】模块引入

模块引入部分除了已经非常熟悉的numpy和matplotlib外，主要就是sklearn模块的数据和常规操作。

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

【2.2】数据集划分

数据计划分直接调用sklearn中成熟的鸢尾花数据datasets。

一共有三种花，每种花分别有50个样本，每个样本都具有4个特征：花萼长度、花萼宽度、花瓣长度、花瓣宽度。

# 加载并预处理数据
iris = datasets.load_iris()
# iris的前 100 个样本的 4 个特征（数值）（花萼长度、花萼宽度、花瓣长度、花瓣宽度）
X = iris.data[:100]  # 只取前两类（线性可分）
# iris.target共有150个样本，每个样本50个，前100个样本对应的标签分别是0和1
y = iris.target[:100]
# 如果y=0,输出-1，否则输出1
y = np.where(y == 0, -1, 1)  # 将标签转换为-1和1

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

这里的iris.data[:100]取了前100个样本的数据，刚好对应前两种花朵的样本数，所以取iris.target[:100]刚好一一对应。

【2.3】归一化处理

数据标准化处理的目的是：将特征缩放至均值为 0、标准差为 1 的标准正态分布

# 数据标准化
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

【2.4】支持向量机函数

为了更好地理解支持向量机的原理，这里详细定义了支持向量机的实现过程：

class SVM:
    def __init__(self, C=1.0, kernel='linear', gamma=0.1, tol=1e-3, max_iter=100):
        # 定义一批self量存储数据
        self.C = C
        self.kernel = kernel
        self.gamma = gamma
        self.tol = tol
        self.max_iter = max_iter

    # 定义 kernel_function函数
    def kernel_function(self, x1, x2):
        # kernel为linear条件下，返回x1x2矩阵相乘后的结果
        if self.kernel == 'linear':
            return np.dot(x1, x2)
        # kernel为rbf条件下，返回多项式计算结果
        elif self.kernel == 'rbf':
            return np.exp(-self.gamma * np.sum(np.square(x1 - x2)))
        else:
            raise ValueError("Unsupported kernel type")

    # 定义select_jrand函数
    def select_jrand(self, i, m):
        # 随机选择一个不等于i的j
        j = i
        while j == i:
            # np.random.randint()函数生成半开区间 [low, high) 内的随机整数。
            # 生成[0.m)之间的随机整数
            j = np.random.randint(0, m)
        return j

    def clip_alpha(self, aj, H, L):
        # 裁剪alpha值到[L, H]区间
        if aj > H:
            aj = H
        if aj < L:
            aj = L
        return aj

    def fit(self, X, y):
        self.X = X
        self.y = y
        # X.shape会输出两个数，X的行数和列数，分别赋值给m和n
        self.m, self.n = X.shape
        # 生成一个纯0矩阵alphas
        self.alphas = np.zeros(self.m)
        self.b = 0
        self.iter_count = 0

        # 缓存核矩阵（优化性能）
        self.K = np.zeros((self.m, self.m))
        for i in range(self.m):
            for j in range(self.m):
                self.K[i, j] = self.kernel_function(X[i], X[j])

        # SMO主循环
        for iter in range(self.max_iter):
            alpha_pairs_changed = 0
            for i in range(self.m):
                # 计算预测值和误差
                fxi = np.sum(self.alphas * y * self.K[:, i]) + self.b
                Ei = fxi - y[i]

                # 检查是否违反KKT条件
                if ((y[i] * Ei < -self.tol) and (self.alphas[i] < self.C)) or \
                        ((y[i] * Ei > self.tol) and (self.alphas[i] > 0)):

                    # 随机选择第二个alpha
                    j = self.select_jrand(i, self.m)

                    # 计算第二个alpha的预测值和误差
                    fxj = np.sum(self.alphas * y * self.K[:, j]) + self.b
                    Ej = fxj - y[j]

                    # 保存旧的alphas
                    alpha_i_old = self.alphas[i].copy()
                    alpha_j_old = self.alphas[j].copy()

                    # 计算L和H
                    if y[i] != y[j]:
                        L = max(0, self.alphas[j] - self.alphas[i])
                        H = min(self.C, self.C + self.alphas[j] - self.alphas[i])
                    else:
                        L = max(0, self.alphas[j] + self.alphas[i] - self.C)
                        H = min(self.C, self.alphas[j] + self.alphas[i])

                    if L == H:
                        continue

                    # 计算核函数的eta值
                    eta = 2.0 * self.K[i, j] - self.K[i, i] - self.K[j, j]
                    if eta >= 0:
                        continue

                    # 更新alpha[j]
                    self.alphas[j] -= y[j] * (Ei - Ej) / eta
                    self.alphas[j] = self.clip_alpha(self.alphas[j], H, L)

                    # 检查alpha[j]的变化是否足够大
                    if abs(self.alphas[j] - alpha_j_old) < 1e-5:
                        continue

                    # 更新alpha[i]
                    self.alphas[i] += y[i] * y[j] * (alpha_j_old - self.alphas[j])

                    # 计算两个常数项b1和b2
                    b1 = self.b - Ei - y[i] * (self.alphas[i] - alpha_i_old) * self.K[i, i] - \
                         y[j] * (self.alphas[j] - alpha_j_old) * self.K[i, j]

                    b2 = self.b - Ej - y[i] * (self.alphas[i] - alpha_i_old) * self.K[i, j] - \
                         y[j] * (self.alphas[j] - alpha_j_old) * self.K[j, j]

                    # 更新b
                    if 0 < self.alphas[i] < self.C:
                        self.b = b1
                    elif 0 < self.alphas[j] < self.C:
                        self.b = b2
                    else:
                        self.b = (b1 + b2) / 2.0

                    alpha_pairs_changed += 1

            # 如果没有更新任何alpha对，增加迭代计数
            if alpha_pairs_changed == 0:
                self.iter_count += 1
            else:
                self.iter_count = 0

            # 如果连续多次迭代没有更新，退出循环
            if self.iter_count >= 10:
                break

        # 保存支持向量
        self.sv_mask = self.alphas > 1e-5
        self.support_vectors = X[self.sv_mask]
        self.support_labels = y[self.sv_mask]
        self.support_alphas = self.alphas[self.sv_mask]

        return self

    def predict(self, X):
        result = np.zeros(len(X))
        for i in range(len(X)):
            prediction = 0
            for a, sv, sl in zip(self.support_alphas, self.support_vectors, self.support_labels):
                prediction += a * sl * self.kernel_function(X[i], sv)
            result[i] = prediction + self.b
        return np.sign(result)

这个类函数体现了支持向量机的运行过程，为了更好地说明这一过程，将在下一篇文章单独详细说明，这里主要介绍这个SVM算法的执行过程。

【2.5】算法应用

之后是对算法进行训练、使用和输出：

# 训练模型
svm = SVM(C=1.0, kernel='linear', max_iter=100)
svm.fit(X_train, y_train)

# 在测试集上评估模型
y_pred = svm.predict(X_test)
accuracy = np.mean(y_pred == y_test)
print(f"Model Accuracy: {accuracy:.2f}")

# 输出支持向量的数量
print(f"Number of Support Vectors: {len(svm.support_vectors)}")

【2.6】测试算例

实际上可以使用任何算例来测试：

# 示例预测
sample = np.array([[5.1, 3.5, 1.4, 0.2]])
sample = scaler.transform(sample)
prediction = svm.predict(sample)
print(f"Predicted Class: {'Setosa' if prediction == -1 else 'Versicolour'}")

实际上整个代码非常简单，就是一个矩阵实例调用了SVM类函数执行了操作。

此时的完整代码为：

import numpy as np
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# 加载并预处理数据
iris = datasets.load_iris()
# iris的前 100 个样本的 4 个特征（数值）（花萼长度、花萼宽度、花瓣长度、花瓣宽度）
X = iris.data[:100]  # 只取前两类（线性可分）
# iris.target共有150个样本，每个样本50个，前100个样本对应的标签分别是0和1
y = iris.target[:100]
# 如果y=0,输出-1，否则输出1
y = np.where(y == 0, -1, 1)  # 将标签转换为-1和1

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# 数据标准化
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)


class SVM:
    def __init__(self, C=1.0, kernel='linear', gamma=0.1, tol=1e-3, max_iter=100):
        # 定义一批self量存储数据
        self.C = C
        self.kernel = kernel
        self.gamma = gamma
        self.tol = tol
        self.max_iter = max_iter

    # 定义 kernel_function函数
    def kernel_function(self, x1, x2):
        # kernel为linear条件下，返回x1x2矩阵相乘后的结果
        if self.kernel == 'linear':
            return np.dot(x1, x2)
        # kernel为rbf条件下，返回多项式计算结果
        elif self.kernel == 'rbf':
            return np.exp(-self.gamma * np.sum(np.square(x1 - x2)))
        else:
            raise ValueError("Unsupported kernel type")

    # 定义select_jrand函数
    def select_jrand(self, i, m):
        # 随机选择一个不等于i的j
        j = i
        while j == i:
            # np.random.randint()函数生成半开区间 [low, high) 内的随机整数。
            # 生成[0.m)之间的随机整数
            j = np.random.randint(0, m)
        return j

    def clip_alpha(self, aj, H, L):
        # 裁剪alpha值到[L, H]区间
        if aj > H:
            aj = H
        if aj < L:
            aj = L
        return aj

    def fit(self, X, y):
        self.X = X
        self.y = y
        # X.shape会输出两个数，X的行数和列数，分别赋值给m和n
        self.m, self.n = X.shape
        # 生成一个纯0矩阵alphas
        self.alphas = np.zeros(self.m)
        self.b = 0
        self.iter_count = 0

        # 缓存核矩阵（优化性能）
        self.K = np.zeros((self.m, self.m))
        for i in range(self.m):
            for j in range(self.m):
                self.K[i, j] = self.kernel_function(X[i], X[j])

        # SMO主循环
        for iter in range(self.max_iter):
            alpha_pairs_changed = 0
            for i in range(self.m):
                # 计算预测值和误差
                fxi = np.sum(self.alphas * y * self.K[:, i]) + self.b
                Ei = fxi - y[i]

                # 检查是否违反KKT条件
                if ((y[i] * Ei < -self.tol) and (self.alphas[i] < self.C)) or \
                        ((y[i] * Ei > self.tol) and (self.alphas[i] > 0)):

                    # 随机选择第二个alpha
                    j = self.select_jrand(i, self.m)

                    # 计算第二个alpha的预测值和误差
                    fxj = np.sum(self.alphas * y * self.K[:, j]) + self.b
                    Ej = fxj - y[j]

                    # 保存旧的alphas
                    alpha_i_old = self.alphas[i].copy()
                    alpha_j_old = self.alphas[j].copy()

                    # 计算L和H
                    if y[i] != y[j]:
                        L = max(0, self.alphas[j] - self.alphas[i])
                        H = min(self.C, self.C + self.alphas[j] - self.alphas[i])
                    else:
                        L = max(0, self.alphas[j] + self.alphas[i] - self.C)
                        H = min(self.C, self.alphas[j] + self.alphas[i])

                    if L == H:
                        continue

                    # 计算核函数的eta值
                    eta = 2.0 * self.K[i, j] - self.K[i, i] - self.K[j, j]
                    if eta >= 0:
                        continue

                    # 更新alpha[j]
                    self.alphas[j] -= y[j] * (Ei - Ej) / eta
                    self.alphas[j] = self.clip_alpha(self.alphas[j], H, L)

                    # 检查alpha[j]的变化是否足够大
                    if abs(self.alphas[j] - alpha_j_old) < 1e-5:
                        continue

                    # 更新alpha[i]
                    self.alphas[i] += y[i] * y[j] * (alpha_j_old - self.alphas[j])

                    # 计算两个常数项b1和b2
                    b1 = self.b - Ei - y[i] * (self.alphas[i] - alpha_i_old) * self.K[i, i] - \
                         y[j] * (self.alphas[j] - alpha_j_old) * self.K[i, j]

                    b2 = self.b - Ej - y[i] * (self.alphas[i] - alpha_i_old) * self.K[i, j] - \
                         y[j] * (self.alphas[j] - alpha_j_old) * self.K[j, j]

                    # 更新b
                    if 0 < self.alphas[i] < self.C:
                        self.b = b1
                    elif 0 < self.alphas[j] < self.C:
                        self.b = b2
                    else:
                        self.b = (b1 + b2) / 2.0

                    alpha_pairs_changed += 1

            # 如果没有更新任何alpha对，增加迭代计数
            if alpha_pairs_changed == 0:
                self.iter_count += 1
            else:
                self.iter_count = 0

            # 如果连续多次迭代没有更新，退出循环
            if self.iter_count >= 10:
                break

        # 保存支持向量
        self.sv_mask = self.alphas > 1e-5
        self.support_vectors = X[self.sv_mask]
        self.support_labels = y[self.sv_mask]
        self.support_alphas = self.alphas[self.sv_mask]

        return self

    def predict(self, X):
        result = np.zeros(len(X))
        for i in range(len(X)):
            prediction = 0
            for a, sv, sl in zip(self.support_alphas, self.support_vectors, self.support_labels):
                prediction += a * sl * self.kernel_function(X[i], sv)
            result[i] = prediction + self.b
        return np.sign(result)


# 训练模型
svm = SVM(C=1.0, kernel='linear', max_iter=100)
svm.fit(X_train, y_train)

# 在测试集上评估模型
y_pred = svm.predict(X_test)
accuracy = np.mean(y_pred == y_test)
print(f"Model Accuracy: {accuracy:.2f}")

# 输出支持向量的数量
print(f"Number of Support Vectors: {len(svm.support_vectors)}")

# 示例预测
sample = np.array([[5.1, 3.5, 1.4, 0.2]])
sample = scaler.transform(sample)
prediction = svm.predict(sample)
print(f"Predicted Class: {'Setosa' if prediction == -1 else 'Versicolour'}")

【3】标准化理解

在对代码通透理解之前，先理解两简单的小函数：X_train = scaler.fit_transform(X_train)和X_test = scaler.transform(X_test)。

这两个小函数都是sklearn中的子函数。

首先要追溯scaler的来历，它是scaler = StandardScaler()运算得到的。StandardScaler()是sklearn中的标准化器，能让不符合标准正态分布的数据转移到标准正态分布的轨道上来。

X_train = scaler.fit_transform(X_train)实际上包括两个动作，先计算X_train的均值和标准差，然后再使用这两个参数来标准化X_train；

X_test = scaler.transform(X_test)。就只包括一个动作，直接使用X_train的均值和标准差来标准化X_test。

给一个测试算例：

from sklearn.preprocessing import StandardScaler
import numpy as np

# 原始数据
X_train = np.array([[1, 2], [3, 4], [5, 6]])  # 训练集
X_test = np.array([[7, 8], [9, 10]])          # 测试集

# 标准化流程
scaler = StandardScaler()
scaler.fit(X_train)  # 计算并存储训练集的统计量

# 手动计算标准化参数
mu = X_train.mean(axis=0)  # 均值 [3, 4]
sigma = X_train.std(axis=0)  # 标准差 [1.633, 1.633]

# 手动标准化测试集
X_test_manual = (X_test - mu) / sigma

# 使用transform方法标准化测试集
X_test_transform = scaler.transform(X_test)

# 使用fit_transform和transform方法标准化测试集
X_test_scaler= scaler.fit_transform(X_train)
X_test_fit_transform= scaler.transform(X_test)
# 验证结果是否一致
print("手动计算结果:\n", X_test_manual)
print("\ntransform方法结果:\n", X_test_transform)
print("\nX_test_scaler方法结果:\n", X_test_scaler)
print("\nfit_transform方法结果:\n", X_test_fit_transform)
print("\n结果是否一致:", np.allclose(X_test_manual, X_test_transform))
print("\n结果是否一致:", np.allclose(X_test_fit_transform, X_test_transform))

这里给出了三种标准化方法：手动、transform和fit_transform+transform。

transform和fit_transform+transform中都用了transform方法，区别在于：transform方法调用了手动计算获得的均值mu和方差sigma；fit_transform+transform方法先通过fit_transform计算获得均值mu和方差sigma，然后transform直接调用这两个参数。

由于手动计算和fit_transform计算得到的均值mu和方差sigma是相等的，所以计算结果相等，完整的输出内容为：

手动计算结果:
[[2.44948974 2.44948974]
[3.67423461 3.67423461]]

transform方法结果:
[[2.44948974 2.44948974]
[3.67423461 3.67423461]]

X_test_scaler方法结果:
[[-1.22474487 -1.22474487]
[ 0. 0. ]
[ 1.22474487 1.22474487]]

fit_transform方法结果:
[[2.44948974 2.44948974]
[3.67423461 3.67423461]]

结果是否一致: True

结果是否一致: True

【4】细节说明

sklearn中的标准化操作是使得数据回归到标准正态分布，不是其他分布。

【5】总结

学习了支持向量机算法的一个简单python示例。

python学智能算法（十六）|机器学习支持向量机简单示例