逻辑回归分类是机器学习中常用的一种分类方法,采用单位阶跃函数的特点进行分类。本文试图用理论推导、python实现来说明该算法。
单位阶跃函数
单位阶跃函数即Sigmoid
函数,其值范围在0~1
之间,并且在x=0处会产生跳跃。表达式如下:
g(z)=11+e−z
函数图形如下。该图形有个特点是,当x<<0时,f(x)=0; 当x>>0时,f(x)=1; 而f(0)=0.5。任何大于0.5的值都会被分为1类,而小于0.5的被归为0类。
由于Sigmoid函数是连续函数,反映的是一种
y=1
概率事件,如f(x)=0.7则表示该x最有可能(有70%的概率)被分为1类,有30%的概率被分为0类,可参考概率论中的极大似然估计。
寻找代价函数
逻辑回归的目的,是在训练集合寻找一个超平面将训练集进行分类,举下例的二元分类进行说明:
如上图所示,在平面中寻找一条线将数据分离。根据几何学可知该条线将数据集分为hθ>0与hθ<0 的两部分,令z=hθ 代入Sigmoid函数中,则完成了数据集的分类。其中hθ(x)
Andrew Ng在斯坦福课程中给出了逻辑回归的代价函数:
J(θ)=1m∑i=1mCost
Cost={−log(hθ(x))1−log(hθ(x))if y=1if y=0
Cost可写为
Cost=−ylog(hθ(x))−(1−y)(1−log(hθ(x))
该方法是根据概率论中的极大似然估计:当
y=1
时,hθ越大越接近1,则表示y=1
的概率越大,Cost越接近0,误差越小;反之对y=0
也是如此。因此只需要求得Cost的极小值就可以得到较好的分类效果
代入得到
J(θ)=−1m∑i=1m(ylog(hθ(x))−(1−y)log(1−hθ(x)))
用梯度下降求J(θ)的极小值
令∂J(θ)∂θ=0 , 由于公式复杂数据量大,在计算中通常不会采用直接求解,而采用梯度下降算法进行求解。J(θ)是关于θ的函数,梯度(也可以看成斜率)定义为∂J(θ)∂θ ,梯度计算为:
θ:=θ−α∂J(θ)∂θ
如下图示例,
α减小“斜率”使得
J(θ)向极小值靠近。
梯度公式推导
用向量表示方程:hθ(x)=g(θTx)
∂J(θ)∂θ=−1m∑i=1m∂(ylog(hθ(x))−+(1−y)log(1−hθ(x)))∂θ=−1m∑i=1my1hθ(x)∂hθ(x)∂θ−(1−y)11−hθ(x)∂hθ(x)∂θ=−1m∑i=1m(yg(θTx)−1−y1−g(θTx))∂∂θg(θTx)=−1m∑i=1m(yg(θTx)−1−y1−g(θTx))g(θTx)(1−g(θTx))∂θTx∂θ=−1m∑i=1m(y−g(θTx))x=−1m∑i=1m(y−hθ(x))x
则梯度下降公式为
θ:=θ+α1m∑i=1m(y−hθ(x))x=θ+α[(y−hθ(x))]x该式为矩阵表达式
求解过程中可忽略1m,将该项看成是α中的一项,即α=1mβ
Python实现
import numpy as np
def sigmoid(z):
return 1.0/(1 + np.exp(-z))
def gradDesc(dataSet, classes, alpha, cycles): #dataSet,classes分别为二维数组、一维数组
firstCol=np.ones((len(dataSet), 1))
trainData=np.concatenate((firstCol, dataSet), axis=1) #在第二维上拼接,即列方向拼接
m,n = trainData.shape
theta=np.ones(n)
for i in range(cycles):
h=sigmoid(np.dot(trainData, theta.T))
gap=classes-h
theta=theta + alpha*np.dot(gap.T, trainData)
return theta
测试算法
data=np.loadtxt(fname='mlslpic/logistci_data1.txt', dtype='float', delimiter=',')
dataSet=data[:,0:2]
classLabel=data[:,2]
weights=gradDesc(dataSet, classLabel,0.002, 1000000)
import matplotlib.pyplot as plt
x0=[]
x1=[]
for i in range(len(dataSet)):
if classLabel[i] == 0:
x0.append(dataSet[i])
else:
x1.append(dataSet[i])
x0=np.array(x0)
x1=np.array(x1)
plt.figure()
plt.plot(x0[:,0],x0[:,1],'r.')
plt.plot(x1[:,0],x1[:,1],'k+')
xp=np.arange(20, 110, 10)
yp=-(weights[0] + weights[1] * xp)/weights[2]
plt.plot(xp,yp)
plt.show()
分类结果
随机梯度下降
梯度下降算法每次更新系数 θ
- python实现
#随机梯度下降
import numpy as np
def randomGradDesc(dataSet, classLabel, cycles):
firstCol=np.ones((len(dataSet), 1))
trainData=np.concatenate((firstCol, dataSet), axis=1) #在第二维上拼接,即列方向拼接
m,n = trainData.shape
theta=np.ones(n)
for i in range(cycles):
for j in range(m):
alpha=0.001 + 4/(400.0 + i + j)
randIndex=int(np.random.uniform(0, m))
h=sigmoid(np.dot(trainData[randIndex], theta.T))
gap=classLabel[randIndex]-h
theta=theta + alpha*gap*trainData[randIndex]
np.delete(trainData,randIndex,axis=0)
return theta
与“批量”梯度下降的区别,每次取一组数据进行训练,计算的Sigmoid是数值而不是向量。
- 测试训练
data=np.loadtxt(fname='mlslpic/logistci_data1.txt', dtype='float', delimiter=',')
dataSet=data[:,0:2]
classLabel=data[:,2]
weights=randomGradDesc(dataSet, classLabel,20000)
import matplotlib.pyplot as plt
x0=[]
x1=[]
for i in range(len(dataSet)):
if classLabel[i] == 0:
x0.append(dataSet[i])
else:
x1.append(dataSet[i])
x0=np.array(x0)
x1=np.array(x1)
plt.figure()
plt.plot(x0[:,0],x0[:,1],'r.')
plt.plot(x1[:,0],x1[:,1],'k+')
xp=np.arange(20, 110, 10)
yp=-(weights[0] + weights[1] * xp)/weights[2]
plt.plot(xp,yp)
plt.show()
- 20000次迭代计算后的结果如下,对比梯度下降算法,分类效果差不多。但是所用的迭代次数及时间相对来说更少。