1 # encoding: utf-8 2 import numpy as np 3 import matplotlib.pyplot as plt 4 import cPickle 5 import gzip 6 7 class SVC(object): 8def __init__(self, c=1.0, delta=0.001): # 初始化 9 self.N = 0 10 self.delta = delta 11 self.X = None 12 self.y = None 13 self.w = None 14 self.wn = 0 15 self.K = np.zeros((self.N, self.N)) 16 self.a = np.zeros((self.N, 1)) 17 self.b = 0 18 self.C = c 19 self.stop=1 20 self.k=0 21 self.cls=0 22 self.train_result=[] 23 24def kernel_function(self,x1, x2): # 核函数 25 return np.dot(x1, x2) 26 27def kernel_matrix(self, x): # 核矩阵 28 for i in range(0, len(x)): 29 for j in range(i, len(x)): 30 self.K[j][i] = self.K[i][j] = self.kernel_function(self.X[i], self.X[j]) 31 32def get_w(self): # 计算更新w 33 ay = self.a * self.y 34 w = np.zeros((1, self.wn)) 35 for i in range(0, self.N): 36 w += self.X[i] * ay[i] 37 return w 38 39def get_b(self, a1, a2, a1_old, a2_old): # 计算更新B 40 y1 = self.y[a1] 41 y2 = self.y[a2] 42 a1_new = self.a[a1] 43 a2_new = self.a[a2] 44 b1_new = -self.E[a1] - y1 * self.K[a1][a1] * (a1_new - a1_old) - y2 * self.K[a2][a1] * ( 45 a2_new - a2_old) + self.b 46 b2_new = -self.E[a2] - y1 * self.K[a1][a2] * (a1_new - a1_old) - y2 * self.K[a2][a2] * ( 47 a2_new - a2_old) + self.b 48 if (0 < a1_new) and (a1_new < self.C) and (0 < a2_new) and (a2_new < self.C): 49 return b1_new[0] 50 else: 51 return (b1_new[0] + b2_new[0]) / 2.0 52 53def gx(self, x): # 判别函数g(x) 54 return np.dot(self.w, x) + self.b 55 56def satisfy_kkt(self, a): # 判断样本点是否满足kkt条件 57 index = a[1] 58 if a[0] == 0 and self.y[index] * self.gx(self.X[index]) > 1: 59 return 1 60 elif a[0] < self.C and self.y[index] * self.gx(self.X[index]) == 1: 61 return 1 62 elif a[0] == self.C and self.y[index] * self.gx(self.X[index]) < 1: 63 return 1 64 return 0 65 66def clip_func(self, a_new, a1_old, a2_old, y1, y2): # 拉格朗日乘子的裁剪函数 67 if (y1 == y2): 68 L = max(0, a1_old + a2_old - self.C) 69 H = min(self.C, a1_old + a2_old) 70 else: 71 L = max(0, a2_old - a1_old) 72 H = min(self.C, self.C + a2_old - a1_old) 73 if a_new < L: 74 a_new = L 75 if a_new > H: 76 a_new = H 77 return a_new 78 79def update_a(self, a1, a2): # 更新a1,a2 80 partial_a2 = self.K[a1][a1] + self.K[a2][a2] - 2 * self.K[a1][a2] 81 if partial_a2 <= 1e-9: 82 print "error:", partial_a2 83 a2_new_unc = self.a[a2] + (self.y[a2] * ((self.E[a1] - self.E[a2]) / partial_a2)) 84 a2_new = self.clip_func(a2_new_unc, self.a[a1], self.a[a2], self.y[a1], self.y[a2]) 85 a1_new = self.a[a1] + self.y[a1] * self.y[a2] * (self.a[a2] - a2_new) 86 if abs(a1_new - self.a[a1]) < self.delta: 87 return 0 88 self.a[a1] = a1_new 89 self.a[a2] = a2_new 90 self.is_update = 1 91 return 1 92 93def update(self, first_a): # 更新拉格朗日乘子 94 for second_a in range(0, self.N): 95 if second_a == first_a: 96 continue 97 a1_old = self.a[first_a] 98 a2_old = self.a[second_a] 99 if self.update_a(first_a, second_a) == 0:100 return101 self.b= self.get_b(first_a, second_a, a1_old, a2_old)102 self.w = self.get_w()103 self.E = [self.gx(self.X[i]) - self.y[i] for i in range(0, self.N)]104 self.stop=0105 106def train(self, x, y, max_iternum=100): # SMO算法107 x_len = len(x)108 self.X = x109 self.N = x_len110 self.wn = len(x[0])111 self.y = np.array(y).reshape((self.N, 1))112 self.K = np.zeros((self.N, self.N))113 self.kernel_matrix(self.X)114 self.b = 0115 self.a = np.zeros((self.N, 1))116 self.w = self.get_w()117 self.E = [self.gx(self.X[i]) - self.y[i] for i in range(0, self.N)]118 self.is_update = 0119 for i in range(0, max_iternum):120 self.stop=1121 data_on_bound = [[x,y] for x,y in zip(self.a, range(0, len(self.a))) if x > 0 and x< self.C]122 if len(data_on_bound) == 0:123 data_on_bound = [[x,y] for x,y in zip(self.a, range(0, len(self.a)))]124 for data in data_on_bound:125 if self.satisfy_kkt(data) != 1:126 self.update(data[1])127 if self.is_update == 0:128 for data in [[x,y] for x,y in zip(self.a, range(0, len(self.a)))]:129 if self.satisfy_kkt(data) != 1:130self.update(data[1])131 if self.stop:132 break133 return self.w, self.b134 135def fit(self,x, y): # 训练模型, 一对一法k(k-1)/2个SVM进行多类分类136 self.cls, y = np.unique(y, return_inverse=True)137 self.k=len(self.cls)138 for i in range(self.k):139 for j in range(i):140 a,b=self.sub_data(x,y,i,j)141 self.train_result.append([i,j,self.train(a,b)])142 143def predict(self,x_new): # 预测144p=np.zeros(self.k)145for i,j,w in self.train_result:146 self.w=w[0]147 self.b=w[1]148 if self.classfy(x_new)==1:149 p[j]+=1150 else:151 p[i]+=1152return self.cls[np.argmax(p)]153 154def sub_data(self,x,y,i,j): # 数据分类155 subx=[]156 suby=[]157 for a,b in zip(x,y):158 if b==i:159 subx.append(a)160 suby.append(-1)161 elif b==j:162 subx.append(a)163 suby.append(1)164 return subx,suby165 166def classfy(self,x_new): # 预测167 y_new=self.gx(x_new)168 cl = int(np.sign(y_new))169 if cl == 0:170 cl = 1171 return cl172 173 174 def load_data():175f = gzip.open('../data/mnist.pkl.gz', 'rb')176training_data, validation_data, test_data = cPickle.load(f)177f.close()178return (training_data, validation_data, test_data)179 180 if __name__ == "__main__":181svc = SVC()182np.random.seed(0)183l=1000184training_data, validation_data, test_data = load_data()185svc.fit(training_data[0][:l],training_data[1][:l])186predictions = [svc.predict(a) for a in test_data[0][:l]]187num_correct = sum(int(a == y) for a, y in zip(predictions, test_data[1][:l]))188print "%s of %s values correct." % (num_correct, len(test_data[1][:l])) #72/100 #808/1000 #8194/10000(较慢）
