机器学习实战4-SVM

1 SVM解释

    

  1. 通俗的讲,就是用来分类的。git

  2. 对于以上二维数据,支持向量就是那三个站在虚线上的点,中间那个粗黑硬实线表示的就是分割面(即超平面),能够将以上扩展到无数维,固然你确定想象不出100维时是个什么样子。算法

  3. 图上两虚线之间的距离称做间隔。数组

  4. 咱们的目标就是使间隔可以最大化,最大化间隔意味着能够更好的分类,即找到合适的超平面(就是那个黑粗硬实线),即求参数w和b。缓存

  5. 利用SMO算法能够高效解决4中问题数据结构

   阅读周志华《机器学习》app

2 SMO理解与简易实现

  关于SMO算法过程能够阅读本篇文章:https://zhuanlan.zhihu.com/p/29212107dom

  关于边界值C:能够理解为SVM对“软间隔的支持度”,C若是为无穷大,则全部的训练机器学习

  样本必须知足SVM的约束条件,C值越小,就容许越多的样本不知足约束条件。ide

  谈下几个主要过程:函数

  1. 初始化,须要数据集及其标签,常数C,容错率,最大循环次数
  2. while(确保α达到最优状态):获得优化后的α和b
    1. for(从1一直循环到m):
      1. 固定αi,计算fxi和ei
      2. 判断αi是否已在边界上,在就不能继续优化
      3. 计算fxj,ej
      4. 计算αj的上下界L,H等相关的值
      5. 更新αj
      6. 更新αi
      7. 更新阈值b
  3. 返回w和b

  看上去好像更糊涂=。=,仍是看上面的算法文章靠谱,如下是代码的简易实现:

 1 from numpy import *
 2 import matplotlib.pyplot as plt  3 
 4 # 导入文件数据
 5 def loadDataSet(fileName):  6     dataMat=[]  7     labelMat=[]  8     fr=open(fileName)  9     for line in fr.readlines():  10         lineArr=line.strip().split('\t')  11         dataMat.append([float(lineArr[0]),float(lineArr[1])])  12         labelMat.append(float(lineArr[2]))  13     return dataMat,labelMat  14 
 15 # m为alpha个数,i为下标,随机输出与i不一样的j
 16 def selectJrand(i,m):  17     j=i  18     while(j==i):  19         j=int(random.uniform(0,m))  20     return j  21 
 22 # 用于调整大于H、或者小于L的alpha
 23 def clipAlpha(aj,H,L):  24     if aj>H:  25         aj=H  26     if L>aj:  27         aj=L  28     return aj  29 
 30 # 简化版本SMO算法
 31 def smoSimple(dataMatIn,classLabels,C,toler,maxIter):  32     # 初始化条件,数据集,类别标签,常数C,容错率和取消前最大的循环次数
 33     dataMatrix=mat(dataMatIn)  34     labelMat=mat(classLabels).transpose()  35     b=0;m,n=shape(dataMatrix)  36     alphas=mat(zeros((m,1)))  37     iter=0  38 
 39     while(iter<maxIter):  40         alphaPairsChanged=0  # alphapaitchanged 是否已进行优化
 41         for i in range(m):  42             # 与周志华机器学习公式有一丝区别:x为m*n矩阵 周志华中xm为列向量,此处为行向量
 43             # 即公式为f(xi) = ∑alpha_mi*y_mi*x_mi*xi.T + b 可理解为距离 f(x)为m维行向量
 44             fXi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T))+b # multiply为对应元素相乘
 45             Ei=fXi-float(labelMat[i]) # 计算误差
 46             # 判断条件:正负间隔均判断,同时检查alpha值,保证其不等于C或者0
 47             if ((labelMat[i]*Ei<-toler)and(alphas[i]<C)) or ((labelMat[i]*Ei>toler)and(alphas[i]>0)):  48                 j=selectJrand(i,m) #取不一样于i的j
 49                 fXj=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b  50                 Ej = fXj - float(labelMat[j])  51                 alphaIold=alphas[i].copy()  52                 alphaJold=alphas[j].copy()  53                 # 计算L、H,即alpha_j的上下界
 54                 if(labelMat[i]!=labelMat[j]):  55                     L=max(0,alphas[j]-alphas[i])  56                     H=min(C,C+alphas[j]-alphas[i])  57                 else:  58                     L=max(0,alphas[j]+alphas[i]-C)  59                     H=min(C,alphas[j]+alphas[i])  60                 if L==H:  61                     # print("L=H")
 62                     continue
 63                 # η=-(K11+K22−2K12)也可去掉负号,但相应更新alpha_j时累减改为累加
 64                 eta=2.0*dataMatrix[i,:]*dataMatrix[j,:].T-\  65                     dataMatrix[i,:]*dataMatrix[i,:].T-\  66                     dataMatrix[j,:]*dataMatrix[j,:].T  67                 if eta>=0:  68                     # print("eta>=0")
 69                     continue
 70                 # 更新alpha_j
 71                 alphas[j]-=labelMat[j]*(Ei-Ej)/eta  72                 alphas[j]=clipAlpha(alphas[j],H,L)  73                 if(abs(alphas[j]-alphaJold)<0.00001):  74                     # print("j not moving enough")
 75                     continue
 76                 # 更新alpha_i
 77                 alphas[i]+=labelMat[j]*labelMat[i]*(alphaJold-alphas[j])  78                 # 更新阈值b
 79                 b1=b-Ei-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T-\  80                     labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T  81                 b2=b-Ej-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T-\  82                     labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T  83                 if (0<alphas[i])and(C>alphas[i]):  84                     b=b1  85                 elif (0<alphas[j])and(C>alphas[j]):  86                     b=b2  87                 else:  88                     b=(b1+b2)/2.0
 89                 alphaPairsChanged+=1
 90                 # print("iter:%d i:%d,pairs changed %d"%(iter,i,alphaPairsChanged))
 91         if(alphaPairsChanged==0):  92             iter+=1
 93         else:  94             iter=0  95         # print("iteration number: %d"%iter)
 96     return b,alphas  97 
 98 # 分类数据点
 99 def classPts(dataArr,labelArr): 100     classifiedPts = {'+1': [], '-1': []} 101     # zip 进行封装,返回的是一个对象,可用list将其转换为列表查看
102     for point, label in zip(dataArr, labelArr): 103         if label == 1.0: 104             classifiedPts['+1'].append(point) 105         else: 106             classifiedPts['-1'].append(point) 107     return classifiedPts 108 
109 # 经过已知数据点和拉格朗日乘子得到分割超平面参数w
110 def get_w(alphas, dataset, labels): 111     import numpy as np 112     alphas, dataset, labels = np.array(alphas), np.array(dataset), np.array(labels) 113     yx = labels.reshape(1, -1).T*np.array([1, 1])*dataset 114     w = np.dot(yx.T, alphas) 115     return w.tolist() 116 
117 # 绘制数据点,分割线,支持向量
118 def drawSvm(classPts,alphas,dataArr,labelArr): 119     # 绘制数据点
120     import numpy as ny 121     for label, pts in classPts.items(): 122         pts = ny.array(pts) 123         ax.scatter(pts[:, 0], pts[:, 1], label=label) 124     # 绘制分割线
125     w = get_w(alphas, dataArr, labelArr) 126     x1, _ = max(dataArr, key=lambda x: x[0]) 127     x2, _ = min(dataArr, key=lambda x: x[0]) 128     a1, a2 = w 129     y1, y2 = (-float(b) - a1[0]*x1)/a2[0], (-float(b) -a1[0]*x2)/a2[0] # b是矩阵类型,a1是列表类型,x1是float,注意不一样形式之间的转换
130  ax.plot([x1, x2], [y1, y2]) 131     # 绘制支持向量
132     for i, alpha in enumerate(alphas): 133         if abs(alpha) > 1e-3: 134             x, y = dataArr[i] 135             ax.scatter([x], [y], s=150, c='none', alpha=0.7, 136                        linewidth=1.5, edgecolor='#AB3319') 137  plt.show() 138 
139 
140 
141 dataArr,labelArr=loadDataSet('testSet.txt') 142 b,alphas=smoSimple(dataArr,labelArr,0.6,0.001,40) 143 classFieldPts=classPts(dataArr,labelArr) 144 ax = plt.figure().add_subplot(111) 145 drawSvm(classFieldPts,alphas,dataArr,labelArr)

  输出结果:

3 完整版Platt SMO算法

  目的:为了进一步优化算法,起到加速做用

两算法区别
  简版SMO        完整Platt SMO
while循环

退出条件:

迭代次数超过指定最大次数

 

退出条件:

1.迭代次数超过指定最大次数

2.遍历整个集合都未对任意alpha对进行修改

iter:

当没有任何alpha值发生变化

时,计一次迭代,即iter+1

iter:

做为一次循环过程,即每一次循环iter+1

 

while循环内

 选择j后会计算错误率Ej

(右边的算法则使用一个全局

缓存保存偏差值,选择Ei-Ej

最大的alpha值)

1.遍历任意可能的alpha

经过inneL选择第二个alpha,可能时进行优化处理

2.遍历全部非边界(即不在边界0或者C上)alpha值

在以上两种方式中进行切换

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

完整代码:

 1 from numpy import *
 2 import matplotlib.pyplot as plt  3 
 4 # 导入文件数据
 5 def loadDataSet(fileName):  6     dataMat=[]  7     labelMat=[]  8     fr=open(fileName)  9     for line in fr.readlines():  10         lineArr=line.strip().split('\t')  11         dataMat.append([float(lineArr[0]),float(lineArr[1])])  12         labelMat.append(float(lineArr[2]))  13     return dataMat,labelMat  14 
 15 # m为alpha个数,i为下标,随机输出与i不一样的j
 16 def selectJrand(i,m):  17     j=i  18     while(j==i):  19         j=int(random.uniform(0,m))  20     return j  21 
 22 # 用于调整大于H、或者小于L的alpha
 23 def clipAlpha(aj,H,L):  24     if aj>H:  25         aj=H  26     if L>aj:  27         aj=L  28     return aj  29 
 30 # 简化版本SMO算法
 31 def smoSimple(dataMatIn,classLabels,C,toler,maxIter):  32     # 初始化条件,数据集,类别标签,常数C,容错率和取消前最大的循环次数
 33     dataMatrix=mat(dataMatIn)  34     labelMat=mat(classLabels).transpose()  35     b=0;m,n=shape(dataMatrix)  36     alphas=mat(zeros((m,1)))  37     iter=0  38 
 39     while(iter<maxIter):  40         alphaPairsChanged=0  # alphapaitchanged 是否已进行优化
 41         for i in range(m):  42             # 与周志华机器学习公式有一丝区别:x为m*n矩阵 周志华中xm为列向量,此处为行向量
 43             # 即公式为f(xi) = ∑alpha_mi*y_mi*x_mi*xi.T + b 可理解为距离 f(x)为m维行向量
 44             fXi=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T))+b # multiply为对应元素相乘
 45             Ei=fXi-float(labelMat[i]) # 计算误差
 46             # 判断条件:正负间隔均判断,同时检查alpha值,保证其不等于C或者0
 47             if ((labelMat[i]*Ei<-toler)and(alphas[i]<C)) or ((labelMat[i]*Ei>toler)and(alphas[i]>0)):  48                 j=selectJrand(i,m) #取不一样于i的j
 49                 fXj=float(multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T))+b  50                 Ej = fXj - float(labelMat[j])  51                 alphaIold = alphas[i].copy()  52                 alphaJold=alphas[j].copy()  53                 # 计算L、H,即alpha_j的上下界
 54                 if(labelMat[i]!=labelMat[j]):  55                     L = max(0, alphas[j]-alphas[i])  56                     H = min(C, C+alphas[j]-alphas[i])  57                 else:  58                     L=max(0, alphas[j]+alphas[i]-C)  59                     H=min(C, alphas[j]+alphas[i])  60                 if L==H:  61                     # print("L=H")
 62                     continue
 63                 # η=-(K11+K22−2K12)也可去掉负号,但相应更新alpha_j时累减改为累加
 64                 eta=2.0*dataMatrix[i,:]*dataMatrix[j,:].T-\  65                     dataMatrix[i,:]*dataMatrix[i,:].T-\  66                     dataMatrix[j,:]*dataMatrix[j,:].T  67                 if eta>=0:  68                     # print("eta>=0")
 69                     continue
 70                 # 更新alpha_j
 71                 alphas[j]-=labelMat[j]*(Ei-Ej)/eta  72                 alphas[j]=clipAlpha(alphas[j],H,L)  73                 if(abs(alphas[j]-alphaJold)<0.00001):  74                     # print("j not moving enough")
 75                     continue
 76                 # 更新alpha_i
 77                 alphas[i]+=labelMat[j]*labelMat[i]*(alphaJold-alphas[j])  78                 # 更新阈值b
 79                 b1=b-Ei-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T-\  80                     labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T  81                 b2=b-Ej-labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T-\  82                     labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T  83                 if (0<alphas[i])and(C>alphas[i]):  84                     b=b1  85                 elif (0<alphas[j])and(C>alphas[j]):  86                     b=b2  87                 else:  88                     b=(b1+b2)/2.0
 89                 alphaPairsChanged+=1
 90                 # print("iter:%d i:%d,pairs changed %d"%(iter,i,alphaPairsChanged))
 91         if(alphaPairsChanged==0):  92             iter+=1
 93         else:  94             iter=0  95         # print("iteration number: %d"%iter)
 96     return b,alphas  97 
 98 # 1. 构建一个数据结构存储全部数据
 99 class optStruct: 100     def __init__(self, dataMatIn, classLabels, C, toler): 101         self.X = dataMatIn 102         self.labelMat = classLabels 103         self.C = C 104         self.tol = toler 105         self.m = shape(dataMatIn)[0] 106         self.alphas = mat(zeros((self.m,1))) 107         self.b = 0 108         self.eCache = mat(zeros((self.m,2))) # 偏差缓存 第一列是eCache是否有效的标志位,第二列是实际E值
109 
110 # 2. 计算k处偏差值
111 def calcEk(oS,k): 112     fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T))+oS.b 113     Ek = fXk-float(oS.labelMat[k]) 114     return Ek 115 
116 # 3. 用于选择第二个α
117 def selectJ(i,oS,Ei): 118     maxK = -1 ; maxDeltaE = 0 ; Ej = 0 119     oS.eCache[i]=[1,Ei] 120     validEcacheList = nonzero(oS.eCache[:,0].A)[0]   # nonzero 返回数组中非零元素的索引值
121     if (len(validEcacheList)) > 1: 122         # 遍历一遍上列表,找出最大差值的Ej
123         for k in validEcacheList: 124             if k == i: 125                 continue
126             Ek = calcEk(oS,k) 127             deltaE = abs(Ei-Ek) 128             if (deltaE>maxDeltaE): 129                 maxK=k 130                 maxDeltaE=deltaE 131                 Ej=Ek 132         return maxK,Ej 133     else: 134         # 第一次,直接随机寻找一个Ej
135         j=selectJrand(i,oS.m) 136         Ej=calcEk(oS,j) 137     return j,Ej 138 
139 # 4. 计算k处偏差值 并存入eCache中
140 def updateEk(oS,k): 141     Ek=calcEk(oS,k) 142     oS.eCache[k]=[1,Ek] 143 
144 # 5. 内部循环,找到配对Ej,则返回1,
145 def innerL(i,oS): 146     Ei=calcEk(oS,i) 147     if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \ 148        ((oS.labelMat[i] * Ei >  oS.tol) and (oS.alphas[i] > 0)): 149         j,Ej=selectJ(i,oS,Ei) 150         alphaIold=oS.alphas[i].copy() 151         alphaJold=oS.alphas[j].copy() 152         if (oS.labelMat[i] != oS.labelMat[j]): 153             L = max(0, oS.alphas[j] - oS.alphas[i]) 154             H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 155         else: 156             L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 157             H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 158         if L == H: 159             print("L=H") 160             return 0 161         eta=2.0*oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T-\ 162             oS.X[j,:]*oS.X[j,:].T 163         if eta >= 0: 164             print("eta>=0") 165             return 0 166         oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta 167         oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) 168  updateEk(oS,j) 169         if (abs(oS.alphas[j] - alphaJold) < 0.00001): 170             print("j not moving enough") 171             return 0 172         oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j]) 173  updateEk(oS, i) 174         b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[i, :].T - oS.labelMat[j] * ( 175                     oS.alphas[j] - alphaJold) * oS.X[i, :] * oS.X[j, :].T 176         b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.X[i, :] * oS.X[j, :].T - oS.labelMat[j] * ( 177                     oS.alphas[j] - alphaJold) * oS.X[j, :] * oS.X[j, :].T 178         if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): 179             oS.b = b1 180         elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): 181             oS.b = b2 182         else: 183             oS.b = (b1 + b2) / 2.0
184         return 1
185     else: 186         return 0 187 
188 # 6. 完整函数
189 def smoP(dataMatIn, classLabels, C, toler, maxIter): 190     # 建立一个 optStruct 对象
191     oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler) 192     iter = 0 193     entireSet = True 194     alphaPairsChanged = 0 195     # 循环遍历:循环maxIter次 而且 (alphaPairsChanged存在能够改变 or 全部行遍历一遍)
196     # 循环迭代结束 或者 循环遍历全部alpha后,alphaPairs仍是没变化
197     while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): 198         alphaPairsChanged = 0 199         # 当entireSet=true or 非边界alpha对没有了;就开始寻找 alpha对,而后决定是否要进行else。
200         if entireSet: 201             # 在数据集上遍历全部可能的alpha
202             for i in range(oS.m): 203                 # 是否存在alpha对,存在就+1
204                 alphaPairsChanged += innerL(i, oS) 205                 print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 206             iter += 1
207         # 对已存在 alpha对,选出非边界的alpha值,进行优化。
208         else: 209             # 遍历全部的非边界alpha值,也就是不在边界0或C上的值。
210             nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 211             for i in nonBoundIs: 212                 alphaPairsChanged += innerL(i, oS) 213                 print("non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 214             iter += 1
215         # 若是找到alpha对,就优化非边界alpha值,不然,就从新进行寻找,若是寻找一遍 遍历全部的行仍是没找到,就退出循环。
216         if entireSet: 217             entireSet = False  # toggle entire set loop
218         elif (alphaPairsChanged == 0): 219             entireSet = True 220         print("iteration number: %d" % iter) 221     return oS.b, oS.alphas 222 
223 # 7. 经过已知数据点和拉格朗日乘子得到分割超平面参数w
224 def calcWs(alphas, dataArr, classLabels): 225     X = mat(dataArr) 226     labelMat = mat(classLabels).transpose() 227     m, n = shape(X) 228     w = zeros((n, 1)) 229     for i in range(m): 230         w += multiply(alphas[i] * labelMat[i], X[i, :].T) 231     return w 232 
233 # 8. 绘制数据点,分割线,支持向量
234 def drawSvm(alphas,dataArr,labelArr): 235     # 对数据点进行分类
236     classifiedPts = {'+1': [], '-1': []} 237     # zip 进行封装,返回的是一个对象,可用list将其转换为列表查看
238     for point, label in zip(dataArr, labelArr): 239         if label == 1.0: 240             classifiedPts['+1'].append(point) 241         else: 242             classifiedPts['-1'].append(point) 243     ax = plt.figure().add_subplot(111) 244     # 绘制数据点
245     import numpy as ny 246     for label, pts in classifiedPts.items(): 247         pts = ny.array(pts) 248         ax.scatter(pts[:, 0], pts[:, 1], label=label) 249     # 绘制分割线
250     w = calcWs(alphas, dataArr, labelArr) 251     x1, _ = max(dataArr, key=lambda x: x[0]) 252     x2, _ = min(dataArr, key=lambda x: x[0]) 253     a1, a2 = w 254     y1, y2 = (-float(b) - a1[0]*x1)/a2[0], (-float(b) -a1[0]*x2)/a2[0] # b是矩阵类型,a1是列表类型,x1是float,注意不一样形式之间的转换
255  ax.plot([x1, x2], [y1, y2]) 256     # 绘制支持向量
257     for i, alpha in enumerate(alphas): 258         if abs(alpha) > 1e-3: 259             x, y = dataArr[i] 260             ax.scatter([x], [y], s=150, c='none', alpha=0.7, 261                        linewidth=1.5, edgecolor='#AB3319') 262  plt.show() 263 
264 dataArr,labelArr=loadDataSet('testSet.txt') 265 b,alphas=smoP(dataArr,labelArr,0.9,0.001,40) 266 drawSvm(alphas,dataArr,labelArr)
View Code

输出结果:

4 利用核函数将数据映射到高维空间

能够将原始空间映射到一个更高维的特征空间,使得样本在这个高维空间中能够线性可分。以下图:

经常使用核函数以下:

核转换函数程序以下:

 1 # 核转换函数
 2 def kernelTrans(X,A,kTup):  3     m,n=shape(X)  4     K=mat(zeros((m,1)))  5     if kTup[0]=='lin':  6         # 线性核函数
 7         K=X*A.T  8     elif kTup[0]=='rbf':  9         # 高斯核函数
10         for j in range(m): 11             deltaRow =X[j,:]-A 12             K[j]=deltaRow*deltaRow.T 13         K=exp(K/(-1*kTup[1]**2)) 14     else: 15         raise NameError('houston we have a problem that kenel is not recognized') 16     return K

使用核函数进行svm分类总程序及结果以下:

 1 from numpy import *
 2 import matplotlib.pyplot as plt  3 
 4 # 导入文件数据
 5 def loadDataSet(fileName):  6     dataMat=[]  7     labelMat=[]  8     fr=open(fileName)  9     for line in fr.readlines():  10         lineArr=line.strip().split('\t')  11         dataMat.append([float(lineArr[0]),float(lineArr[1])])  12         labelMat.append(float(lineArr[2]))  13     return dataMat,labelMat  14 
 15 # m为alpha个数,i为下标,随机输出与i不一样的j
 16 def selectJrand(i,m):  17     j=i  18     while(j==i):  19         j=int(random.uniform(0,m))  20     return j  21 
 22 # 用于调整大于H、或者小于L的alpha
 23 def clipAlpha(aj,H,L):  24     if aj>H:  25         aj=H  26     if L>aj:  27         aj=L  28     return aj  29 
 30 # 核转换函数
 31 def kernelTrans(X,A,kTup):  32     m,n=shape(X)  33     K=mat(zeros((m,1)))  34     if kTup[0]=='lin':  35         # 线性核函数
 36         K=X*A.T  37     elif kTup[0]=='rbf':  38         # 高斯核函数
 39         for j in range(m):  40             deltaRow =X[j,:]-A  41             K[j]=deltaRow*deltaRow.T  42         K=exp(K/(-1*kTup[1]**2))  43     else:  44         raise NameError('houston we have a problem that kenel is not recognized')  45     return K  46 
 47 # 1.加核新结构
 48 class optStruct:  49     def __init__(self, dataMatIn, classLabels, C, toler,kTup):  50         self.X = dataMatIn  51         self.labelMat = classLabels  52         self.C = C  53         self.tol = toler  54         self.m = shape(dataMatIn)[0]  55         self.alphas = mat(zeros((self.m,1)))  56         self.b = 0  57         self.eCache = mat(zeros((self.m,2))) # 偏差缓存 第一列是eCache是否有效的标志位,第二列是实际E值
 58         self.K=mat(zeros((self.m,self.m)))  59         # 计算K
 60         for i in range(self.m):  61             self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup)  62 
 63 # 2.加核新calcEk
 64 def calcEk(oS,k):  65     fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k])+oS.b  66     Ek = fXk-float(oS.labelMat[k])  67     return Ek  68 
 69 # 3. 用于选择第二个α
 70 def selectJ(i,oS,Ei):  71     maxK = -1 ; maxDeltaE = 0 ; Ej = 0  72     oS.eCache[i]=[1,Ei]  73     validEcacheList = nonzero(oS.eCache[:,0].A)[0]   # nonzero 返回数组中非零元素的索引值
 74     if (len(validEcacheList)) > 1:  75         # 遍历一遍上列表,找出最大差值的Ej
 76         for k in validEcacheList:  77             if k == i:  78                 continue
 79             Ek = calcEk(oS,k)  80             deltaE = abs(Ei-Ek)  81             if (deltaE>maxDeltaE):  82                 maxK=k  83                 maxDeltaE=deltaE  84                 Ej=Ek  85         return maxK,Ej  86     else:  87         # 第一次,直接随机寻找一个Ej
 88         j=selectJrand(i,oS.m)  89         Ej=calcEk(oS,j)  90     return j,Ej  91 
 92 # 4. 计算k处偏差值 并存入eCache中
 93 def updateEk(oS,k):  94     Ek=calcEk(oS,k)  95     oS.eCache[k]=[1,Ek]  96 
 97 # 5.加核新innerL
 98 def innerL(i,oS):  99     Ei=calcEk(oS,i) 100     if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \ 101        ((oS.labelMat[i] * Ei >  oS.tol) and (oS.alphas[i] > 0)): 102         j,Ej=selectJ(i,oS,Ei) 103         alphaIold=oS.alphas[i].copy() 104         alphaJold=oS.alphas[j].copy() 105         if (oS.labelMat[i] != oS.labelMat[j]): 106             L = max(0, oS.alphas[j] - oS.alphas[i]) 107             H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 108         else: 109             L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 110             H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 111         if L == H: 112             print("L=H") 113             return 0 114         eta=2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j] 115         if eta >= 0: 116             print("eta>=0") 117             return 0 118         oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta 119         oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) 120  updateEk(oS,j) 121         if (abs(oS.alphas[j] - alphaJold) < 0.00001): 122             print("j not moving enough") 123             return 0 124         oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j]) 125  updateEk(oS, i) 126         b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - \ 127                          oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.K[i,j] 128         b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - \ 129                          oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.K[j,j] 130         if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): 131             oS.b = b1 132         elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): 133             oS.b = b2 134         else: 135             oS.b = (b1 + b2) / 2.0
136         return 1
137     else: 138         return 0 139 
140 # 6. 完整函数
141 def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup): 142     # 建立一个 optStruct 对象
143     oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler,kTup) 144     iter = 0 145     entireSet = True 146     alphaPairsChanged = 0 147     # 循环遍历:循环maxIter次 而且 (alphaPairsChanged存在能够改变 or 全部行遍历一遍)
148     # 循环迭代结束 或者 循环遍历全部alpha后,alphaPairs仍是没变化
149     while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): 150         alphaPairsChanged = 0 151         # 当entireSet=true or 非边界alpha对没有了;就开始寻找 alpha对,而后决定是否要进行else。
152         if entireSet: 153             # 在数据集上遍历全部可能的alpha
154             for i in range(oS.m): 155                 # 是否存在alpha对,存在就+1
156                 alphaPairsChanged += innerL(i, oS) 157                 print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 158             iter += 1
159         # 对已存在 alpha对,选出非边界的alpha值,进行优化。
160         else: 161             # 遍历全部的非边界alpha值,也就是不在边界0或C上的值。
162             nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 163             for i in nonBoundIs: 164                 alphaPairsChanged += innerL(i, oS) 165                 print("non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 166             iter += 1
167         # 若是找到alpha对,就优化非边界alpha值,不然,就从新进行寻找,若是寻找一遍 遍历全部的行仍是没找到,就退出循环。
168         if entireSet: 169             entireSet = False  # toggle entire set loop
170         elif (alphaPairsChanged == 0): 171             entireSet = True 172         print("iteration number: %d" % iter) 173     return oS.b, oS.alphas 174 
175 # 7. 经过已知数据点和拉格朗日乘子得到分割超平面参数w
176 def calcWs(alphas, dataArr, classLabels): 177     X = mat(dataArr) 178     labelMat = mat(classLabels).transpose() 179     m, n = shape(X) 180     w = zeros((n, 1)) 181     for i in range(m): 182         w += multiply(alphas[i] * labelMat[i], X[i, :].T) 183     return w 184 
185 # 8. 绘制数据点,分割线,支持向量
186 def drawSvm(alphas,dataArr,labelArr): 187     # 对数据点进行分类
188     classifiedPts = {'+1': [], '-1': []} 189     # zip 进行封装,返回的是一个对象,可用list将其转换为列表查看
190     for point, label in zip(dataArr, labelArr): 191         if label == 1.0: 192             classifiedPts['+1'].append(point) 193         else: 194             classifiedPts['-1'].append(point) 195     ax = plt.figure().add_subplot(111) 196     # 绘制数据点
197     import numpy as ny 198     for label, pts in classifiedPts.items(): 199         pts = ny.array(pts) 200         ax.scatter(pts[:, 0], pts[:, 1], label=label) 201     # 绘制分割线
202     w = calcWs(alphas, dataArr, labelArr) 203     x1, _ = max(dataArr, key=lambda x: x[0]) 204     x2, _ = min(dataArr, key=lambda x: x[0]) 205     a1, a2 = w 206     y1, y2 = (-float(b) - a1[0]*x1)/a2[0], (-float(b) -a1[0]*x2)/a2[0] # b是矩阵类型,a1是列表类型,x1是float,注意不一样形式之间的转换
207  ax.plot([x1, x2], [y1, y2]) 208     # 绘制支持向量
209     for i, alpha in enumerate(alphas): 210         if abs(alpha) > 1e-3: 211             x, y = dataArr[i] 212             ax.scatter([x], [y], s=150, c='none', alpha=0.7, 213                        linewidth=1.5, edgecolor='#AB3319') 214  plt.show() 215 
216 # 9.测试函数
217 def testRbf(k1=1.5): 218     dataArr, labelArr = loadDataSet('testSetRBF.txt') 219     b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000,('rbf',k1)) 220     datMat=mat(dataArr) 221     labelMat=mat(labelArr).transpose() 222     svInd=nonzero(alphas.A>0)[0] 223     sVs=datMat[svInd] 224     labelSV=labelMat[svInd] 225     print("there are %d Support Vectors" %shape(sVs)[0]) 226     # 如何利用核函数进行分类
227     m,n=shape(datMat) 228     errorCount=0 229     for i in range(m): 230         kenelEval=kernelTrans(sVs,datMat[i,:],('rbf',k1)) 231         predict=kenelEval.T*multiply(labelSV,alphas[svInd])+b 232         if sign(predict)!=sign(labelArr[i]): 233             errorCount+=1
234     print("the training error rate is %f"%(float(errorCount)/m)) 235     dataArr,labelArr=loadDataSet('testSetRBF2.txt') 236     errorCount=0 237     datMat=mat(dataArr);labelMat=mat(labelArr).transpose() 238     m,n=shape(datMat) 239     for i in range(m): 240         kenelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1)) 241     predict = kenelEval.T * multiply(labelSV, alphas[svInd] )+ b 242     if sign(predict) != sign(labelArr[i]): 243         errorCount += 1
244     print("the training error rate is %f" % (float(errorCount) / m)) 245 
246 testRbf()
View Code

5 实例:手写数据集

能够直接采用前面的函数,不过输入数集须要稍微进行处理。本例有两个文件夹,一个存放训练集,一个存放测试集。每一个文件夹中一个文件表明一个一张手写图片,是一个32*32向量,咱们须要利用一个函数将其转化为1*1024行向量。将文件夹中全部文件放入矩阵中,这样能够获得训练集矩阵和标签,一样获得测试集矩阵。程序以下:

 1 # 将32x32的二进制图像转换为1x1024向量。
 2 def img2vector(filename):  3     import numpy as np  4     returnVect = np.zeros((1,1024))  5     fr = open(filename)  6     for i in range(32):  7         lineStr = fr.readline()  8         for j in range(32):  9             returnVect[0,32*i+j] = int(lineStr[j]) 10     return returnVect 11 
12 # 导入数据集和训练集函数
13 def loadImages(dirName): 14     from os import listdir 15     hwLabels = [] 16     trainingFileList = listdir(dirName)  # listdir 返回指定文件夹中全部文件及文件夹的名称列表
17     m = len(trainingFileList) 18     trainingMat = zeros((m, 1024))  # 存放训练数据矩阵
19     for i in range(m): 20         fileNameStr = trainingFileList[i] 21         fileStr = fileNameStr.split('.')[0] 22         classNumStr = int(fileStr.split('_')[0]) 23         # hwlabels存放数集标签,1的标签为1,9的标枪为-1
24         if classNumStr == 9:  # 若是为9,输出-1
25             hwLabels.append(-1) 26         else:  # 若是不为9,输出1
27             hwLabels.append(1) 28         trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))  # 调用img2vector()函数
29     return trainingMat, hwLabels

识别总代码

 1 # 手写数字集svm识别 1和9
 2 import random  3 from numpy import *
 4 
 5 # 将32x32的二进制图像转换为1x1024向量。
 6 def img2vector(filename):  7     import numpy as np  8     returnVect = np.zeros((1,1024))  9     fr = open(filename)  10     for i in range(32):  11         lineStr = fr.readline()  12         for j in range(32):  13             returnVect[0,32*i+j] = int(lineStr[j])  14     return returnVect  15 
 16 # 导入数据集和训练集函数
 17 def loadImages(dirName):  18     from os import listdir  19     hwLabels = []  20     trainingFileList = listdir(dirName)  # listdir 返回指定文件夹中全部文件及文件夹的名称列表
 21     m = len(trainingFileList)  22     trainingMat = zeros((m, 1024))  # 存放训练数据矩阵
 23     for i in range(m):  24         fileNameStr = trainingFileList[i]  25         fileStr = fileNameStr.split('.')[0]  26         classNumStr = int(fileStr.split('_')[0])  27         # hwlabels存放数集标签,1的标签为1,9的标枪为-1
 28         if classNumStr == 9:  # 若是为9,输出-1
 29             hwLabels.append(-1)  30         else:  # 若是不为9,输出1
 31             hwLabels.append(1)  32         trainingMat[i, :] = img2vector('%s/%s' % (dirName, fileNameStr))  # 调用img2vector()函数
 33     return trainingMat, hwLabels  34 
 35 # m为alpha个数,i为下标,随机输出与i不一样的j
 36 def selectJrand(i,m):  37     j=i  38     while(j==i):  39         j=int(random.uniform(0,m))  40     return j  41 
 42 # 用于调整大于H、或者小于L的alpha
 43 def clipAlpha(aj,H,L):  44     if aj>H:  45         aj=H  46     if L>aj:  47         aj=L  48     return aj  49 
 50 # 核转换函数
 51 def kernelTrans(X,A,kTup):  52     m,n=shape(X)  53     K=mat(zeros((m,1)))  54     if kTup[0]=='lin':  55         # 线性核函数
 56         K=X*A.T  57     elif kTup[0]=='rbf':  58         # 高斯核函数
 59         for j in range(m):  60             deltaRow =X[j,:]-A  61             K[j]=deltaRow*deltaRow.T  62         K=exp(K/(-1*kTup[1]**2))  63     else:  64         raise NameError('houston we have a problem that kenel is not recognized')  65     return K  66 
 67 # 1.加核新结构
 68 class optStruct:  69     def __init__(self, dataMatIn, classLabels, C, toler,kTup):  70         self.X = dataMatIn  71         self.labelMat = classLabels  72         self.C = C  73         self.tol = toler  74         self.m = shape(dataMatIn)[0]  75         self.alphas = mat(zeros((self.m,1)))  76         self.b = 0  77         self.eCache = mat(zeros((self.m,2))) # 偏差缓存 第一列是eCache是否有效的标志位,第二列是实际E值
 78         self.K=mat(zeros((self.m,self.m)))  79         # 计算K
 80         for i in range(self.m):  81             self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup)  82 
 83 # 2.加核新calcEk
 84 def calcEk(oS,k):  85     fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k])+oS.b  86     Ek = fXk-float(oS.labelMat[k])  87     return Ek  88 
 89 # 3. 用于选择第二个α
 90 def selectJ(i,oS,Ei):  91     maxK = -1 ; maxDeltaE = 0 ; Ej = 0  92     oS.eCache[i]=[1,Ei]  93     validEcacheList = nonzero(oS.eCache[:,0].A)[0]   # nonzero 返回数组中非零元素的索引值
 94     if (len(validEcacheList)) > 1:  95         # 遍历一遍上列表,找出最大差值的Ej
 96         for k in validEcacheList:  97             if k == i:  98                 continue
 99             Ek = calcEk(oS,k) 100             deltaE = abs(Ei-Ek) 101             if (deltaE>maxDeltaE): 102                 maxK=k 103                 maxDeltaE=deltaE 104                 Ej=Ek 105         return maxK,Ej 106     else: 107         # 第一次,直接随机寻找一个Ej
108         j=selectJrand(i,oS.m) 109         Ej=calcEk(oS,j) 110     return j,Ej 111 
112 # 4. 计算k处偏差值 并存入eCache中
113 def updateEk(oS,k): 114     Ek=calcEk(oS,k) 115     oS.eCache[k]=[1,Ek] 116 
117 # 5.加核新innerL
118 def innerL(i,oS): 119     Ei=calcEk(oS,i) 120     if ((oS.labelMat[i] * Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or \ 121        ((oS.labelMat[i] * Ei >  oS.tol) and (oS.alphas[i] > 0)): 122         j,Ej=selectJ(i,oS,Ei) 123         alphaIold=oS.alphas[i].copy() 124         alphaJold=oS.alphas[j].copy() 125         if (oS.labelMat[i] != oS.labelMat[j]): 126             L = max(0, oS.alphas[j] - oS.alphas[i]) 127             H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) 128         else: 129             L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) 130             H = min(oS.C, oS.alphas[j] + oS.alphas[i]) 131         if L == H: 132             print("L=H") 133             return 0 134         eta=2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j] 135         if eta >= 0: 136             print("eta>=0") 137             return 0 138         oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta 139         oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) 140  updateEk(oS,j) 141         if (abs(oS.alphas[j] - alphaJold) < 0.00001): 142             print("j not moving enough") 143             return 0 144         oS.alphas[i] += oS.labelMat[j] * oS.labelMat[i] * (alphaJold - oS.alphas[j]) 145  updateEk(oS, i) 146         b1 = oS.b - Ei - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i,i] - \ 147                          oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.K[i,j] 148         b2 = oS.b - Ej - oS.labelMat[i] * (oS.alphas[i] - alphaIold) * oS.K[i,j] - \ 149                          oS.labelMat[j] * (oS.alphas[j] - alphaJold) * oS.K[j,j] 150         if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): 151             oS.b = b1 152         elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): 153             oS.b = b2 154         else: 155             oS.b = (b1 + b2) / 2.0
156         return 1
157     else: 158         return 0 159 
160 # 6. 完整函数
161 def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup): 162     # 建立一个 optStruct 对象
163     oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(), C, toler,kTup) 164     iter = 0 165     entireSet = True 166     alphaPairsChanged = 0 167     # 循环遍历:循环maxIter次 而且 (alphaPairsChanged存在能够改变 or 全部行遍历一遍)
168     # 循环迭代结束 或者 循环遍历全部alpha后,alphaPairs仍是没变化
169     while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)): 170         alphaPairsChanged = 0 171         # 当entireSet=true or 非边界alpha对没有了;就开始寻找 alpha对,而后决定是否要进行else。
172         if entireSet: 173             # 在数据集上遍历全部可能的alpha
174             for i in range(oS.m): 175                 # 是否存在alpha对,存在就+1
176                 alphaPairsChanged += innerL(i, oS) 177                 print("fullSet, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 178             iter += 1
179         # 对已存在 alpha对,选出非边界的alpha值,进行优化。
180         else: 181             # 遍历全部的非边界alpha值,也就是不在边界0或C上的值。
182             nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] 183             for i in nonBoundIs: 184                 alphaPairsChanged += innerL(i, oS) 185                 print("non-bound, iter: %d i:%d, pairs changed %d" % (iter, i, alphaPairsChanged)) 186             iter += 1
187         # 若是找到alpha对,就优化非边界alpha值,不然,就从新进行寻找,若是寻找一遍 遍历全部的行仍是没找到,就退出循环。
188         if entireSet: 189             entireSet = False  # toggle entire set loop
190         elif (alphaPairsChanged == 0): 191             entireSet = True 192         print("iteration number: %d" % iter) 193     return oS.b, oS.alphas 194 
195 # 7. 经过已知数据点和拉格朗日乘子得到分割超平面参数w
196 def calcWs(alphas, dataArr, classLabels): 197     X = mat(dataArr) 198     labelMat = mat(classLabels).transpose() 199     m, n = shape(X) 200     w = zeros((n, 1)) 201     for i in range(m): 202         w += multiply(alphas[i] * labelMat[i], X[i, :].T) 203     return w 204 
205 # 8. 绘制数据点,分割线,支持向量
206 def drawSvm(alphas,dataArr,labelArr): 207     # 对数据点进行分类
208     classifiedPts = {'+1': [], '-1': []} 209     # zip 进行封装,返回的是一个对象,可用list将其转换为列表查看
210     for point, label in zip(dataArr, labelArr): 211         if label == 1.0: 212             classifiedPts['+1'].append(point) 213         else: 214             classifiedPts['-1'].append(point) 215     ax = plt.figure().add_subplot(111) 216     # 绘制数据点
217     import numpy as ny 218     for label, pts in classifiedPts.items(): 219         pts = ny.array(pts) 220         ax.scatter(pts[:, 0], pts[:, 1], label=label) 221     # 绘制分割线
222     w = calcWs(alphas, dataArr, labelArr) 223     x1, _ = max(dataArr, key=lambda x: x[0]) 224     x2, _ = min(dataArr, key=lambda x: x[0]) 225     a1, a2 = w 226     y1, y2 = (-float(b) - a1[0]*x1)/a2[0], (-float(b) -a1[0]*x2)/a2[0] # b是矩阵类型,a1是列表类型,x1是float,注意不一样形式之间的转换
227  ax.plot([x1, x2], [y1, y2]) 228     # 绘制支持向量
229     for i, alpha in enumerate(alphas): 230         if abs(alpha) > 1e-3: 231             x, y = dataArr[i] 232             ax.scatter([x], [y], s=150, c='none', alpha=0.7, 233                        linewidth=1.5, edgecolor='#AB3319') 234  plt.show() 235 
236 # 9.测试函数
237 def testRbf(k1=1.5): 238     dataArr, labelArr = loadImages('trainingDigits') 239     b, alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000,('rbf',k1)) 240     datMat=mat(dataArr) 241     labelMat=mat(labelArr).transpose() 242     svInd=nonzero(alphas.A>0)[0] 243     sVs=datMat[svInd] 244     labelSV=labelMat[svInd] 245     print("there are %d Support Vectors" %shape(sVs)[0]) 246     # 如何利用核函数进行分类
247     m,n=shape(datMat) 248     errorCount=0 249     for i in range(m): 250         kenelEval=kernelTrans(sVs,datMat[i,:],('rbf',k1)) 251         predict=kenelEval.T*multiply(labelSV,alphas[svInd])+b 252         if sign(predict)!=sign(labelArr[i]): 253             errorCount+=1
254     print("the training error rate is %f"%(float(errorCount)/m)) 255     dataArr,labelArr=loadImages('testDigits') 256     errorCount=0 257     datMat=mat(dataArr);labelMat=mat(labelArr).transpose() 258     m,n=shape(datMat) 259     for i in range(m): 260         kenelEval = kernelTrans(sVs, datMat[i, :], ('rbf', k1)) 261     predict = kenelEval.T * multiply(labelSV, alphas[svInd] )+ b 262     if sign(predict) != sign(labelArr[i]): 263         errorCount += 1
264     print("the training error rate is %f" % (float(errorCount) / m)) 265 
266 testRbf()
View Code

能够尝试改变k1的值来观察速度和准确率的变化,下图是k1=1.3结果。

 本文相关数据集:https://pan.baidu.com/s/1EMODwE2qTtxKEVa6jcrVSw    y3xs

相关文章
相关标签/搜索