1、生成数据集
class moon_data_class(object):
def __init__(self,N,d,r,w):
self.N=N
self.w=w
self.d=d
self.r=r
def 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.N
N = self.N
r = self.r
w2 = self.w/2
d = self.d
done = True
data = np.empty(0)
while done:
#generate Rectangular data
tmp_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 ---upper
idx = 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 ----down
data_t = np.empty([N, 2])
data_t[:, 0] = data[0:N, 0] + r
data_t[:, 1] = -data[0:N, 1] - d
db_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 = 0
self.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.0
for i in range(m):
diff = train_x[i, :] - train_x_i
kernel_value[i] = np.exp(diff * diff.T / (sigma))
else: # 不使用核函数
kernel_value = train_x * train_x_i.T
return kernel_value
def 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_matrix
def 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_k
def 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_i
output: 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 = 0
alpha_j = 0
error_j = 0
if len(candidateAlphaList) > 1:
for alpha_k in candidateAlphaList:
if alpha_k == alpha_i:
continue
error_k = cal_error(svm, alpha_k)
if abs(error_k - error_i) > maxStep:
maxStep = abs(error_k - error_i)
alpha_j = alpha_k
error_j = error_k
else: # 随机选择
alpha_j = alpha_i
while 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_j
def 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、计算eta
eta = 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_j
svm.alphas[alpha_j] -= svm.train_y[alpha_j] * (error_i - error_j) / eta
# 5、确定最终的alpha_j
if svm.alphas[alpha_j] > H:
svm.alphas[alpha_j] = H
if 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_i
svm.alphas[alpha_i] += svm.train_y[alpha_i] * svm.train_y[alpha_j] \
* (alpha_j_old - svm.alphas[alpha_j])
# 8、更新b
b1 = 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 = b1
elif (0 < svm.alphas[alpha_j]) and (svm.alphas[alpha_j] < svm.C):
svm.b = b2
else:
svm.b = (b1 + b2) / 2.0
# 9、更新error
update_error_tmp(svm, alpha_j)
update_error_tmp(svm, alpha_i)
return 1
else:
return 0
def 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 = True
alpha_pairs_changed = 0
iteration = 0
while (iteration < max_iter) and ((alpha_pairs_changed > 0) or entireSet):
print("\t iterration: ", iteration)
alpha_pairs_changed = 0
if entireSet:
# 对所有的样本
for x in range(svm.n_samples):
alpha_pairs_changed += choose_and_update(svm, x)
iteration += 1
else:
# 非边界样本
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 = False
elif alpha_pairs_changed == 0:
entireSet = True
return svm
def 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.b
return predict
def 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.0
for i in range(n_samples):
# 对每一个样本得到预测值
predict=svm_predict(svm, test_x[i, :])
# 判断每一个样本的预测值与真实值是否一致
if np.sign(predict) == np.sign(test_y[i]):
correct += 1
accuracy = correct / n_samples
return accuracy
3、运行测试
if __name__ == "__main__":
# 1、导入训练数据
#dataSet, labels = load_data_libsvm("heart_scale")
step=0
color=['.r','.g','.b','.y']#颜色种类
dcolor=['*r','*g','*b','*y']#颜色种类
frames = []
N = 400
d = -5
r = 10
width = 6
data_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_pos
dataSet = np.mat(input_cells)
labels = np.mat(label)
# 2、训练SVM模型
C = 0.001
toler = 0.1
maxIter = 1000
kernel_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 = 0
y_p_old = 0
for 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 2018
@author: ASUS
"""
import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
from sklearn import datasets
class moon_data_class(object):
def __init__(self,N,d,r,w):
self.N=N
self.w=w
self.d=d
self.r=r
def 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.N
N = self.N
r = self.r
w2 = self.w/2
d = self.d
done = True
data = np.empty(0)
while done:
#generate Rectangular data
tmp_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 ---upper
idx = 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 ----down
data_t = np.empty([N, 2])
data_t[:, 0] = data[0:N, 0] + r
data_t[:, 1] = -data[0:N, 1] - d
db_moon = np.concatenate((db_moon, data_t), axis=0)
return db_moon
sess = 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 = 200
d = -4
r = 10
width = 6
data_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 = 250
x_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() +1
y_min, y_max = x_vals[:,1].min() - 1, x_vals[:,1].max() +1
xx, 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、运行结果