【数模修炼之旅】08 支持向量机模型 深度解析（教程+代码）
接下来 C君将会用至少30个小节来为大家深度解析数模领域常用的算法，大家可以关注这个专栏，持续学习哦，对于大家的能力提高会有极大的帮助。
1 支持向量机模型介绍及应用
支持向量机（SVM, Support Vector Machines）是一种强大的监督学习算法，用于分类和回归任务。它在高维空间中寻找一个超平面，以最大化不同类别之间的边界。下面将详细介绍SVM的基本原理以及它在数学建模中的应用。
1.1 支持向量机的基本原理
- 分类器定义： 支持向量机通过一个超平面将数据分为两类。这个超平面是由数据中最难区分的点，即支持向量，确定的。
- 最大化间隔： SVM的目标是找到一个最大化类别间隔的超平面。间隔定义为到最近数据点的最小距离，这些点即为支持向量。
- 核技巧： 当数据不是线性可分的时候，SVM通过一个转换函数（核函数）将数据映射到一个更高维的空间，在这个新空间中寻找超平面。常用的核函数包括线性核、多项式核、径向基函数核（RBF）等。
- 软间隔与正则化： 为了提高模型的泛化能力，SVM引入了软间隔的概念，允许一些数据点在超平面的错误一侧，通过引入惩罚参数C来平衡间隔宽度与分类误差。
1.2 支持向量机在数学建模中的应用
- 模式识别：
图像识别：SVM可以用于识别图像中的对象，如车辆、人脸等。
手写识别：利用SVM对手写数字或字母进行分类。
- 生物信息学：
蛋白质分类：使用SVM根据蛋白质的特征对其进行功能分类。
基因表达数据分析：通过SVM分析基因表达数据，帮助识别疾病相关的基因。
- 金融分析：
信用评分：SVM用于评估个人的信用等级，预测贷款违约的可能性。
股市预测：利用历史数据，通过SVM模型预测股票价格的走势。
- 语音识别：
SVM能够用于区分不同的语音模式，实现有效的语音识别功能。
- 环境科学：
气候变化预测：应用SVM模型在气候数据上，预测未来的气候变化。
2 支持向量机模型的基本步骤
支持向量机（SVM）的基本步骤可以从它的原理出发，详细解释如何构建模型和进行数据分类。这个过程通常包括几个关键步骤，从数据准备到模型训练，再到最终的模型评估和应用。
步骤1: 数据准备
- 数据收集：首先需要收集足够的数据来训练模型。这些数据应该代表了你想要模型解决的问题。
- 特征选择：根据问题的具体性质，选择最能代表数据特性的特征。
- 数据预处理：包括缩放（如归一化或标准化）、处理缺失值、去除噪声等，以提高模型的性能和准确度。
步骤2: 选择核函数
- 线性可分问题：如果数据是线性可分的，通常选择线性核。
- 非线性问题：对于非线性数据，可以选择多项式核、径向基函数（RBF）核或其他更复杂的核函数。核函数的选择依赖于数据的分布和特性。
步骤3: 构建SVM模型
- 模型配置：设置SVM参数，如正则化参数C（控制误差和决策面间隔的权衡）、核函数参数（如RBF核的γ）等。
- 训练模型：使用训练数据集训练SVM模型。在这个过程中，SVM算法会尝试找到最优的超平面，即能够最大化正负样本间隔的平面。
步骤4: 模型优化
- 交叉验证：使用交叉验证（如k-fold交叉验证）来评估不同参数的模型性能，从而选择最佳的参数组合。
- 参数调整：基于交叉验证结果调整模型参数，如调整C值或核参数，以达到更好的分类效果。
步骤5: 模型评估
- 测试集评估：使用独立的测试数据集评估模型的准确度、召回率、F1分数等性能指标。
- 性能分析：分析误差类型，查看模型在哪些类型的数据上表现不好，可能需要回到数据准备阶段调整特征或数据处理方式。
步骤6: 模型部署与应用
- 模型部署：将训练好的模型部署到实际应用中，如在线预测系统、产品推荐系统等。
- 持续监控：监控模型的表现，确保模型在新数据上依然有效。根据需要进行模型更新。
3 支持向量机模型代码（matlab+python）
3.1 python
from numpy import *
import random
import matplotlib.pyplot as plt
import numpydef kernelTrans(X,A,kTup): # 核函数(此例未使用)m,n=shape(X)K = mat(zeros((m,1)))if kTup[0] =='lin':K=X*A.Telif kTup[0]=='rbf':for j in range(m):deltaRow = X[j,:]-AK[j]=deltaRow*deltaRow.T # ||w||^2 = w^T * wK =exp(K/(-1*kTup[1]**2)) # K = e^(||x-y||^2 / (-2*sigma^2))else:raise NameError("Houston we Have a problem --")return Kclass optStruct:def __init__(self,dataMain,classLabel,C,toler,kTup):self.X = dataMain # 样本矩阵self.labelMat = classLabelself.C = C # 惩罚因子self.tol = toler # 容错率self.m = shape(dataMain)[0] # 样本点个数self.alphas = mat(zeros((self.m,1))) # 产生m个拉格郎日乘子，组成一个m×1的矩阵self.b =0 # 决策面的截距self.eCache = mat(zeros((self.m,2))) # 产生m个误差 E=f(x)-y ,设置成m×2的矩阵，矩阵第一列是标志位，标志为1就是E计算好了，第二列是误差E# self.K = mat(zeros((self.m,self.m)))# for i in range(self.m): # K[,]保存的是任意样本之间的相似度(用高斯核函数表示的相似度)# self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup)def loadDataSet(filename): # 加载数据dataMat = []labelMat = []fr = open(filename)for line in fr.readlines():lineArr = line.split()dataMat.append([float(lineArr[0]),float(lineArr[1])])labelMat.append(float(lineArr[2])) # 一维列表return dataMat, labelMatdef selectJrand(i, m): # 随机选择一个不等于i的下标j =iwhile(j==i):j = int(random.uniform(0,m))return jdef clipAlpha(aj, H,L):if aj>H: # 如果a^new 大于上限值，那么就把上限赋给它aj = Hif L>aj: # 如果a^new 小于下限值，那么就把下限赋给它aj = Lreturn ajdef calcEk(oS, k): # 计算误差E, k代表第k个样本点，它是下标，oS是optStruct类的实例# fXk = float(multiply(oS.alphas,oS.labelMat).T * oS.K[:,k] + oS.b) # 公式f(x)=sum(ai*yi*xi^T*x)+bfXk = float(multiply(oS.alphas,oS.labelMat).T * (oS.X*oS.X[k,:].T)) +oS.bEk = fXk - float(oS.labelMat[k]) # 计算误差 E=f(x)-yreturn Ekdef selectJ(i, oS, Ei): # 选择两个拉格郎日乘子,在所有样本点的误差计算完毕之后，寻找误差变化最大的那个样本点及其误差maxK = -1 # 最大步长的因子的下标maxDeltaE = 0 # 最大步长Ej = 0 # 最大步长的因子的误差oS.eCache[i] = [1,Ei]valiEcacheList = nonzero(oS.eCache[:,0].A)[0] # nonzero结果是两个array数组，第一个数组是不为0的元素的x坐标，第二个数组是该位置的y坐标# 此处寻找误差矩阵第一列不为0的数的下标print("valiEcacheList is {}".format(valiEcacheList))if (len(valiEcacheList))>1:for k in valiEcacheList: # 遍历所有计算好的Ei的下标，valiEcacheLIst保存了所有样本点的E，计算好的有效位置是1,没计算好的是0if k == i:continueEk = calcEk(oS,k)deltaE = abs(Ei-Ek) # 距离第一个拉格朗日乘子a1绝对值最远的作为第二个朗格朗日乘子a2if deltaE>maxDeltaE:maxK = k # 记录选中的这个乘子a2的下标maxDeltaE = deltaE # 记录他俩的绝对值Ej = Ek # 记录a2此时的误差return maxK, Ejelse: # 如果是第一次循环，随机选择一个alphasj = selectJrand(i, oS.m)# j = 72Ej = calcEk(oS, j)return j,Ejdef updateEk(oS, k):Ek = calcEk(oS, k)oS.eCache[k] = [1,Ek] # 把第k个样本点的误差计算出来，并存入误差矩阵，有效位置设为1def innerL(i, oS):Ei = calcEk(oS, i) # KKT条件， 若yi*(w^T * x +b)-1<0 则 ai=C 若yi*(w^T * x +b)-1>0 则 ai=0print("i is {0},Ei is {1}".format(i,Ei))if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):j,Ej = selectJ(i,oS,Ei)print("第二个因子的坐标{}".format(j))alphaIold = oS.alphas[i].copy() # 用了浅拷贝， alphaIold 就是old a1,对应公式alphaJold = oS.alphas[j].copy()if oS.labelMat[i] != oS.labelMat[j]: # 也是根据公式来的，y1 不等于 y2时L = max(0,oS.alphas[j] - oS.alphas[i])H = min(oS.C, oS.C+oS.alphas[j]-oS.alphas[i])else:L = max(0,oS.alphas[j]+oS.alphas[i]-oS.C)H = min(oS.C,oS.alphas[j]+oS.alphas[i])if L==H: # 如果这个j让L=H,i和j这两个样本是同一类别，且ai=aj=0或ai=aj=C，或者不同类别，aj=C且ai=0# 当同类别时 ai+aj = 常数 ai是不满足KKT的，假设ai=0,需增大它，那么就得减少aj，aj已经是0了，不能最小了，所以此情况不允许发生# 当不同类别时 ai-aj=常数，ai是不满足KKT的，ai=0,aj=C,ai需增大，它则aj也会变大，但是aj已经是C的不能再大了，故此情况不允许print("L=H")return 0# eta = 2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j] # eta=K11+K22-2*K12eta = 2.0*oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].Tif eta >= 0: # 这里跟公式正好差了一个负号，所以对应公式里的 K11+K22-2*K12 <=0，即开口向下，或为0成一条直线的情况不考虑print("eta>=0")return 0oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta # a2^new = a2^old+y2(E1-E2)/etaprint("a2 归约之前是{}".format(oS.alphas[j]))oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) # 根据公式，看看得到的a2^new是否在上下限之内print("a2 归约之后is {}".format(oS.alphas[j]))# updateEk(oS,j) # 把更新后的a2^new的E更新一下if abs(oS.alphas[j]-alphaJold)<0.00001:print("j not moving enough")return 0oS.alphas[i] +=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j]) # 根据公式a1^new = a1^old+y1*y2*(a2^old-a2^new)print("a1更新之后是{}".format(oS.alphas[i]))# updateEk(oS,i)# b1^new = b1^old+(a1^old-a1^new)y1*K11+(a2^old-a2^new)y2*K12-E1# b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i]-oS.labelMat[j]*\# (oS.alphas[j]-alphaJold)*oS.K[i,j]b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T-oS.labelMat[j]* \(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T# b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]-oS.labelMat[j]* \# (oS.alphas[j]-alphaJold)*oS.K[j,j]b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T-oS.labelMat[j]* \(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].TupdateEk(oS,j) # 个人认为更新误差应在更新b之后，因为公式算出的b的公式使用的是以前的EiupdateEk(oS,i)# b2^new=b2^old+(a1^old-a1^new)y1*K12+(a2^old-a2^new)y2*K22-E2if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]):oS.b = b1.A[0][0]elif (0<oS.alphas[j]) and (oS.C > oS.alphas[j]):oS.b = b2.A[0][0]else:oS.b = (b1+b2)/2.0print("b is {}".format(oS.b))return 1else:return 0def smoP(dataMatIn, classLabels, C,toler,maxIter,kTup=('lin',)):oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(),C,toler,kTup)iter = 0entireSet = True # 两种遍历方式交替alphaPairsChanged = 0while (iter<maxIter) and ((alphaPairsChanged>0) or (entireSet)):alphaPairsChanged = 0if entireSet:for i in range(oS.m):alphaPairsChanged += innerL(i,oS)print("fullSet, iter:%d i: %d pairs changed %d"%(iter,i ,alphaPairsChanged))iter+=1print("第一种遍历alphaRairChanged is {}".format(alphaPairsChanged))print("-----------eCache is {}".format(oS.eCache))print("***********alphas is {}".format(oS.alphas))print("---------------------------------------")else:nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] # 这时数组相乘，里面其实是True 和False的数组，得出来的是# 大于0并且小于C的alpha的下标for i in nonBoundIs:alphaPairsChanged += innerL(i,oS)print("non-bound, iter: %d i:%d, pairs changed %d"%(iter,i,alphaPairsChanged))print("第二种遍历alphaPairChanged is {}".format(alphaPairsChanged))iter+=1if entireSet:entireSet = False # 当第二种遍历方式alpha不再变化，那么继续第一种方式扫描，第一种方式不再变化，此时alphachanged为0且entireSet为false，退出循环elif (alphaPairsChanged==0):entireSet=Trueprint("iteration number: %d"%iter)return oS.b,oS.alphasdef calcWs(alphas,dataArr,classLabels): # 通过alpha来计算wX = mat(dataArr)labelMat = mat(classLabels).transpose()m,n = shape(X)w = zeros((n,1))for i in range(m):w += multiply(alphas[i]*labelMat[i], X[i,:].T) # w = sum(ai*yi*xi)return wdef draw_points(dataArr,classlabel, w,b,alphas):myfont = FontProperties(fname='/usr/share/fonts/simhei.ttf') # 显示中文plt.rcParams['axes.unicode_minus'] = False # 防止坐标轴的‘-’变为方块m = len(classlabel)red_points_x=[]red_points_y =[]blue_points_x=[]blue_points_y =[]svc_points_x =[]svc_points_y =[]# print(type(alphas))svc_point_index = nonzero((alphas.A>0) * (alphas.A <0.8))[0]svc_points = array(dataArr)[svc_point_index]svc_points_x = [x[0] for x in list(svc_points)]svc_points_y = [x[1] for x in list(svc_points)]print("svc_points_x",svc_points_x)print("svc_points_y",svc_points_y)for i in range(m):if classlabel[i] ==1:red_points_x.append(dataArr[i][0])red_points_y.append(dataArr[i][1])else:blue_points_x.append(dataArr[i][0])blue_points_y.append(dataArr[i][1])fig = plt.figure() # 创建画布ax = fig.add_subplot(111)ax.set_title("SVM-Classify") # 设置图片标题ax.set_xlabel("x") # 设置坐标名称ax.set_ylabel("y")ax1=ax.scatter(red_points_x, red_points_y, s=30,c='red', marker='s') #s是shape大小，c是颜色，marker是形状，'s'代表是正方形,默认'o'是圆圈ax2=ax.scatter(blue_points_x, blue_points_y, s=40,c='green')# ax.set_ylim([-6,5])print("b",b)print("w",w)x = arange(-4.0, 4.0, 0.1) # 分界线x范围,步长为0.1# x = arange(-2.0,10.0)if isinstance(b,numpy.matrixlib.defmatrix.matrix):b = b.A[0][0]y = (-b-w[0][0]*x)/w[1][0] # 直线方程 Ax + By + C = 0ax3,=plt.plot(x,y, 'k')ax4=plt.scatter(svc_points_x,svc_points_y,s=50,c='orange',marker='p')plt.legend([ax1, ax2,ax3,ax4], ["red points","blue points", "decision boundary","support vector"], loc='lower right') # 标注plt.show()dataArr,labelArr = loadDataSet('/home/zhangqingfeng/test/svm_test_data')
b,alphas = smoP(dataArr,labelArr,0.8,0.001,40)
w=calcWs(alphas,dataArr,labelArr)
draw_points(dataArr,labelArr,w,b,alphas)可参考数据集
-0.397822 8.058397 -1
0.824839 13.730343 -1
1.507278 5.027866 1
0.099671 6.835839 1
-0.344008 10.717485 -1
1.785928 7.718645 1
-0.918801 11.560217 -1
-0.364009 4.747300 1
-0.841722 4.119083 1
0.490426 1.960539 1
-0.007194 9.075792 -1
0.356107 12.447863 -1
0.342578 12.281162 -1
-0.810823 -1.466018 1
2.530777 6.476801 1
1.296683 11.607559 -1
0.475487 12.040035 -1
-0.783277 11.009725 -1
0.074798 11.023650 -1
-1.337472 0.468339 1
-0.102781 13.763651 -1
-0.147324 2.874846 1
0.518389 9.887035 -1
1.015399 7.571882 -1
-1.658086 -0.027255 1
1.319944 2.171228 1
2.056216 5.019981 1
-0.851633 4.375691 1
-1.510047 6.061992 -1
-1.076637 -3.181888 1
1.821096 10.283990 -1
3.010150 8.401766 1
-1.099458 1.688274 1
-0.834872 -1.733869 1
-0.846637 3.849075 1
1.400102 12.628781 -1
1.752842 5.468166 1
0.078557 0.059736 1
0.089392 -0.715300 1
1.825662 12.693808 -1
0.197445 9.744638 -1
0.126117 0.922311 1
-0.679797 1.220530 1
0.677983 2.556666 1
0.761349 10.693862 -1
-2.168791 0.143632 1
1.388610 9.341997 -1
0.317029 14.739025 -13.2 matlab
假设我们有一个简单的数据集，其中包含两个特征和一个二元标签。我们将使用MATLAB的fitcsvm函数来训练一个带径向基函数（RBF）核的SVM模型，并使用该模型对测试数据进行分类。
% 假设数据已经准备好，以下是模拟数据生成的代码
% 生成一些随机数据作为示例
rng(1); % 为了结果可重现
X = [randn(20,2)+1; randn(20,2)-1]; % 生成40个样本，每个样本2个特征
Y = [ones(20,1); -ones(20,1)]; % 生成标签% 将数据分为训练集和测试集
cv = cvpartition(size(X,1),'HoldOut',0.2);
idx = cv.test;% 训练数据
XTrain = X(~idx,:);
YTrain = Y(~idx);
% 测试数据
XTest = X(idx,:);
YTest = Y(idx);% 使用径向基核函数创建并训练SVM模型
SVMModel = fitcsvm(XTrain, YTrain, 'KernelFunction', 'rbf', 'BoxConstraint', 1);% 对测试数据进行预测
[label, score] = predict(SVMModel, XTest);% 计算准确率
accuracy = sum(label == YTest) / numel(YTest) * 100;% 显示结果
fprintf('准确率: %.2f%%\n', accuracy);% 可视化
figure;
gscatter(X(:,1), X(:,2), Y);
hold on;
sv = SVMModel.SupportVectors;
plot(sv(:,1),sv(:,2),'ko','MarkerSize',10);
title('SVM 分类结果与支持向量');
legend('类别 1', '类别 -1', '支持向量');
hold off;
需要参加数模竞赛的同学，可以看下面的名片，会有最新的助攻哦：（大型比赛前会对名片进行更新）
