支持向量机 

  1 import numpy as np
  2 
  3 def clipalpha(aj,H,L):
  4     """
  5     剪枝
  6     """
  7     if aj>H:
  8         aj=H
  9     if aj<L:
 10         aj=L
 11     return aj
 12     
 13 def selectJrand(i,m):
 14     """
 15     取得一个不等于i的0-m之间的随机值
 16     """
 17     j=i
 18     while (j==i):
 19         j=np.random.randint(0,m)
 20     return j
 21 
 22 class optstruct:
 23     def __init__(self,datamat,classlabels,C,toler):
 24         self.x=datamat
 25         self.labelmat=classlabels
 26         self.C=C
 27         self.tol=toler
 28         self.m=np.shape(datamat)[0]
 29         self.alphas=np.mat(np.zeros((self.m,1)))
 30         self.b=0
 31         self.ecache=np.mat(np.zeros((self.m,2)))
 32         
 33 def calcek(obj,k):
 34     """
 35     计算Ek，Ek为预测值与真实输入之差
 36     """
 37     fxk=float(np.multiply(obj.alphas,obj.labelmat).T*(obj.x*obj.x[k,:].T))\
 38     +obj.b
 39     #g(xi)-yi=∑(ajyjxixj+b)-yi 数乘后矩阵乘
 40     Ek=fxk-float(obj.labelmat[k])
 41     return Ek
 42     
 43 def select(i,obj,Ei):
 44     """
 45     找出|Ej-Ei|最大的值，并返回Ej
 46     """
 47     maxK=-1;maxdeltae=0;Ej=0
 48     obj.ecache[i]=[1,Ei]
 49     valecache=np.nonzero(obj.ecache[:,0].A)[0]
 50     #将ecache内的alpha值，即第一列取出，并取出其中的非零值
 51     if (len(valecache))>1:
 52         for k in valecache:
 53             if k==i:
 54                 continue
 55             Ek=calcek(obj,k)
 56             deltaE=np.abs(Ei-Ek)
 57             if (deltaE>maxdeltae):
 58                 maxK=k
 59                 maxdeltae=deltaE
 60                 Ej=Ek
 61         return maxK,Ej
 62     else:
 63         j=np.selectJrand(i,obj.m)
 64         Ej=calcek(obj,j)
 65     return j,Ej
 66 
 67 def update(obj,k):
 68     """
 69     更新完alpha后，需要更新Ek
 70     """
 71     Ek=calcek(obj,k)
 72     obj.ecache=[1,Ek]
 73         
 74 def innerl(i,obj):
 75     """
 76     优化
 77     """
 78     Ei=calcek(obj,i)
 79     if ((obj.labelmat[i]*Ei<-obj.tol) and (obj.alphas[i]<obj.C) or ((\
 80         obj.label.mat[i]*Ei>obj.tol and obj.alphas[i]>0))):
 81         #由于KKT条件，alpha<C时，松弛变量必为0，alpha>0时，对应的向量必为
 82         #支持向量，即yi(w*xi+b)<=1。因此找到满足上述条件的alpha值就是违反
 83         #KKT条件的alpha值
 84         j,Ej=select(i,obj,Ei)
 85         alphaIold=obj.alphas[i]
 86         alphaJold=obj.alphas[j]
 87         if (obj.labelmat[i] != obj.labelmat[j]):
 88             L=max(0,obj.alphas[j]-obj.alphas[i])
 89             H=min(obj.C,obj.C+obj.alphas[j]-obj.alphas[i])
 90         else:
 91             L=max(0,obj.alphas[j]+obj.alphas[i]-obj.C)
 92             H=min(obj.C,obj.alphas[j]+obj.alphas[i])
 93             #线性规划的另一种表达方式，需要考虑alphas[j]在不同情况下的取值
 94         if L==H:
 95             return 0
 96         #此时alphas[i]-alphas[j]=C,没有优化空间
 97         eta=2*obj.x[i,:]*obj.x[j,:].T-obj.x[i,:]*obj.x[i,:].T-obj.x[j,:]\
 98         *obj.x[j,:].T
 99         if eta>=0:
100             return 0
101         #-eta＜0，即K11+K22-2K12>0:说明核函数非正定，一般不可能出现。
102         obj.alphas[j]-=obj.labelmat[j]*(Ei-Ej)/eta
103         obj.alphas[j]=clipalpha(obj.alphas[j],H,L)
104         update(obj,j)
105         if (abs(obj.alphas[j]-alphaJold)<0.00001):
106             return 0
107         obj.alphas[i]+=obj.labelmat[j]*obj.labelmat[i]*(alphaJold-\
108                   obj.alphas[j])
109         #保持总量不变，alphas[j]增加多少，alphas[i]就减少多少
110         update(obj,i)
111         b1=obj.b-Ei-obj.labelmat[i]*(obj.alphas[i]-alphaIold)*obj.x[i,:]*\
112         obj.x[i,:].T-obj.labelmat[j]*(obj.alphas[j]-alphaJold)*obj.x[i,:]\
113         *obj.x[j,:].T
114         b2=obj.b-Ej-obj.labelmat[i]*(obj.alphas[i]-alphaIold)*obj.x[i,:]*\
115         obj.x[j,:].T-obj.labelmat[j]*(obj.alphas[j]-alphaJold)*obj.x[j,:]\
116         *obj.x[j,:].T
117         if (0<obj.alphas[i]) and (obj.C>obj.alphas[i]):
118             obj.b=b1
119         elif (0<obj.alphas[j]) and (obj.C>obj.alphas[j]):
120             obj.b=b2
121         else:
122             obj.b=(b1+b2)/2
123         return 1
124     else:
125         return 0
126     
127 def smosvm(datamatin,classlabels,c,toler,maxiter,ktup=('lin',0)):
128     obj=optstruct(np.mat(datamatin),np.mat(classlabels).transpose(),c,toler)
129     svmiter=0
130     entireset=True
131     alphapaieschanged=0
132     while(svmiter<maxiter) and ((alphapaieschanged>0) or (entireset)):
133         alphapaieschanged=0
134         if entireset:
135             for i in range(obj.m):
136                 alphapaieschanged+=innerl(i,obj)
137             svmiter+=1
138         else:
139             nonboundis=np.nonzero((obj.alphas.A>0)*(obj.alphas.A<c))[0]
140             for i in nonboundis:
141                alphapaieschanged+=innerl(i,obj)
142             svmiter+=1
143         if entireset:
144             entireset=False
145         elif (alphapaieschanged==0):
146             entireset=True
147     return obj.b,obj.alphas

 

posted @ 2017-10-21 10:40  尤格索托斯  阅读(166)  评论(0)    收藏  举报