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):
8 def __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
24 def kernel_function(self,x1, x2): # 核函数
25 return np.dot(x1, x2)
26
27 def 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
32 def 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
39 def 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
53 def gx(self, x): # 判别函数g(x)
54 return np.dot(self.w, x) + self.b
55
56 def 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
66 def 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
79 def 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
93 def 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 return
101 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=0
105
106 def train(self, x, y, max_iternum=100): # SMO算法
107 x_len = len(x)
108 self.X = x
109 self.N = x_len
110 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 = 0
115 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 = 0
119 for i in range(0, max_iternum):
120 self.stop=1
121 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:
130 self.update(data[1])
131 if self.stop:
132 break
133 return self.w, self.b
134
135 def 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
143 def predict(self,x_new): # 预测
144 p=np.zeros(self.k)
145 for 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]+=1
150 else:
151 p[i]+=1
152 return self.cls[np.argmax(p)]
153
154 def 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,suby
165
166 def classfy(self,x_new): # 预测
167 y_new=self.gx(x_new)
168 cl = int(np.sign(y_new))
169 if cl == 0:
170 cl = 1
171 return cl
172
173
174 def load_data():
175 f = gzip.open('../data/mnist.pkl.gz', 'rb')
176 training_data, validation_data, test_data = cPickle.load(f)
177 f.close()
178 return (training_data, validation_data, test_data)
179
180 if __name__ == "__main__":
181 svc = SVC()
182 np.random.seed(0)
183 l=1000
184 training_data, validation_data, test_data = load_data()
185 svc.fit(training_data[0][:l],training_data[1][:l])
186 predictions = [svc.predict(a) for a in test_data[0][:l]]
187 num_correct = sum(int(a == y) for a, y in zip(predictions, test_data[1][:l]))
188 print "%s of %s values correct." % (num_correct, len(test_data[1][:l])) #72/100 #808/1000 #8194/10000(较慢）
浙公网安备 33010602011771号