机器学习概论 SMO算法实现（实验三）

一、实验目的
1、熟悉SMO算法的求解过程；
2、掌握SMO算法的代码实现；
3、掌握SVM的分类方法；
4、熟悉numpy的使用。
二、实验设备
计算机：CPU四核i7 6700处理器；内存8G； SATA硬盘2TB硬盘； Intel芯片主板；集成声卡、千兆网卡、显卡； 20寸液晶显示器。
编译环境：python解释器、Pycharm编辑器
三、实验内容
1、新建项目和文件，并导入numpy
（1）打开Pycharm，新建项目，并在该项目下新建文件svm_smo.py。
（2）导入一些编程中需要的包。

import numpy as np
import matplotlib.pyplot as plt
import math
import random

2、构造一个函数data_producer
（1）定义函数头data_producer，该函数用于生成二维正态分布的点，共有三个参数，miu，sigma和sample_no，分别表示二维正态分布的均值，标准差和产生样本个数。

def data_producer(miu, sigma, sample_no):
    data = np.random.multivariate_normal(miu, sigma, sample_no)
    return data

3、编写一个类SVMSMO
（1）构造一个类SVMSMO，并添加构造函数，需要训练数据集、测试数据集、最大迭代次数和核函数的参数作为输入；
（2）构造类方法get_j用于获取第二个alpha的索引；
（3）构造类方法clip_alpha用于切割更新的alpha_j；
（4）构造类方法kernel_function用于返回核函数的值；
（5）构造类方法smo用于计算alpha和b；
（6）构造类方法svm用于计算分类结果。

class SVMSMO:
    def __init__(self, data_train, data_test, maxcihsu):
        self.data_train = data_train
        self.data_test = data_test
        self.maxcihsu = maxcihsu   #迭代次数
        self.alpha = None  #用于保存alpha值
        self.b=0        #用于保存b的值


    def get_j(self, i):
        while True:
            j = random.choice(range(160)) #随机选择第二个变量
            if j != i:
                break
        return j


    def clip_alpha(self, a_j, l, h):   #用于剪辑的解
        if a_j > h:
            a_j = h
        elif a_j < l:
            a_j = l
        else:
            a_j = a_j
        return a_j


    def kernel_function(self, x_i, x_j):  #线性核函数
        return np.dot(x_i, x_j)


    def gaisidange(self, x_i, x_j):   #高斯核函数
        a = 1.1
        c = np.array(x_i)-np.array(x_j)
        if len(x_j.shape) == 1:      #因为是向量，所以根据传入的参数要运用不同的计算方式
            kern = np.exp((-1*np.dot(c,c.T))/(2*a**2))
        else:
            kern = np.exp((-1*np.diag(np.dot(c,c.T)))/(2*a**2))  #np.diag表示对角线，意思是取对角线上的值
        return kern


    def smo(self):
        alphas = np.zeros(int(2*num*0.8))    #a为行向量形式(1600,)
        b = 0; C = 0.8; iter = 0; renew = 1      #初始化 偏置b，松弛量c,循环次数iter, 更新次数renew
        while (iter <= self.maxcihsu) and renew > 0: #迭代次数
            renew = 0
            for i in range(int(2*num*0.8)):
                fxi=np.dot(np.multiply(alphas, data_train[:, 2]), self.gaisidange(data_train[i, :-1], data_train[:, :-1])) + b
                Ei=fxi - float(data_train[i, 2])      #计算E(i)
                if (alphas[i] < C and np.sum(np.multiply(data_train[:, 2], fxi)) <= 1) or  (alphas[i] > 0 and np.sum(np.multiply(data_train[:, 2], fxi)) >= 1) or \
                        ((alphas[i]==C) or np.sum(np.multiply(data_train[:,2],fxi)) == 1 and alphas[i] == 0):
                    j= self.get_j(i)         #选择a_j
                    Ej = np.dot(np.multiply(alphas, data_train[:, 2]), self.gaisidange(data_train[j, :-1], data_train[:, : -1]))+b - float(data_train[j, 2])# 预测值和真实输出之差
                    alphaIold = alphas[i].copy()
                    alphaJold = alphas[j].copy()     #计算e(j)
                    if data_train[i, 2] != data_train[j, 2]:# 这里是对SMO最优化问题的子问题的约束条件的分析
                        #当y1等于y2时，即两个都为正1或者都为负1，可以得到：
                        #alpha[i]+alpha[j]=s
                        #L=max(0,s-C)
                        #H=min(c,s)
                        #若y1!=y2,
                        L = max(0, alphas[j] - alphas[i]) # L和H分别是alpha所在的对角线端点的界
                        H = min(C, C + alphas[j] - alphas[i])  # 调整alphas[j]位于0到c之间
                    else:
                        #若y1=y2
                        L = max(0, alphas[j] + alphas[i] - C)
                        H = min(C, alphas[j] + alphas[i])
                    if L==H:
                        continue   # L=H停止本次循环
                    # 是一个中间变量：eta=2xi*xi-xixi-xjxj，是alphas[j]的最优修改量
                    eta =self.gaisidange(data_train[i,:-1],data_train[i,:-1])+self.gaisidange(data_train[j,:-1],data_train[j,:-1])\
                    -2*self.gaisidange(data_train[i,:-1],data_train[j,:-1])                       #计算eta
                    alphaJold = alphas[j]                                                       #更新alpha值
                    alphaJnew =alphaJold + data_train[j, 2]*(Ei-Ej)/eta
                    alphaJnew_clipped = self.clip_alpha(alphaJnew, L, H)
                    alphas[j] = alphaJnew_clipped
                    alphasIold = alphas[i]
                    alphasInew = alphasIold+data_train[i, 2]*data_train[j,2]*(alphaJold-alphaJnew_clipped)
                    alphas[i] = alphasInew
                    b_old = b          #更新b值
                    b_i_new = b_old-Ej-data_train[i, 2]*(alphasInew-alphasIold)*self.gaisidange(data_train[i, :-1], data_train[i, :-1])\
                    -data_train[j, 2]*(alphaJnew-alphaJold)*self.gaisidange(data_train[j, :-1], data_train[i, :-1])
                    b_j_new = b_old-Ej-data_train[i, 2]*(alphasInew-alphasIold)*self.gaisidange(data_train[i, :-1], data_train[j, :-1])\
                    -data_train[j, 2]*(alphaJnew-alphaJold)*self.gaisidange(data_train[j, :-1], data_train[j, :-1])
                    if 0 < alphas[i] < C:
                        b = b_i_new
                    elif 0 < alphas[i] < C:
                        b = b_j_new
                    else:
                        b = (b_i_new+b_j_new)/2
                    renew += 1
            iter += 1
            print(iter)
        self.b = b
        self.alpha = alphas.copy()
        print("b的值",self.b)
    def svm(self):
        pred = np.zeros(int(2*num*0.2))
        for i in range(int(2*num*0.2)):
            fx=np.dot(np.multiply(self.alpha,data_train[:,2]),self.gaisidange(data_test[i,:-1],data_train[:,:-1]))+self.b   #根据x的值求fx的值，fx的值有正有负，可对应label的值
            if fx>=0:
                pred[i]=1
            if fx<0:
                pred[i]=-1
        self.data_test=np.dot(pred,self.data_test)  #根据label的值计算横纵坐标，主要是让其更加有序化
        ccc=data_test[:,2]
        sss=pred
        accuracy=sum(data_test[:,2]==pred)/(int(2*num*0.2))#计算准确率
        return accuracy

4、编写if name”main”:
（1）构造if name”main”:并在其中调用data_producer获取数据。
（2）声明SVMSMO的一个对象my_svmsmo，并使用该对象的类方法对上一步产生的数据进行分类，打印分类结果。

if __name__ == "__main__":
    num = 100
    data_p = data_producer([2, 2], [[1, -0.1], [-.1, 1]], num)
    data_n = data_producer([-2, -2], [[1, -0.1], [-.1, 1]], num)      #num个负例
    lable_p = np.ones((num, 1))
    data_p = np.hstack((data_p, lable_p))      #横向添加标签
    lable_n = -1*np.ones((num, 1))           #生成100行1列的二维矩阵，并把其赋值为-1
    data_n = np.hstack((data_n, lable_n))
    data_train = np.vstack((data_p[0:int(num*0.8), :], data_n[0:int(num*0.8), :])) #将样本数据分为测试集和训练集并进行纵向合并
    data_test = np.vstack((data_p[int(num*0.8):, :], data_n[int(num*0.8):, :]))
    my_svmsmo = SVMSMO(data_train, data_test, 100)
    my_svmsmo.smo()
    accuracy = my_svmsmo.svm()
    print(accuracy)
    plt.figure()
    plt.subplot(121)
    plt.plot(data_train[:int(num*0.8), 0], data_train[:int(num*0.8), 1], ".", color='y')
    plt.plot(data_train[int(num*0.8):, 0], data_train[int(num*0.8):, 1], "*", color='b')
    plt.subplot(122)
    plt.plot(data_test[0:int(num*0.2), 0], data_test[0:int(num*0.2), 1], ".", color='y')
    plt.plot(data_test[int(num*0.2):, 0], data_test[int(num*0.2):, 1], "*", color='b')
    plt.show()

实验截图：

五、实验总结
本次实验又一次练习了numpy和matplotlib的使用方法，并且在一定程度上熟悉了SMO算法的求解过程，掌握SMO算法的代码实现和SVM的分类方法；虽然掌握的不够熟练，但是基本明白了实现的过程，有待进一步去学习，渗透。