AnFany_SVM_Regression.py

# -*- coding：utf-8 -*-
# &Author AnFany

# SMO算法实现支持向量机回归

"""
第一部分：引入库
"""
importnumpyasnp
importmatplotlib.pyplotasplt
frompylabimportmpl

mpl.rcParams['font.sans-serif'] = ['FangSong'] # 中文字体名称
mpl.rcParams['axes.unicode_minus'] =False# 显示负号

# 引入数据
importSVM_Regression_Dataasrdata

"""
第二部分：构建核函数以及SVM的结构
"""

# 构建核函数
classKERNEL:
"""
 linear：线性 rbf：高斯 sigmoid：Sigmoid型 poly：多项式
 核函数：注意输入数据的shape以及输出数据的shape。
 xVSy包括3种情况：单样本VS单样本 单样本VS多样本 多样本VS多样本
 """

def__init__(self, polyd=3, rbfsigma=0.2, tanhbeta=0.6, tanhtheta=-0.6):
self.polyd=polyd
self.rbfsigma=rbfsigma
self.tanhbeta=tanhbeta
self.tanhtheta=tanhtheta

deftrans(self, x):
x=np.array(x)
ifx.ndim==1:
x=np.array([x])
returnx

# 线性核函数
deflinear(self, x, y): # 输出的结果shape=(len(y), len(x))
x, y=self.trans(x), self.trans(y)
iflen(x) ==1:
return (x*y).sum(axis=1, keepdims=True)
else:
sx=x.reshape(x.shape[0], -1, x.shape[1])
return (sx*y).sum(axis=2).T

# Singmoid型核函数
defsigmoid(self, x, y): # 输出的结果shape=(len(y), len(x))
x, y=self.trans(x), self.trans(y)
iflen(x) ==1:
returnnp.tanh(self.tanhbeta* ((x*y).sum(axis=1, keepdims=True)) +self.tanhtheta)
else:
sx=x.reshape(x.shape[0], -1, x.shape[1])
returnnp.tanh(self.tanhbeta* (sx*y).sum(axis=2).T+self.tanhtheta)

# 多项式核函数
defpoly(self, x, y): # 输出的结果shape=(len(y), len(x))
x, y=self.trans(x), self.trans(y)
iflen(x) ==1:
return (x*y).sum(axis=1, keepdims=True) **self.polyd
else:
sx=x.reshape(x.shape[0], -1, x.shape[1])
return (sx*y).sum(axis=2).T**self.polyd

# 高斯核函数
defrbf(self, x, y): # 输出的结果shape=(len(y), len(x))
x, y=self.trans(x), self.trans(y)
iflen(x) ==1andlen(y) ==1:
returnnp.exp(self.linear((x-y), (x-y)) / (-2*self.rbfsigma**2))
eliflen(x) ==1andlen(y) !=1:
returnnp.exp((np.power(x-y, 2)).sum(axis=1, keepdims=True) / (-2*self.rbfsigma**2))
else:
sx=x.reshape(x.shape[0], -1, x.shape[1])
returnnp.exp((np.power(sx-y, 2)).sum(axis=2).T/ (-2*self.rbfsigma**2))


# 构建SVM的结构
classSVR:
def__init__(self, feature, labels, kernel='rbf', C=0.8, toler=0.001, epsilon=0.001, times=100, eps=0.0001):
# 训练样本的属性数据、标签数据
self.feature=feature
self.labels=labels

# SMO算法变量
self.C=C
self.toler=toler

self.alphas=np.zeros(len(self.feature))
self.alphas_star=np.zeros(len(self.feature))

self.b=0
self.eps=eps# 选择拉格朗日因子
self.epsilon=epsilon

# 核函数
self.kernel=eval('KERNEL().'+kernel)

# 拉格朗日误差序列
self.errors= [self.get_error(i) foriinrange(len(self.feature))]

# 循环的最大次数
self.times=times

# 计算分割线的值,x为单样本
defline_num(self, x):
ks=self.kernel(x, self.feature)
wx=np.matrix(self.alphas-self.alphas_star) *ks
num=np.array(wx+self.b)
returnnum[0][0]

# 获得编号为i的样本对应的误差
defget_error(self, i):
x, y=self.feature[i], self.labels[i]
error=self.line_num(x) -y
returnerror

# 更改拉格朗日因子后，更新所有样本对应的误差
defupdate_errors(self):
self.errors= [self.get_error(i) foriinrange(len(self.feature))]


# 预测函数
defpredictall(self, prefeat):
ks=self.kernel(prefeat, self.feature)
wx=np.matrix(self.alphas-self.alphas_star) *ks
num=np.array(wx+self.b)
returnnum[0]


"""
第三部分：构建SMO算法需要的函数
"""


deftakestep(svr, i1, i2):
ifi1==i2:
return0
alpha1, alphas1=svr.alphas[i1], svr.alphas_star[i1]

alpha2, alphas2=svr.alphas[i2], svr.alphas_star[i2]

phi1, phi2=svr.errors[i1], svr.errors[i2]

x1, x2=svr.feature[i1], svr.feature[i2]

y1, y2=svr.labels[i1], svr.labels[i2]

k11, k12, k22=svr.kernel(x1, x1)[0][0], svr.kernel(x1, x2)[0][0], svr.kernel(x2, x2)[0][0]

eta=2*k12-k11-k22

gamma=alpha1-alphas1+alpha2-alphas2

case1=case2=case3=case4=finished=0

a1old, a1olds=alpha1, alphas1

delta_phi=phi1-phi2

whilenotfinished:
if (case1==0) and (alpha1>0or (alphas1==0anddelta_phi>0)) \
and (alpha2>0or (alphas2==0anddelta_phi<0)):
# 计算L和H (a1, a2)
L=max(0, gamma-svr.C)
H=min(svr.C, gamma)
ifL<H:
ifeta>0:
a2=alpha2-delta_phi/eta
a2=min(a2, H)
a2=max(L, a2)
a1=alpha1-a2+alpha2
else:
a2=L
a1=alpha1- (a2-alpha2)
object1=-0.5*a1*a1*eta+a1* (delta_phi+ (a1old-a1olds) *eta)
a2=H
a1=alpha1- (a2-alpha2)
object2=-0.5*a1*a1*eta+a1* (delta_phi+ (a1old-a1olds) *eta)
ifobject1>object2:
a2=L
else:
a2=H
a1=alpha1- (a2-alpha2)
ifabs(delta_phi-eta* (a1-alpha1)) >svr.epsilon:
svr.alphas[i1] =a1
svr.alphas[i2] =a2
else:
finished=1
case1=1

elif (case2==0) and (alpha1>0or (alphas1==0anddelta_phi>2*svr.epsilon)) \
and (alphas2>0or (alpha2==0anddelta_phi>2*svr.epsilon)):
# 计算L和H (a1, a2*)
L=max(0, gamma)
H=min(svr.C, svr.C+gamma)

ifL<H:
ifeta>0:
a2=alphas2+ (delta_phi-2*svr.epsilon) /eta
a2=min(a2, H)
a2=max(L, a2)
a1=alpha1+a2-alphas2
else:
a2=L
a1=alpha1+ (a2-alphas2)
object1=-0.5*a1*a1*eta-2*svr.epsilon*a1+a1* (delta_phi+ (a1old-a1olds) *eta)
a2=H
a1=alpha1+ (a2-alpha2)
object2=-0.5*a1*a1*eta-2*svr.epsilon*a1+a1* (delta_phi+ (a1old-a1olds) *eta)
ifobject1>object2:
a2=L
else:
a2=H
a1=alpha1+ (a2-alphas2)

# 判断变化大小
ifabs(delta_phi-eta* (a1-alpha1)) >svr.epsilon:
svr.alphas[i1] =a1
svr.alphas_star[i2] =a2
else:
finished=1
case2=1


elif (case3==0) and (alphas1>0or (alpha1==0anddelta_phi<2*svr.epsilon)) \
and (alpha2>0or (alphas2==0anddelta_phi<2*svr.epsilon)):
# 计算L和H (a1*, a2)
L=max(0, -gamma)
H=min(svr.C, svr.C-gamma)

ifL<H:
ifeta>0:
a2=alpha2- (delta_phi-2*svr.epsilon) /eta
a2=min(a2, H)
a2=max(L, a2)
a1=alphas1+a2-alpha2
else:
a2=L
a1=alphas1+ (a2-alpha2)
object1=-0.5*a1*a1*eta-2*svr.epsilon*a1-a1* (delta_phi+ (a1old-a1olds) *eta)
a2=H
a1=alphas1+ (a2-alpha2)
object2=-0.5*a1*a1*eta-2*svr.epsilon*a1-a1* (delta_phi+ (a1old-a1olds) *eta)
ifobject1>object2:
a2=L
else:
a2=H
a1=alphas1+ (a2-alpha2)

# 判断变化大小
ifabs((delta_phi-eta* (-a1+alphas1))) >svr.epsilon:
svr.alphas_star[i1] =a1
svr.alphas[i2] =a2
else:
finished=1
case3=1


elif (case4==0) and (alphas1>0or (alpha1==0anddelta_phi<0)) \
and (alphas2>0or (alpha2==0anddelta_phi>0)):
# 计算L和H (a1*, a2*)
L=max(0, -gamma-svr.C)
H=min(svr.C, -gamma)

ifL<H:
ifeta>0:
a2=alphas2+delta_phi/eta
a2=min(a2, H)
a2=max(L, a2)
a1=alphas1-a2+alphas2
else:
a2=L
a1=alphas1- (a2-alphas2)
object1=-0.5*a1*a1*eta-a1* (delta_phi+ (a1old-a1olds) *eta)
a2=H
a1=alphas1- (a2-alphas2)
object2=-0.5*a1*a1*eta-a1* (delta_phi+ (a1old-a1olds) *eta)
ifobject1>object2:
a2=L
else:
a2=H
a1=alphas1- (a2-alphas2)
# 判断变化大小
ifabs((delta_phi-eta* (-a1+alphas1))) >svr.epsilon:
svr.alphas_star[i1] =a1
svr.alphas_star[i2] =a2
else:
finished=1
case4=1
else:
finished=1

delta_phi=delta_phi-eta* ((alpha1-alphas1) - (a1old-a1olds))


# 更新b值
b1, b2=0, 0
foriiinrange(len(svr.feature)):
b1+= (svr.alphas[ii] -svr.alphas_star[ii]) *svr.kernel(svr.feature[ii], x1)[0][0]
b2+= (svr.alphas[ii] -svr.alphas_star[ii]) *svr.kernel(svr.feature[ii], x2)[0][0]

b1=y1-b1
b2=y2-b2

b12= (b1+b2) /2
svr.b=b12


# 更新误差
svr.update_errors()

ifabs(delta_phi) >svr.epsilon:
return1
else:
return0

defexamineExample(svr, i2):
alpha2=svr.alphas[i2]
alphas2=svr.alphas_star[i2]
phi2=svr.errors[i2]

if (phi2>svr.epsilonandalphas2<svr.C) or (phi2<svr.epsilonandalphas2>0) \
or(-phi2>svr.epsilonandalpha2<svr.C) or (-phi2>svr.epsilonandalpha2>0):
manzu= [ioforioinrange(len(svr.feature)) ifsvr.alphas[io] !=0andsvr.alphas[io] !=svr.C]
iflen(manzu) >1:
delta_phi=0
i1=0
forhiinrange(len(svr.feature)):
# 选择最大的启发式选择
yphi=svr.errors[hi]
ifabs(yphi-phi2) >delta_phi:
delta_phi=abs(yphi-phi2)
i1=hi
iftakestep(svr, i1, i2):
return1
# 随机边界样本选择
forjjinrange(len(svr.feature)):
i1=np.random.choice(np.arange(len(svr.feature)), 1)[0]
ifsvr.alphas[i1] !=0andsvr.alphas[i1] !=svr.C:
iftakestep(svr, i1, i2):
return1
# 随机全部样本选择
forjjinrange(len(svr.feature)):
i1=np.random.choice(np.arange(len(svr.feature)), 1)[0]
iftakestep(svr, i1, i2):
return1
return0



# 主要函数
defmainfun(svr):
numChanged=0
examineAll=1
SigFig=-100
Loopcounter=0

while ((numChanged>0orexamineAll) orSigFig<3) andLoopcounter<10000:
Loopcounter+=1
numChanged=0
ifexamineAll:
print('全部样本')
foriginrange(len(svr.feature)):
numChanged+=examineExample(svr, ig)
else:
print('边界样本')
manzu= [ipoforipoinrange(len(svr.feature)) ifsvr.alphas[ipo] !=0andsvr.alphas[ipo] !=svr.C]
forgiinmanzu:
numChanged+=examineExample(svr, gi)

ifLoopcounter%2==0:
MinimumNumChanged=max(1, 0.1*len(svr.alphas))
else:
MinimumNumChanged=1

ifexamineAll==1:
examineAll=0
else:
ifnumChanged<MinimumNumChanged:
examineAll=1
dobject=0
pobject=0

forghiinrange(len(svr.feature)):
p1=svr.feature[ghi]
dobject+=max(0, svr.line_num(p1) -svr.labels[ghi] -svr.epsilon) * (svr.C-svr.alphas_star[ghi])
dobject-=min(0, svr.line_num(p1) -svr.labels[ghi] -svr.epsilon) *svr.alphas_star[ghi]
dobject+=min(0, svr.labels[ghi] -svr.epsilon-svr.line_num(p1)) * (svr.C-svr.alphas[ghi])
dobject-=max(0, svr.labels[ghi] -svr.epsilon-svr.line_num(p1)) *svr.alphas[ghi]


p1=svr.feature[ghi]
pobject+= (0.5* (svr.alphas[ghi] -svr.alphas_star[ghi]) * (svr.line_num(p1) -svr.b))
pobject-= (svr.epsilon* (svr.alphas[ghi] +svr.alphas_star[ghi]))
pobject+= (svr.labels[ghi] * (svr.alphas[ghi] -svr.alphas_star[ghi]))

print('gggggggggggggggggggggggggggggggggggg', pobject, dobject)
pobject+=dobject

SigFig=np.log10(dobject/ (abs(pobject) +1))
print('SigFig = %.9f'%SigFig)

print('结束训练')
returnsvr.alphas, svr.alphas_star, svr.b


# 返回训练数据、预测数据的输出值
defpl(trfe, trla, prfe, maxnum, minnum):
# 建立结构
svre=SVR(feature=trfe, labels=trla)
mainfun(svre)
traiout=svre.predictall(trfe)
preout=svre.predictall(prfe)

# 首先是数据转化范围
traiout=traiout* (maxnum-minnum) +minnum
preout=preout* (maxnum-minnum) +minnum

returntraiout, preout


# 绘图的函数
defhuitu(suout, shiout, c=['b', 'k'], sign='训练', cudu=3):
print(suout)
print(shiout)
# 绘制原始数据和预测数据的对比
plt.subplot(2, 1, 1)
plt.plot(list(range(len(suout))), suout, c=c[0], linewidth=cudu, label='%s：算法输出'%sign)
plt.plot(list(range(len(shiout))), shiout, c=c[1], linewidth=cudu, label='%s：实际值'%sign)
plt.legend()
plt.title('对比')

# 绘制误差和0的对比图
plt.subplot(2, 1, 2)
plt.plot(list(range(len(suout))), suout-shiout, c='r', linewidth=cudu, label='%s：误差'%sign)
plt.plot(list(range(len(suout))), list(np.zeros(len(suout))), c='r', linewidth=cudu, label='0值')
plt.legend()
plt.title('误差')
# 需要添加一个误差的分布图

# 显示
plt.show()


'''第四部分：最终的运行程序'''
if__name__=="__main__":
datasvr=rdata.model_data
outtri, poupre=pl(datasvr[0], datasvr[1].T[0], datasvr[2], datasvr[4][0], datasvr[4][1])

trii=datasvr[1].T[0] * (datasvr[4][0] -datasvr[4][1]) +datasvr[4][1]
huitu(outtri, trii, c=['b', 'k'], sign='训练', cudu=3)

prii=datasvr[3].T[0] * (datasvr[4][0] -datasvr[4][1]) +datasvr[4][1]
huitu(poupre, prii, c=['b', 'k'], sign='预测', cudu=3)