机器学习笔记（四）Logistic回归

我们都知道，如果预测值y是个连续的值，我们通常用回归的方法去预测，但如果预测值y是个离散的值，也就是所谓的分类问题，用线性回归肯定是不合理的，因为你预测的值没有一个合理的解释啊。比如对于二分类问题，我用{0，1}表示两类，那超出这个区间的预测值的意义是什么呢？所以我们就要引入一种特殊的回归方法，但通常用于分类问题，这就是Logistic回归，我们通过处理，将所有的预测值都限制在0到1以内，其形式如下：



这里要先介绍一下sigmoid函数，定义如下



它有一个特性如下



我们在Python中把它画出来看一下

import numpy as np
from matplotlib import pyplot as plt

def sigmoid(x):
return 1/(1+np.exp(-x))

x=np.linspace(-6,6,100)
y=sigmoid(x)
plt.plot(x,y,'r-')
plt.title('Sigmoid')
plt.show()

从图中可以看出sigmoid函数取值在[0,1]之间，适合做一个概率分布函数。那么我们的目标就是选择一组合适的参数theta，满足样本属于正类时thet*x尽可能大，样本属于负类时，theta*x就尽可能小，这样通过sigmoid函数，我们就可以确定该样本属于正类的概率有多大了，那么有



那么每个样本为正类的概率就是



假设有m个样本，似然函数为



对数似然函数为



我们要让似然函数取最大值，那么可以定义Logistic回归的损失函数为-l(theta)，那么就转变成求损失函数的最小值了，我们可以用线性回归中所提到的梯度下降


巧合的是我们发现这个梯度的形式和线性回归完全是一样的，我们同样可以利用批量梯度下降或者随机梯度下降来求这个最值。

接下来我们构建一个实例来运用Logistic回归。

首先是构造数据。我们以y=x为边界，曲线上方称为正例，曲线下方称为负例，区域范围为[0,5]

import numpy as np
from matplotlib import pyplot as plt

np.random.seed(0)
data=[]
for i in np.arange(100):
a=np.random.rand()*5
b=np.random.rand()*5
if a>b:
data.append([a,b,1])
else:
data.append([a,b,0])
data=np.array(data)
posdata=data[data[:,2]==1]
negdata=data[data[:,2]==0]
print '正例数量为',len(posdata),'负例数量为',len(negdata)
plt.plot(posdata[:,0],posdata[:,1],'r^',label='Pos')
plt.plot(negdata[:,0],negdata[:,1],'g*',label='Neg')
plt.legend(loc='upper right')
plt.title('Row data')
plt.grid()
plt.show()

正例数量为 54 负例数量为 46



好像正例负例和我之前说的弄反了，但无所谓啦，这个不影响我们测试的……

我们首先用sklearn库里自带的Logistic回归试一下，模型训练完成之后，我们在区域内取500*500个点来进行预测，看结果是否满足我们的要求

log=LogisticRegression()
model=log.fit(x,y)
print model.coef_,model.intercept_
N, M = 500, 500
x1_min, x1_max = x[:, 0].min(), x[:, 0].max()
x2_min, x2_max = x[:, 1].min(), x[:, 1].max()
t1 = np.linspace(x1_min, x1_max, N)
t2 = np.linspace(x2_min, x2_max, M)
x1, x2 = np.meshgrid(t1, t2)
x_test = np.stack((x1.flat, x2.flat), axis=1)
y_hat=model.predict(x_test).reshape(x1.shape)
plt.pcolormesh(x1, x2, y_hat, cmap= mpl.colors.ListedColormap(['r', 'g']))
plt.show()

[[ 2.32453508 -2.39690698]] [ 0.0221314]



从图可以看到，预测的结果非常好，这和我们一开始造数据的规则十分符合，说明Logistic回归确实有很好的分类效果。

接下来我们再自己用梯度下降法来实现以下看是否能达到同样的效果。

def sigmoid(x):
return 1/(1+np.exp(-x))
alpha=0.5
theta=np.array([0,0,0])
x0=np.ones(x.shape[0]).reshape([-1,1])
x=np.c_[x0,x]
for i in xrange(5000):
y_hat=sigmoid(np.dot(x,theta.T))
grad0=(y-y_hat).reshape(-1,1)*x0
grad1=(y-y_hat).reshape(-1,1)*x[:,1].reshape(-1,1)
grad2=(y-y_hat).reshape(-1,1)*x[:,2].reshape(-1,1)
theta[0]=theta[0]+alpha*sum(grad0)
theta[1]=theta[1]+alpha*sum(grad1)
theta[2]=theta[2]+alpha*sum(grad2)
print theta
x0test=np.ones(x_test.shape[0]).reshape([-1,1])
x_test=np.c_[x0test,x_test]
y_hat=sigmoid(np.dot(x_test,theta.T))
y_hat[y_hat>0.5]=1
y_hat[y_hat<=0.5]=0
y_hat=y_hat.reshape(x1.shape)

plt.pcolormesh(x1, x2, y_hat, cmap= mpl.colors.ListedColormap(['r', 'g']))
plt.show()

系数为[  0  68 -69]



分类效果同样十分卓著，甚至感觉比sklearn还要好一些……不过有个明显的差异就是拟合系数的不同sklearn的结果[ 0.0221314 2.32453508 -2.39690698]，自己写的梯度下降算法的结果[ 0 68-69]。差别好像还挺大的，不过参数间的比例关系倒是挺一致的，这个地方我就有点懵逼了，虽然两种结果都挺不错的，但为何最后参数相差这么大，这个问题我还要好好思考一下，如果有看到的同学想明白了希望可以指点我一下。

最后就是Logistic回归的推广啦。前面我们讨论的都是二分类问题，那么针对多分类的问题如何处理呢？

感觉最直观的想法就是针对每个类做一个一对多的Logistic回归，然后每个样本在所有的模型里都算一下，将其归于得分最高的那一类。

当然，我们也可以直接用Logistic回归的推广，即softmax回归



然后接下来推导的思路和上面完全一样，求似然函数，对数似然函数，损失函数为对数似然函数的负值，最后利用梯度下降求极值，寻求最佳参数组合。

在实际的运用中，特别是一些封装好的库中，它把Logistic和softmax回归封装在了一起，它可以自动根据你分类的数量决定用哪一种回归，所以完全不用你多操心。下次再写个demo演示一下softmax回归吧，今天就到这儿，打这些推导公式心好累，周末要结束了，唉/(ㄒoㄒ)/~~