數(shù)學(xué)建模——支持向量機(jī)模型詳解Python代碼

from numpy import * import random import matplotlib.pyplot as plt import numpydef kernelTrans(X,A,kTup): # 核函數(shù)(此例未使用)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] # 樣本點個數(shù)self.alphas = mat(zeros((self.m,1))) # 產(chǎn)生m個拉格郎日乘子，組成一個m×1的矩陣self.b =0 # 決策面的截距self.eCache = mat(zeros((self.m,2))) # 產(chǎn)生m個誤差 E=f(x)-y ,設(shè)置成m×2的矩陣，矩陣第一列是標(biāo)志位，標(biāo)志為1就是E計算好了，第二列是誤差E# self.K = mat(zeros((self.m,self.m)))# for i in range(self.m): # K[,]保存的是任意樣本之間的相似度(用高斯核函數(shù)表示的相似度)# self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup)def loadDataSet(filename): # 加載數(shù)據(jù)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): # 隨機(jī)選擇一個不等于i的下標(biāo)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個樣本點，它是下標(biāo)，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 # 最大步長的因子的下標(biāo)maxDeltaE = 0 # 最大步長Ej = 0 # 最大步長的因子的誤差oS.eCache[i] = [1,Ei]valiEcacheList = nonzero(oS.eCache[:,0].A)[0] # nonzero結(jié)果是兩個array數(shù)組，第一個數(shù)組是不為0的元素的x坐標(biāo)，第二個數(shù)組是該位置的y坐標(biāo)# 此處尋找誤差矩陣第一列不為0的數(shù)的下標(biāo)print("valiEcacheList is {}".format(valiEcacheList))if (len(valiEcacheList))>1:for k in valiEcacheList: # 遍歷所有計算好的Ei的下標(biāo)，valiEcacheLIst保存了所有樣本點的E，計算好的有效位置是1,沒計算好的是0if k == i:continueEk = calcEk(oS,k)deltaE = abs(Ei-Ek) # 距離第一個拉格朗日乘子a1絕對值最遠(yuǎn)的作為第二個朗格朗日乘子a2if deltaE>maxDeltaE:maxK = k # 記錄選中的這個乘子a2的下標(biāo)maxDeltaE = deltaE # 記錄他倆的絕對值Ej = Ek # 記錄a2此時的誤差return maxK, Ejelse: # 如果是第一次循環(huán)，隨機(jī)選擇一個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個樣本點的誤差計算出來，并存入誤差矩陣，有效位置設(shè)為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("第二個因子的坐標(biāo){}".format(j))alphaIold = oS.alphas[i].copy() # 用了淺拷貝， alphaIold 就是old a1,對應(yīng)公式alphaJold = oS.alphas[j].copy()if oS.labelMat[i] != oS.labelMat[j]: # 也是根據(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# 當(dāng)同類別時 ai+aj = 常數(shù) ai是不滿足KKT的，假設(shè)ai=0,需增大它，那么就得減少aj，aj已經(jīng)是0了，不能最小了，所以此情況不允許發(fā)生# 當(dāng)不同類別時 ai-aj=常數(shù)，ai是不滿足KKT的，ai=0,aj=C,ai需增大，它則aj也會變大，但是aj已經(jīng)是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: # 這里跟公式正好差了一個負(fù)號，所以對應(yīng)公式里的 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) # 根據(jù)公式，看看得到的a2^new是否在上下限之內(nèi)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]) # 根據(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) # 個人認(rèn)為更新誤差應(yīng)在更新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] # 這時數(shù)組相乘，里面其實是True 和False的數(shù)組，得出來的是# 大于0并且小于C的alpha的下標(biāo)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 # 當(dāng)?shù)诙N遍歷方式alpha不再變化，那么繼續(xù)第一種方式掃描，第一種方式不再變化，此時alphachanged為0且entireSet為false，退出循環(huán)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 # 防止坐標(biāo)軸的‘-’變?yōu)榉綁Km = 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() # 創(chuàng)建畫布ax = fig.add_subplot(111)ax.set_title("SVM-Classify") # 設(shè)置圖片標(biāo)題ax.set_xlabel("x") # 設(shè)置坐標(biāo)名稱ax.set_ylabel("y")ax1=ax.scatter(red_points_x, red_points_y, s=30,c='red', marker='s') #s是shape大小，c是顏色，marker是形狀，'s'代表是正方形,默認(rèn)'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') # 標(biāo)注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)

可參考數(shù)據(jù)集

-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 -1

總結(jié)

以上是生活随笔為你收集整理的数学建模——支持向量机模型详解Python代码的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。