假设我们有一些数据点,我们使用一条直线对这些点进行拟合,这条线称为最佳拟合直线,这个拟合过程称为回归。利用Logistic回归进行分类的主要思想是:根据现有数据对分类边界线建立回归公式,以此进行分类。我们想要得到一个函数,能够接受所有的输入然后预测出类别。例如在两个类的情况下,函数输出0或1。该函数称为海维塞德阶跃函数(Heaviside step function),或者直接称为单位阶跃函数。但是海维塞德阶跃函数的问题在于:该函数在跳跃点上从0瞬间跳跃到1,这个瞬间跳跃过程有时很难处理。不过另外一个函数也有类似的性质,而且数学上更容易处理,它就是Sigmoid函数,具体公式如下:
σ(z)=11+e −z
确定了分类器的函数形式后,就需要确定最佳回归系数了。
Sigmoid函数输入记为z,公式为:
z=w 0 x 0 +w 1 x 1 +w 2 x 2 +⋯+w n x n
也可以写成向量的形式:
Z=W T X
其中x是分类器的输入数据,向量w就是要找到的最佳系数。
梯度上升
基本思想:要找到某函数的最大值,最好的方法是沿着该函数的梯度方向探寻。如果梯度记为
∇
,则函数f(x,y)的梯度由下式表示:
∇f(x,y)=⎛ ⎝ ⎜ ∂f(x,y)∂x ∂f(x,y)∂y ⎞ ⎠ ⎟
这个梯度意味着要沿着x的方向移动 ∂f(x,y)∂x ,沿着y的方向移动 ∂f(x,y)∂y
梯度算子总是指向函数值增长最快的方向,移动量大小称为步长,记为 α ,梯度上升算法的迭代公式如下:
w:=w+α∇ w f(w)
它将一直被迭代执行,直到达到某个停止条件,比如迭代次数达到某个特定值,或者算法某个可允许的误差范围内。下面使用100个样本点,每个点包含两个数值型特征,我们将通过梯度上升找到最佳回归系数,也就是拟合出Logistic回归模型的最佳参数。
from numpy import *
def loadDataSet():
dataMat = []; labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
def sigmoid(inX):
return 1.0/(1+exp(-inX))
def gradAscent(dataMatIn, classLabels):
dataMatrix = mat(dataMatIn) #convert to NumPy matrix
labelMat = mat(classLabels).transpose() #convert to NumPy matrix
m,n = shape(dataMatrix)
alpha = 0.001
maxCycles = 500
weights = ones((n,1))
for k in range(maxCycles): #heavy on matrix operations
h = sigmoid(dataMatrix*weights) #matrix mult
error = (labelMat - h) #vector subtraction
weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
return weights
通过以下代码运行:
dataArr,labelMat=loadDataSet()
print gradAscent(dataArr, labelMat)
得到回归系数如下:
[[ 4.12414349]
[ 0.48007329]
[-0.6168482 ]]
画出决策边界:
def plotBestFit(weights):
import matplotlib.pyplot as plt
dataMat,labelMat=loadDataSet()
dataArr = array(dataMat)
n = shape(dataArr)[0]
xcord1 = []; ycord1 = []
xcord2 = []; ycord2 = []
for i in range(n):
if int(labelMat[i])== 1:
xcord1.append(dataArr[i,1]); ycord1.append(dataArr[i,2])
else:
xcord2.append(dataArr[i,1]); ycord2.append(dataArr[i,2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1, ycord1, s=30, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X1'); plt.ylabel('X2');
plt.show()
通过以下代码运行:
weights = gradAscent(dataArr,labelMat)
plotBestFit(weights.getA())
结果如下:
分类结果还不错,但是需要大量的计算,所以用以下算法改进:
def stocGradAscent0(dataMatrix, classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
weights = ones(n) #initialize to all ones
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
weights = weights + alpha * error * dataMatrix[i]
return weights
运行:
weights = stocGradAscent0(array(dataArr), labelMat)
plotBestFit(weights)
结果如下:
结果不像之前那样好,是因为这次迭代次数减少了。
随机梯度算法改进。
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = shape(dataMatrix)
weights = ones(n) #initialize to all ones
for j in range(numIter):
dataIndex = range(m)
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001 #apha decreases with iteration, does not
randIndex = int(random.uniform(0,len(dataIndex)))#go to 0 because of the constant
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
运行:
weights = stocGradAscent1(array(dataArr), labelMat)
plotBestFit(weights)
结果如下: