1、生成数据集
class moon_data_class(object):def __init__(self,N,d,r,w):self.N=Nself.w=wself.d=dself.r=rdef sgn(self,x):if(x>0):return 1;else:return -1;def sig(self,x):return 1.0/(1+np.exp(x))def dbmoon(self):N1 = 10*self.NN = self.Nr = self.rw2 = self.w/2d = self.ddone = Truedata = np.empty(0)while done:#generate Rectangular datatmp_x = 2*(r+w2)*(np.random.random([N1, 1])-0.5)tmp_y = (r+w2)*np.random.random([N1, 1])tmp = np.concatenate((tmp_x, tmp_y), axis=1)tmp_ds = np.sqrt(tmp_x*tmp_x + tmp_y*tmp_y)#generate double moon data ---upperidx = np.logical_and(tmp_ds > (r-w2), tmp_ds < (r+w2))idx = (idx.nonzero())[0]if data.shape[0] == 0:data = tmp.take(idx, axis=0)else:data = np.concatenate((data, tmp.take(idx, axis=0)), axis=0)if data.shape[0] >= N:done = False#print (data)db_moon = data[0:N, :]#print (db_moon)#generate double moon data ----downdata_t = np.empty([N, 2])data_t[:, 0] = data[0:N, 0] + rdata_t[:, 1] = -data[0:N, 1] - ddb_moon = np.concatenate((db_moon, data_t), axis=0)return db_moon
2、SVM算法
class SVM:def __init__(self, dataSet, labels, C, toler, kernel_option):self.train_x = dataSet # 训练特征self.train_y = labels # 训练标签self.C = C # 惩罚参数self.toler = toler# 迭代的终止条件之一self.n_samples = np.shape(dataSet)[0] # 训练样本的个数self.alphas = np.mat(np.zeros((self.n_samples, 1))) # 拉格朗日乘子self.b = 0self.error_tmp = np.mat(np.zeros((self.n_samples, 2))) # 保存E的缓存self.kernel_opt = kernel_option # 选用的核函数及其参数self.kernel_mat = calc_kernel(self.train_x, self.kernel_opt) # 核函数的输出def cal_kernel_value(train_x, train_x_i, kernel_option):'''样本之间的核函数的值input: train_x(mat):训练样本train_x_i(mat):第i个训练样本kernel_option(tuple):核函数的类型以及参数output: kernel_value(mat):样本之间的核函数的值'''kernel_type = kernel_option[0] # 核函数的类型，分为rbf和其他m = np.shape(train_x)[0] # 样本的个数kernel_value = np.mat(np.zeros((m, 1)))if kernel_type == 'rbf': # rbf核函数sigma = kernel_option[1]if sigma == 0:sigma = 1.0for i in range(m):diff = train_x[i, :] - train_x_ikernel_value[i] = np.exp(diff * diff.T / (sigma))else: # 不使用核函数kernel_value = train_x * train_x_i.Treturn kernel_valuedef calc_kernel(train_x, kernel_option):'''计算核函数矩阵input: train_x(mat):训练样本的特征值kernel_option(tuple):核函数的类型以及参数output: kernel_matrix(mat):样本的核函数的值'''m = np.shape(train_x)[0] # 样本的个数kernel_matrix = np.mat(np.zeros((m, m))) # 初始化样本之间的核函数值for i in range(m):kernel_matrix[:, i] = cal_kernel_value(train_x, train_x[i, :], kernel_option)return kernel_matrixdef cal_error(svm, alpha_k):'''误差值的计算input: svm:SVM模型alpha_k(int):选择出的变量output: error_k(float):误差值'''output_k = float(np.multiply(svm.alphas, svm.train_y).T * svm.kernel_mat[:, alpha_k] + svm.b)error_k = output_k - float(svm.train_y[alpha_k])return error_kdef update_error_tmp(svm, alpha_k):'''重新计算误差值input: svm:SVM模型alpha_k(int):选择出的变量output: 对应误差值'''error = cal_error(svm, alpha_k)svm.error_tmp[alpha_k] = [1, error]def select_second_sample_j(svm, alpha_i, error_i):'''选择第二个样本input: svm:SVM模型alpha_i(int):选择出的第一个变量error_i(float):E_ioutput: alpha_j(int):选择出的第二个变量error_j(float):E_j'''# 标记为已被优化svm.error_tmp[alpha_i] = [1, error_i]candidateAlphaList = np.nonzero(svm.error_tmp[:, 0].A)[0]maxStep = 0alpha_j = 0error_j = 0if len(candidateAlphaList) > 1:for alpha_k in candidateAlphaList:if alpha_k == alpha_i: continueerror_k = cal_error(svm, alpha_k)if abs(error_k - error_i) > maxStep:maxStep = abs(error_k - error_i)alpha_j = alpha_kerror_j = error_kelse: # 随机选择alpha_j = alpha_iwhile alpha_j == alpha_i:alpha_j = int(np.random.uniform(0, svm.n_samples))error_j = cal_error(svm, alpha_j)return alpha_j, error_jdef choose_and_update(svm, alpha_i):'''判断和选择两个alpha进行更新input: svm:SVM模型alpha_i(int):选择出的第一个变量'''error_i = cal_error(svm, alpha_i) # 计算第一个样本的E_i# 判断选择出的第一个变量是否违反了KKT条件if (svm.train_y[alpha_i] * error_i < -svm.toler) and (svm.alphas[alpha_i] < svm.C) or\(svm.train_y[alpha_i] * error_i > svm.toler) and (svm.alphas[alpha_i] > 0):# 1、选择第二个变量alpha_j, error_j = select_second_sample_j(svm, alpha_i, error_i)alpha_i_old = svm.alphas[alpha_i].copy()alpha_j_old = svm.alphas[alpha_j].copy()# 2、计算上下界if svm.train_y[alpha_i] != svm.train_y[alpha_j]:L = max(0, svm.alphas[alpha_j] - svm.alphas[alpha_i])H = min(svm.C, svm.C + svm.alphas[alpha_j] - svm.alphas[alpha_i])else:L = max(0, svm.alphas[alpha_j] + svm.alphas[alpha_i] - svm.C)H = min(svm.C, svm.alphas[alpha_j] + svm.alphas[alpha_i])if L == H:return 0# 3、计算etaeta = 2.0 * svm.kernel_mat[alpha_i, alpha_j] - svm.kernel_mat[alpha_i, alpha_i] \- svm.kernel_mat[alpha_j, alpha_j]if eta >= 0:return 0# 4、更新alpha_jsvm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta# 5、确定最终的alpha_jif svm.alphas[alpha_j] > H:svm.alphas[alpha_j] = Hif svm.alphas[alpha_j] < L:svm.alphas[alpha_j] = L# 6、判断是否结束if abs(alpha_j_old - svm.alphas[alpha_j]) < 0.00001:update_error_tmp(svm, alpha_j)return 0# 7、更新alpha_isvm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \* (alpha_j_old - svm.alphas[alpha_j])# 8、更新bb1 = svm.b - error_i - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \* svm.kernel_mat[alpha_i, alpha_i] \- svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \* svm.kernel_mat[alpha_i, alpha_j]b2 = svm.b - error_j - svm.train_y[alpha_i] * (svm.alphas[alpha_i] - alpha_i_old) \* svm.kernel_mat[alpha_i, alpha_j] \- svm.train_y[alpha_j] * (svm.alphas[alpha_j] - alpha_j_old) \* svm.kernel_mat[alpha_j, alpha_j]if (0 < svm.alphas[alpha_i]) and (svm.alphas[alpha_i] < svm.C):svm.b = b1elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):svm.b = b2else:svm.b = (b1 + b2) / 2.0# 9、更新errorupdate_error_tmp(svm, alpha_j)update_error_tmp(svm, alpha_i)return 1else:return 0def SVM_training(train_x, train_y, C, toler, max_iter, kernel_option = ('rbf', 0.431029)):'''SVM的训练input: train_x(mat):训练数据的特征train_y(mat):训练数据的标签C(float):惩罚系数toler(float):迭代的终止条件之一max_iter(int):最大迭代次数kerner_option(tuple):核函数的类型及其参数output: svm模型'''# 1、初始化SVM分类器svm = SVM(train_x, train_y, C, toler, kernel_option)# 2、开始训练entireSet = Truealpha_pairs_changed = 0iteration = 0while (iteration < max_iter) and ((alpha_pairs_changed > 0) or entireSet):print("\t iterration: ", iteration)alpha_pairs_changed = 0if entireSet:# 对所有的样本for x in range(svm.n_samples):alpha_pairs_changed += choose_and_update(svm, x)iteration += 1else:# 非边界样本bound_samples = []for i in range(svm.n_samples):if svm.alphas[i,0] > 0 and svm.alphas[i,0] < svm.C:bound_samples.append(i)for x in bound_samples:alpha_pairs_changed += choose_and_update(svm, x)iteration += 1# 在所有样本和非边界样本之间交替if entireSet:entireSet = Falseelif alpha_pairs_changed == 0:entireSet = Truereturn svmdef svm_predict(svm, test_sample_x):'''利用SVM模型对每一个样本进行预测input: svm:SVM模型test_sample_x(mat):样本output: predict(float):对样本的预测'''# 1、计算核函数矩阵kernel_value = cal_kernel_value(svm.train_x, test_sample_x, svm.kernel_opt)# 2、计算预测值predict = kernel_value.T * np.multiply(svm.train_y, svm.alphas) #+ svm.breturn predictdef cal_accuracy(svm, test_x, test_y):'''计算预测的准确性input: svm:SVM模型test_x(mat):测试的特征test_y(mat):测试的标签output: accuracy(float):预测的准确性'''n_samples = np.shape(test_x)[0] # 样本的个数correct = 0.0for i in range(n_samples):# 对每一个样本得到预测值predict=svm_predict(svm, test_x[i, :])# 判断每一个样本的预测值与真实值是否一致if np.sign(predict) == np.sign(test_y[i]):correct += 1accuracy = correct / n_samplesreturn accuracy
3、运行测试
if __name__ == "__main__":# 1、导入训练数据#dataSet, labels = load_data_libsvm("heart_scale")step=0color=['.r','.g','.b','.y']#颜色种类dcolor=['*r','*g','*b','*y']#颜色种类frames = []N = 400d = -5r = 10width = 6data_source = moon_data_class(N, d, r, width)data = data_source.dbmoon()# x0 = [1 for x in range(1,401)]input_cells = np.array([np.reshape(data[0:2*N, 0], len(data)), np.reshape(data[0:2*N, 1], len(data))]).transpose()labels_pre = [[-1.] for y in range(1, 401)]labels_pos = [[1. ] for y in range(1, 401)]label=labels_pre+labels_posdataSet = np.mat(input_cells)labels = np.mat(label)# 2、训练SVM模型C = 0.001toler = 0.1maxIter = 1000kernel_option = ('rbf', -10)svm_model = SVM_training(dataSet, labels, C, toler, maxIter,kernel_option)# 3、计算训练的准确性accuracy = cal_accuracy(svm_model, dataSet, labels) print("The training accuracy is: %.3f%%" % (accuracy * 100))# 4、保存最终的SVM模型print("------------ 4、save model ----------------")#svm.save_svm_model(svm_model, "model_file")test_x = []test_y = []test_p = []predict = 0y_p_old = 0for x in np.arange(-15.,25.,1):for y in np.arange(-15.,25.,1):predict = svm_predict(svm_model, np.array([x, y]))#print(predict)y_p = np.sign(predict)[0, 0]#y_p =get_prediction(np.array([x, y]),svm_model)#y_p =float(y_p)if(y_p_old > 0 and y_p < 0):test_x.append(x)test_y.append(y)test_p.append([y_p_old,y_p])y_p_old = y_p#画决策边界plt.plot( test_x, test_y, 'g--') plt.plot(data[0:N, 0], data[0:N, 1], 'r*', data[N:2*N, 0], data[N:2*N, 1], 'b*')plt.show()
4、运行结果
5、利用tensorflow实现
# -*- coding: utf-8 -*-"""Created on Fri Nov 9 22:00:44 @author: ASUS"""import matplotlib.pyplot as pltimport numpy as npimport tensorflow as tffrom sklearn import datasetsclass moon_data_class(object):def __init__(self,N,d,r,w):self.N=Nself.w=wself.d=dself.r=rdef sgn(self,x):if(x>0):return 1;else:return -1;def sig(self,x):return 1.0/(1+np.exp(x))def dbmoon(self):N1 = 10*self.NN = self.Nr = self.rw2 = self.w/2d = self.ddone = Truedata = np.empty(0)while done:#generate Rectangular datatmp_x = 2*(r+w2)*(np.random.random([N1, 1])-0.5)tmp_y = (r+w2)*np.random.random([N1, 1])tmp = np.concatenate((tmp_x, tmp_y), axis=1)tmp_ds = np.sqrt(tmp_x*tmp_x + tmp_y*tmp_y)#generate double moon data ---upperidx = np.logical_and(tmp_ds > (r-w2), tmp_ds < (r+w2))idx = (idx.nonzero())[0]if data.shape[0] == 0:data = tmp.take(idx, axis=0)else:data = np.concatenate((data, tmp.take(idx, axis=0)), axis=0)if data.shape[0] >= N:done = False#print (data)db_moon = data[0:N, :]#print (db_moon)#generate double moon data ----downdata_t = np.empty([N, 2])data_t[:, 0] = data[0:N, 0] + rdata_t[:, 1] = -data[0:N, 1] - ddb_moon = np.concatenate((db_moon, data_t), axis=0)return db_moonsess = tf.Session()(x_vals, y_vals) = datasets.make_circles(n_samples=500, factor=.5,noise=.1)y_vals = np.array([1 if y==1 else -1 for y in y_vals])N = 200d = -4r = 10width = 6data_source = moon_data_class(N, d, r, width)data = data_source.dbmoon()# x0 = [1 for x in range(1,401)]x_vals = np.array([np.reshape(data[0:2*N, 0], len(data)), np.reshape(data[0:2*N, 1], len(data))]).transpose()labels_pre = [-1 for y in range(1, 201)]labels_pos = [1 for y in range(1, 201)]y_vals = np.array(labels_pre+labels_pos)class1_x = data[0:N, 0]class1_y = data[0:N, 1]class2_x = data[N:2*N, 0]class2_y = data[N:2*N, 1]batch_size = 250x_data = tf.placeholder(shape = [None,2],dtype = tf.float32)y_target = tf.placeholder(shape = [None,1],dtype = tf.float32)prediction_grid = tf.placeholder(shape = [None, 2],dtype = tf.float32)b = tf.Variable(tf.random_normal(shape = [1, batch_size]))gamma = tf.constant(-0.05)dist = tf.reduce_sum(tf.square(x_data),1)dist = tf.reshape(dist,[-1,1])sq_dists = tf.add(tf.subtract(dist, tf.multiply(2., tf.matmul(x_data, tf.transpose(x_data)))),tf.transpose(dist))my_kernel = tf.exp(tf.multiply(gamma, tf.abs(sq_dists)))model_output = tf.matmul(b,my_kernel)first_term= tf.reduce_sum(b)b_vec_cross = tf.matmul(tf.transpose(b),b)y_target_cross = tf.matmul(y_target,tf.transpose(y_target))second_term = tf.reduce_sum(tf.multiply(my_kernel, tf.multiply(b_vec_cross,y_target_cross)))loss = tf.negative(tf.subtract(first_term, second_term))rA = tf.reshape(tf.reduce_sum(tf.square(x_data),1),[-1,1])rB = tf.reshape(tf.reduce_sum(tf.square(prediction_grid),1),[-1,1])pred_sq_dist = tf.add(tf.subtract(rA, tf.multiply(2., tf.matmul(x_data, tf.transpose(prediction_grid)))),tf.transpose(rB))pred_kernel = tf.exp(tf.multiply(gamma, tf.abs(pred_sq_dist)))prediction_output = tf.matmul(tf.multiply(tf.transpose(y_target),b), pred_kernel)prediction = tf.sign(prediction_output - tf.reduce_mean(prediction_output))accuracy = tf.reduce_mean(tf.cast(tf.equal(tf.squeeze(prediction), tf.squeeze(y_target)), tf.float32))my_opt = tf.train.GradientDescentOptimizer(0.01)train_step = my_opt.minimize(loss)init = tf.global_variables_initializer()sess.run(init)loss_vec = []batch_accuracy = []for i in range(5000):rand_index = np.random.choice(len(x_vals),size=batch_size)rand_x = x_vals[rand_index]rand_y = np.transpose([y_vals[rand_index]])sess.run(train_step, feed_dict={x_data:rand_x, y_target:rand_y})temp_loss = sess.run(loss, feed_dict={x_data:rand_x, y_target:rand_y})loss_vec.append(temp_loss)acc_temp = sess.run(accuracy,feed_dict ={x_data:rand_x, y_target:rand_y,prediction_grid:rand_x})batch_accuracy.append(acc_temp)if (i+1)%100==0:print('Step # ' + str(i+1))print('Loss = ' + str(temp_loss))x_min, x_max = x_vals[:,0].min() - 1, x_vals[:,0].max() +1y_min, y_max = x_vals[:,1].min() - 1, x_vals[:,1].max() +1xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02), np.arange(y_min, y_max, 0.02))grid_points = np.c_[xx.ravel(), yy.ravel()][grid_predictions] = sess.run(prediction,feed_dict ={x_data:rand_x, y_target:rand_y,prediction_grid:grid_points})grid_predictions = grid_predictions.reshape(xx.shape)plt.contourf(xx,yy,grid_predictions, cmap=plt.cm.Paired,alpha=0.8)plt.plot(class1_x,class1_y, 'ro',label='I. setosa')plt.plot(class2_x,class2_y, 'rx',label='Non setosa')plt.legend(loc='lower right')plt.ylim([-15,15])plt.xlim([-15,25])plt.show()
6、运行结果
