[机器学习][源码]机器学习实战ch6 SVM

把代码保存于此，python3实现，详解就参考《机器学习实战》（Peter Harrington）啦...

1. svmMLiA.py : 

'''ref:
http://blog.csdn.net/puqutogether/article/details/44565093
'''
from numpy import *
# 6-1
def loadDataSet(filename):
dataMat=[];labelMat=[]
fr=open(filename)
for line in fr.readlines():
lineArr=line.strip().split('\t')
dataMat.append( [float(lineArr[0]),float(lineArr[1])] )
labelMat.append( float(lineArr[2]) )
return dataMat,labelMat


#随机选择随机数，不等于i
def selectJrand(i,m):
j=i
while(j==i):
j=int(random.uniform(0,m)) # 一直在挑选随机数j，直到不等于i，随机数的范围在0~m
return j


#门限函数, 最大不能超过H，最小不能低于L
def clipAlpha(aj,H,L):
if aj>H:
aj=H
if aj<L:
aj=L
return aj


# 6-2
#简化版的SMO函数
#随机选取第2个alpha
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
# 输入数据，标记，常数C，容错率，最大迭代次数
dataMatrix = mat(dataMatIn)  # 转换成矩阵
labelMat = mat(classLabels).transpose() # 转换成矩阵，并转置，标记成为一个列向量，每一行和数据矩阵对应
m,n = shape(dataMatrix) # 行(样本数)，列 
b = 0 # 参数b的初始化
alphas = mat( zeros((m,1)) )# 参数alphas是个list，初始化也是全0，大小等于样本数
iter = 0 # 当前迭代次数，maxIter是最大迭代次数 


while(iter<maxIter):
alphaPairsChanged=0 # 表示alpha对是否更新，=0 表示未更新
for i in range(m): # 【外循环】，选第一个alpha
# 第i样本的预测类别, f(xi)=w xi+b
# multiply(a,b)就是个乘法,如果a,b是两个长度为N的数组,那么对应元素相乘,得到的仍是一个长度为N的数组 
# dataMatrix: N*p, 
# multiply(alphas,labelMat).T: 1*N
# dataMatrix*dataMatrix[i,:].T: N*1
fXi = float( multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[i,:].T) )  +  b
Ei = fXi - float(labelMat[i])  # 预测的误差


#是否可以继续优化--alpha(i)违反KKT条件，继续优化
if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or \
((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
j=selectJrand(i,m) # 【内循环】，选择第j个样本
# 第j个样本的预测类别
fXj=float( multiply(alphas,labelMat).T * (dataMatrix*dataMatrix[j,:].T) )  +  b
Ej = fXj - float(labelMat[j])  # 误差


alphaIold = alphas[i].copy() # 拷贝,分配新的内存,旧的alpha[i]的值
alphaJold = alphas[j].copy() # 拷贝,分配新的内存,旧的alpha[j]的值


#更新L和H, alpha(j)的范围
if(labelMat[i]!=labelMat[j]):
L=max(0,alphas[j]-alphas[i])
H=min(C,C+alphas[j]-alphas[i])
else:
L=max(0,alphas[j]+alphas[i]-C)
H=min(C,alphas[j]+alphas[i])
if L==H:
print("L=H")
continue  #结束当前循环


#更新alpha(j)
eta=2.0*dataMatrix[i,:]*dataMatrix[j,:].T - \
dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
if eta>=0: 
print("eta>=0")
continue  #结束当前循环
alphas[j] -= labelMat[j]*(Ei-Ej)/eta
alphas[j] = clipAlpha(alphas[j],H,L) # 门限函数将alpha(j)调整到[L,H]之间


if(abs(alphas[j]-alphaJold)<0.00001): #如果alpha(j)修改量很微小
print("j not moving enough")
continue  #结束当前循环


alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j]) #更新alpha[i]


# 重新计算b
b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
if (0 < alphas[i]) and (C > alphas[i]): b = b1
elif (0 < alphas[j]) and (C > alphas[j]): b = b2
else: b = (b1 + b2)/2.0


alphaPairsChanged += 1 # 发生改变的 alpha对 的数目
print("iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsCh