In [127]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
In [166]:
def kernelTrans(X, A, kTup):
'''
Parameters:
X: np.array
A: a row of np.array
kTup: a tuple like ('rbf', 0.4)
return:
K: np.array(a column)
'''
m, n = X.shape
K = np.zeros((m, ))
if kTup[0] == 'lin':
K = X*A.T
elif kTup[0] == 'rbf':
for j in range(m):
deltaRow = X[j] - A
K[j] = np.dot(deltaRow, deltaRow.T)
K = np.exp(K/(-1*kTup[1]**2))
else:
raise NameError('Houston We hace a problem -- That kernel is not recognized')
return K
In [129]:
def clipAlpha(alpha, High, Low):
if alpha > High:
alpha = High
elif Low > alpha:
alpha = Low
return alpha
In [130]:
class optStruct:
def __init__(self, data, labels, C, toler, kTup):
self.X = data
self.labels = labels
self.C = C
self.tol = toler
self.m = data.shape[0]
self.alphas = np.zeros((self.m, ))
self.b = 0
self.errCache = np.zeros((self.m, 2)) # store error cache
self.K = np.zeros((self.m, self.m)) # store the kernel function's values
for i in range(self.m):
self.K[:, i] = kernelTrans(self.X, self.X[i], kTup)
In [131]:
def calculateEk(oS, k):
fXk = np.dot(oS.alphas*oS.labels, oS.K[:, k])+ oS.b
Ek = fXk - float(oS.labels[k])
return Ek
In [132]:
def select_random(index, m):
j = index
while j == index:
j = int(np.random.uniform(0, m))
return j
In [133]:
def updateEk(oS, k):
Ek = calculateEk(oS, k)
oS.errCache[k] = [1, Ek]
In [134]:
def selectJ(i, oS, Ei):
maxK = -1
maxDeltaE = 0
Ej = 0
oS.errCache[i] = [1, Ei]
validEcacheList = np.nonzero(oS.errCache[:, 0])[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList:
if k == i:
continue
Ek = calculateEk(oS, k)
deltaE = abs(Ei - Ek)
# pick up the max one
if (deltaE > maxDeltaE):
maxK = k
maxDeltaE = deltaE
Ej = Ek
return maxK, Ej
else:
# This is for the first time to choose an alpha
j = select_random(i, oS.m)
Ej = calculateEk(oS, j)
return j, Ej
In [135]:
def innerLoop(i, oS):
'''
return:
0: no alpha Pairs are Changed
1: Changed
'''
Ei = calculateEk(oS, i)
if ((oS.labels[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or\
((oS.labels[i]*Ei > -oS.tol) and (oS.alphas[i] > 0)):
j, Ej = selectJ(i, oS, Ei) # different from simpleSMO
alphaIold = oS.alphas[i].copy()
alphaJold = oS.alphas[j].copy()
if (oS.labels[i] != oS.labels[j]):
Low = max(0, oS.alphas[j] - oS.alphas[i])
High = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
Low = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
High = min(oS.C, oS.alphas[j] + oS.alphas[i])
if Low == High:
print "L==H"
return 0
eta = 2.0 * oS.K[i, j] - oS.K[i, i] - oS.K[j, j]
if eta >= 0:
print"eta>=0"
return 0
oS.alphas[j] -= oS.labels[j]*(Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j], High, Low)
updateEk(oS, j)
if (abs(oS.alphas[j] - alphaJold) < 0.00001):
print "j not moving enough"
return 0
oS.alphas[i] += oS.labels[j]*oS.labels[i]*(alphaJold - oS.alphas[j])
updateEk(oS, i)
b1 = oS.b - Ei - oS.labels[i]*(oS.alphas[i] - alphaIold)*oS.K[i, i] -\
oS.labels[j]*(oS.alphas[j] - alphaJold)*oS.K[i, j]
b2 = oS.b - Ej - oS.labels[i]*(oS.alphas[i] - alphaIold)*oS.K[i, j] -\
oS.labels[j]*(oS.alphas[j] - alphaJold)*oS.K[j, j]
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[j]):
oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]):
oS.b = b2
else:
oS.b = (b1 + b2)/2.0
return 1
else:
return 0
In [167]:
def SMO(data, labels, C, toler, maxIter, kTup):
'''
returns:
b, alphas
'''
oS = optStruct(data, labels, C, toler, kTup) # kernel matrix calculated here
iteration = 0
entireSet = True
alphaPairsChanged = 0
while (iteration < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged = 0
if entireSet:
for i in range(oS.m):
alphaPairsChanged += innerLoop(i, oS)
print "fullset, iter: %d, i: %d, pairs changed: %d" %(iteration, i , alphaPairsChanged)
iteration += 1
else:
nonBoundIs = np.nonzero((oS.alphas > 0) * (oS.alphas < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerLoop(i ,oS)
print "non-bound, iter: %d, i: %d, pairs changed: %d" %(iteration, i , alphaPairsChanged)
iteration += 1
if entireSet:
entireSet = False
elif (alphaPairsChanged == 0):
entireSet = True
print "Iteraion number: %d" %iteration
return oS.b, oS.alphas
In [157]:
def testRBF(k1=1.3):
# read data using pandas
df = pd.read_csv('testSetRBF.txt', delimiter='\t', names=['x1', 'x2', 'y'])
data = df.loc[:, ['x1', 'x2']].values
labels = df.loc[:, 'y'].values
# using SMO with racial bias function to obtain b and alphas
b, alphas = SMO(data, labels, 200, 0.0001, 10000, ('rbf', k1))
# find support vectors
svInd = np.nonzero(alphas > 0)[0]
sVs = data[svInd]
labelSV = labels[svInd]
print "there are %d Support Vectors" % sVs.shape[0]
# calculate the precision of training set
m, n = data.shape
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, data[i], ('rbf', k1))
predict = np.dot(kernelEval.T, labelSV*alphas[svInd]) + b
if np.sign(predict) != np.sign(labels[i]):
errorCount += 1
print "the training error rate is: %f" %(float(errorCount)/m)
# read the test data and calculate the precision of test data
df = pd.read_csv('testSetRBF2.txt', delimiter='\t', names=['x1', 'x2', 'y'])
data = df.loc[:, ['x1', 'x2']].values
labels = df.loc[:, 'y'].values
errorCount = 0
for i in range(m):
kernelEval = kernelTrans(sVs, data[i], ('rbf', k1))
predict = np.dot(kernelEval, labelSV*alphas[svInd]) + b
if np.sign(predict) != np.sign(labels[i]):
errorCount += 1
print "the test error rate is: %f" % (float(errorCount)/m)
In [158]:
K = np.zeros((100,100))
In [159]:
df = pd.read_csv('testSetRBF.txt', delimiter='\t', names=['x1', 'x2', 'y'])
data = df.loc[:, ['x1', 'x2']].values
labels = df.loc[:, 'y'].values
In [165]:
testRBF(k1=0.4)
In [ ]: