• 实验环境

使用python编程实现，并在Mindspore框架下实现。

• 实验内容

题目： 鸢尾花数据集（Iris data set）是模式识别中著名的数据集。本实验通过花萼（sepal）和花瓣（petal）的长和宽，建立SVM分类器来判断样本属于山鸢尾（Iris Setosa）、变色鸢尾（Iris Versicolor）还是维吉尼亚鸢尾（Iris Virginica）。请按要求完成实验。

数据集： 文件iris.txt为该实验的数据集，包含150个样本，对应数据集的每行数据。每行数据包含每个样本的四个特征（按顺序分别为花萼长度、花萼宽度、花瓣长度、花瓣宽度）和样本的类别信息（Setosa、Versicolor、 Virginica中的一种）。

实验要求： 

（1）建立SVM分类器并用交叉验证法进行分析。 

（2）利用PCA降维，将数据转化为二维，然后绘制出分类决策边界。

iris.txt数据集：（只截取一部分）

运行程序结果图：

解释关键代码：

1. loadData 函数：读取数据文件，将每行数据分割并转换为浮点数，根据类别标签设置不同的标签值。

• 初始化变量：初始化五个空列表，用于存储数据点的特征值、三个类别的标签值和类别名称。
• 读取文件：打开名为 'iris.txt' 的文件，并将其文件对象存储在变量 fr 中。
• 逐行读取和解析数据：逐行读取文件内容，使用 strip() 方法去除每行末尾的换行符，然后使用 split(',') 方法按逗号分割每行数据，得到一个包含该行数据的列表 lineArr。
• 提取特征数据：将每行的前四个值（鸢尾花的四个特征）转换为浮点数，并添加到 dataMat 列表中。
• 提取标签数据：根据第五个值（类别名称），如果类别是'Iris-setosa'，则在 labelMat1 中添加1，否则添加-1。这里使用1和-1来表示两个类别。类似地，为'Iris-versicolor'类别在labelMat2中添加1或-1。为'Iris-virginica'类别在labelMat3中添加1或-1。将类别名称添加到ylabel列表中。

2.pca 函数：pca 函数用于数据降维，以减少计算复杂度并提高算法效率。

• 从数据矩阵中减去每列的平均值，以中心化数据。
• 计算中心化后数据的协方差矩阵。参数 rowvar=0 表示每行是一个变量。
• 计算协方差矩阵的特征值和特征向量。
• 对特征值进行排序，得到从大到小的索引。
• 选择最大的 topNfeat 个特征值对应的特征向量。
• 根据选择的特征值索引，获取对应的特征向量。
• 将原始数据投影到这些特征向量上，得到降维后的数据矩阵。
• 返回降维后的数据矩阵。

3.selectJrand 函数：随机选择一个与当前索引不同的索引，用于SMO算法中的优化。

4.clipAlpha 函数：确保拉格朗日乘子的值在允许的范围内。

5.smoSimple 函数：实现SMO算法，通过迭代优化拉格朗日乘子，求解SVM的优化问题。

• 初始化变量：将类别标签转换为矩阵形式，并转置。初始化偏置项 b，并获取数据矩阵的行数m（样本数量）和列数n（特征数量）。初始化拉格朗日乘子矩阵alphas，初始值都为0。初始化迭代次数计数器。
• 迭代优化：

1. 初始化一个标志，用于记录在当前迭代中是否有拉格朗日乘子被优化。
2. 选择拉格朗日乘子对进行优化：遍历每个样本，尝试找到可以优化的拉格朗日乘子对。计算第 i 个样本的预测值和误差。检查第 i 个样本的误差是否超过了容忍度，并且是否还有优化空间。随机选择另一个样本 j，用于与样本 i 配对优化。
3. 计算新的拉格朗日乘子值：计算第 j 个样本的预测值和误差。保存旧的拉格朗日乘子值，以便后续更新。根据类别标签和当前的拉格朗日乘子值计算新的拉格朗日乘子的上下界。如果上下界相同，则跳过当前的迭代。计算 eta，它是两个样本之间的核函数值的差。如果 eta 非负，则跳过当前的迭代。更新第 j 个样本的拉格朗日乘子值，并确保它在上下界之间。如果拉格朗日乘子的变化非常小，则跳过当前的迭代。更新第 i 个样本的拉格朗日乘子值。计算新的偏置项 b。根据拉格朗日乘子的值选择新的偏置项。如果找到了可以优化的拉格朗日乘子对，则增加变化计数器。

• 检查迭代是否结束：如果当前迭代中没有找到可以优化的拉格朗日乘子对，则增加迭代次数；否则重置迭代次数

6.calcWs 函数：根据优化后的拉格朗日乘子和数据点计算权重向量。

使用shape函数获取dataMatrix的行数m（样本数量）和列数n（特征数量）。初始化的零向量w，用于存储最终计算出的权重。遍历每个样本，i从0到m-1。对于每个样本，计算其对应的权重贡献。

7.show 函数：绘制PCA降维后的数据点和决策边界，用不同颜色区分不同类别。

• 初始化六个空列表，用于存储三个类别的数据点坐标。
• 生成一个从 -4 到 5 的数组，步长为 0.1，用于绘制决策边界。
• 创建一个新的图形和子图，并把 dataxArr 转换为矩阵形式。
• 遍历所有样本，根据类别标签分组存储坐标。
• 计算两个SVM决策边界的值。调整 z 和 z1 的形状以匹配 xx 和 yy 的网格形状。
• 根据决策边界的值给网格上的每个点分配类别标签。
• 使用等高线填充函数 contourf 绘制决策边界。
• 绘制每个类别的数据点。
• 设置图形的标签、标题、坐标轴范围，并显示图形。

8.主程序：调用 loadData 函数加载鸢尾花数据集。使用PCA将数据降维到2维。使用SMO算法求解两个类别的SVM参数。计算两个类别的权重向量。调用 show 函数显示分类结果。

#!/usr/bin/python
#-*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
import matplotlib.animation as ai
import numpy as np
import time

def loadData():    #加载数据函数
    dataMat = []
    labelMat1 = []
    labelMat2 = []
    labelMat3 = []
    ylabel = []
    fr = open('iris.txt')
    for line in fr.readlines():
        lineArr = line.strip().split(',')
        dataMat.append([float(lineArr[0]), float(lineArr[1]), float(lineArr[2]), float(lineArr[3])])
        if(lineArr[4]=='Iris-setosa'):
          labelMat1.append(float(1))
        else:
          labelMat1.append(float(-1))
        if(lineArr[4]=='Iris-versicolor'):
          labelMat2.append(float(1))
        else:
          labelMat2.append(float(-1))
        if(lineArr[4]=='Iris-virginica'):
          labelMat3.append(float(1))
        else:
          labelMat3.append(float(-1))
        ylabel.append(lineArr[4])
    return dataMat,labelMat1,labelMat2,labelMat3,ylabel


def pca(dataMat, topNfeat):              #pca降维
    re = dataMat - mean(dataMat, axis = 0)   #去除平均值 
    covMat = cov(re,rowvar=0) #计算协防差矩阵  
    eVals, eVects = linalg.eig(mat(covMat))  
    eValInd = argsort(eVals)  
    #从小到大对N个值排序  
    eValInd = eValInd[: -(topNfeat + 1) : -1]  
    reVects = eVects[:, eValInd]
    lowDataMat = re * reVects   #转换到降维空间
    return lowDataMat

def selectJrand(i,m):                         #随机选择alpha
    j=i             #排除i
    while (j==i):
          j = int(random.uniform(0,m))
    return j
     
def clipAlpha(aj,H,L):                       #规范alpha的值
    if aj > H:
       aj = H
    if L > aj:
       aj = L
    return aj

def smoSimple(dataMatrix, classLabels, C, toler, maxIter):       #简单SMO算法求解
    labelMat = mat(classLabels).T
    b = -1; m,n = shape(dataMatrix) 
    alphas = mat(zeros((m,1)))
    iter = 0
    while (iter < maxIter):
        alphaPairsChanged = 0   #alpha是否已经进行了优化
        for i in range(m):
            #   w = alpha * y * x;  f(x_i) = w^T * x_i + b
            # 预测的类别
            fXi = float(multiply(alphas,labelMat).T*dataMatrix*dataMatrix[i,:].T) + b    
            Ei = fXi - float(labelMat[i])   #得到误差，如果误差太大，检查是否可能被优化
            #必须满足约束
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)): 
                j = selectJrand(i,m)
                fXj = float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy()                
                if (labelMat[i] != labelMat[j]):                                          
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: 
                   continue
                # Eta = -(2 * K12 - K11 - K22)，且Eta非负，此处eta = -Eta则非正
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0:
                   continue
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j],H,L)
                  #如果内层循环通过以上方法选择的α_2不能使目标函数有足够的下降，那么放弃α_1
                if (abs(alphas[j] - alphaJold) < 0.00001): 
                    continue
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                alphaPairsChanged += 1
        if (alphaPairsChanged == 0): iter += 1
        else: iter = 0
    return b,alphas

def calcWs(alphas,dataMatrix, labelMat):      #求出参数w
    m,n = shape(dataMatrix) 
    w = zeros((n,1))
    for i in range(m):
        w += multiply(alphas[i]*labelMat[i],dataMatrix[i,:].T)
    return w


def show(dataxArr, ydata,b1, b3):           #显示函数
    n = len(dataxArr)
    xcord1 = [];ycord1 = []
    xcord2 = [];ycord2 = []
    xcord3 = [];ycord3 = []
    x = arange(-4,5,0.1)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    dataxArr = mat(dataxArr)
  
    #原始数据
    for i in range(n):
        if ydata[i]=='Iris-setosa':
            xcord1.append(dataxArr[i,0]);ycord1.append(dataxArr[i,1])
        elif ydata[i]=='Iris-versicolor':
            xcord2.append(dataxArr[i,0]);ycord2.append(dataxArr[i,1])
        else:
            xcord3.append(dataxArr[i,0]);ycord3.append(dataxArr[i,1])
    
    #分类边界
    xx,yy= np.meshgrid(x,x)
    z = xx * w1[0][0] + yy * w1[1][0] + float(b1)
    z1 = xx * w3[0][0] + yy * w3[1][0] + float(b3)
    z.reshape(xx.shape)
    z1.reshape(xx.shape)
    m ,n = shape(z)
    for i in range(m):
        for j in range(n):
            if(z[i][j]>0):
               z[i][j]=1
            elif(z1[i][j]>0):
               z[i][j]=26
            else:
               z[i][j]=50

    #图形的显示
    plt.contourf(xx, yy, z, 10, cmap = plt.cm.Spectral)
    ax.scatter(xcord1,ycord1,s=30,c='b',marker='o',label="Iris-setosa")
    ax.scatter(xcord2,ycord2,s=30,c='g',marker='o',label="Iris-versicolor")
    ax.scatter(xcord3,ycord3,s=30,c='r',marker='s',label='Iris-virginica')
    plt.xlabel('X1',size=25)                  # 横坐标
    plt.ylabel('X2',size=25)                  # 纵坐标
    plt.xlim(xx.min(), xx.max())
    plt.ylim(yy.min(), yy.max())
    plt.title ('SVM_PCA',size=30)             # 标题
    plt.xticks(())
    plt.yticks(())
    plt.show()  
  
xdata,ydata1,ydata2,ydata3,ylabe = loadData() #加载数据
X = pca(xdata,2) #将数据降为2维
b1 , alphas1 = smoSimple(X,ydata1,0.8,0.0001,40)  #求解第一个分隔平面的参数
b3 , alphas3 = smoSimple(X,ydata3,0.8,0.0001,40)  #求解第二个分隔平面的参数
w1 = calcWs(alphas1,X,ydata1) #求解w1
w3 = calcWs(alphas3,X,ydata3) #求解w3
show(X, ylabe,b1,b3)  #显示图形

源码：实现了一个支持向量机（SVM）分类器，用于对鸢尾花数据集进行分类。代码使用了主成分分析（PCA）进行降维，并应用了序列最小优化（SMO）算法来求解SVM的优化问题。

• 实验收获

1. 理解SVM分类器：学习如何建立支持向量机（SVM）分类器，了解其工作原理，包括边际最大化和核技巧。
2. 掌握交叉验证：通过交叉验证，学会如何评估模型的泛化能力，以及如何选择合适的模型参数。
3. PCA降维技术：学习主成分分析（PCA）降维技术，理解其在数据预处理中的重要性，特别是在高维数据集中。
4. 决策边界的理解：通过绘制分类决策边界，直观地理解SVM如何区分不同类别的数据。