机器学习必备算法之(一)逻辑回归(logistics regression)及Python实现

笔者为数学系的一个小白，最近系统的在复习机器学习以及一些深度学习的内容，准备开个博记录一下这个有趣又痛苦的过程~hiahiahia，主要记录机器学习的几大经典算法的理论以及Python的实现。非计算机的大神！！！代码不会太牛(求轻喷，暴风哭泣)。话不多说，就从逻辑回归开始吧。

线性回归

首先我们来说说什么是线性回归。给定一个数据集D={(x1,y1),…,(xm,ym)}D=\{(x_1,y_1),\dots,(x_m,y_m)\}D={(x1​,y1​),…,(xm​,ym​)}，其中xi=(xi1;…;xid},yi∈Rx_i=(x_{i1};\dots;x_{id}\},y_i\in Rxi​=(xi1​;…;xid​},yi​∈R。线性回归是试图寻找一个线性模型尽可能的准确预测输出标记yyy。我们可以把模型归结为

f(xi)=wxi+bif(x_i)=wx_i+b_i f(xi​)=wxi​+bi​

尽可能地使得f(xi)≈yif(x_i)\approx y_if(xi​)≈yi​。那么如何如何让描述这个近似相等的程度呢，这里有很多种数学意义上的“距离”可以公式化的表述，我们以均方误差为例，把上述线性模型问题写成一个标准的优化问题(实际上就是最小二乘问题啦)：

(w∗,b∗)=argminw,b∑i=1m(f(xi)−yi)2(w^*,b^*)=argmin_{w,b}{\sum^m_{i=1}(f(x_i)-y_i)^2} (w∗,b∗)=argminw,b​i=1∑m​(f(xi​)−yi​)2

为了方便理解我们举个栗子，(如下图所示,图片来源于网络)，假设在二维空间中，数据集就是平面上的众多点的集合，我们想找到一条直线尽可能多的穿过这些点，直线上的点到直线距离为0，但总会有一些点不在这条直线上，与直线有一定的距离，我们的目的就是寻找到最优的那对(w∗,b∗)(w^*,b^*)(w∗,b∗)，使得这个距离的平方和最小。

推广到多元线性回归上，我们可以得到下式的优化问题，这里不推显示解啦，求导等于0及可。

w^∗=arg min⁡w^(y−Xw^)T(y−Xw^)\hat{w}^* =\argmin_{\hat{w}}{(y-X\hat{w})^T(y-X\hat{w})} w^∗=w^argmin​(y−Xw^)T(y−Xw^)

其中w^=(w;b)\hat{w}=(w;b)w^=(w;b)，X=(X;1)X=(X;1)X=(X;1)

逻辑回归原理

接下来我们来说回我们的主题，刚刚的线性回归得到的是连续的预测标记的值，是RRR上的实值。那么我们如果想做分类问题怎么办呢？

也就是说，我们想让这个输出变成离散的，比如输出一个病人通过几项指标最后预测出患病或者是健康，再比如我们通过一些天气条件判断明天下雪或者不下雪…一个不难想到的办法就是进行一个阈值的判断。我们得到一个预测后是实数，如果这个预测值大于阈值，我们把它归结在一类中，如果小于就放在另一类中，这就是简单的二分类。二分类任务中，最理想的比较预测后数值与阈值的函数就是阶跃函数(unit-step function)：

y={0z<00.5z=01z>0y=\left\{ \begin{array}{rcl} 0& & {z<0}\\ 0.5 & & {z=0}\\ 1 & & {z>0}\\ \end{array} \right. y=⎩⎨⎧​00.51​​z<0z=0z>0​

但阶跃函数不好用的一点就在于其在0点处不连续不可微。那么用来替换阶跃函数的一个函数就是一种sigmoid函数，逻辑函数(logistics function)（当当~主角终于登场了）：

y=11+e−zy=\frac{1}{1+e^{-z}} y=1+e−z1​

通过这个函数我们可以把预测出来的实值映射到[0,1]之间。并且它有一个非常重要的意义：其输出的值代表着预测为正类的概率。为什么呢？我们来看看下图，因为比如我们的得到的数值越接近于1，也就是zzz越接近于∞\infty∞，即使阈值向右移动（增大），它也一定大于阈值，它还是正类，也就是说，它成为正类的可能性更大。反之同理。

所以说，逻辑回归叫做回归，但它并不是回归算法，它是一种分类算法。

对二分类问题的最大似然估计

根据上面我们分析得到，输出f(xi)f(x_i)f(xi​)表现的是预测为正类的概率，则有：

p(y=1∣x)=ew^Tx1+ew^Tx=f(xi)p(y=1|x)=\frac{e^{\hat{w}^Tx}}{1+e^{\hat{w}^Tx}}=f(x_i) p(y=1∣x)=1+ew^Txew^Tx​=f(xi​)

p(y=0∣x)=11+ew^Tx=1−f(xi)p(y=0|x)=\frac{1}{1+e^{\hat{w}^Tx}}=1-f(x_i) p(y=0∣x)=1+ew^Tx1​=1−f(xi​)

这里的w^=(w;b)\hat{w}=(w;b)w^=(w;b)。

（接下来是求极大似然估计，忘了的同学们可以去查阅一下茆时松那本概率论与数理统计，还记得本科上课的那个老师每节课都用塑料普通话反复强调要好好学茆si松 的那本书hahahahaha跑题了，咳咳）

对于给定的含有m个样本点的数据集{(xi,yi)}i=1m\{(x_i,y_i)\}^m_{i=1}{(xi​,yi​)}i=1m​，我们求其似然函数：

L(w^)=Πi=1mP(yi∣xi;w^)=Πi=1mf(xi)yi(1−f(xi))1−yiL(\hat{w})=\Pi^{m}_{i=1}{P(y_i|x_i;\hat{w})}=\Pi^{m}_{i=1}{f(x_i)^{y_i}(1-f(x_i))^{1-y_i}} L(w^)=Πi=1m​P(yi​∣xi​;w^)=Πi=1m​f(xi​)yi​(1−f(xi​))1−yi​

对其两边同时求对数，得到对数似然函数：

logL(w^)=∑i=1m(yilog(f(xi))+(1−yi)log(1−f(xi)))=∑i=1m(yilog(f(xi)1−f(xi))+log(1−f(xi)))=∑i=1m(yiw^xi−log(1+ew^Txi))\begin{aligned} logL(\hat{w}) &=\sum^m_{i=1}(y_ilog(f(x_i))+(1-y_i)log(1-f(x_i)))\\ &=\sum^m_{i=1}(y_ilog(\frac{f(x_i)}{1-f(x_i)})+log(1-f(x_i)))\\ &=\sum^m_{i=1}(y_i\hat{w}x_i-log(1+e^{\hat{w}^Tx_i}))\\ \end{aligned} logL(w^)​=i=1∑m​(yi​log(f(xi​))+(1−yi​)log(1−f(xi​)))=i=1∑m​(yi​log(1−f(xi​)f(xi​)​)+log(1−f(xi​)))=i=1∑m​(yi​w^xi​−log(1+ew^Txi​))​

接下来求其极大即可。常用梯度上升法和牛顿法，更新规则为：

梯度上升法：

w^j+1=w^j+α∂logL(w^)∂w^∂logL(w^)∂w^=∑i=1mxi(yi−f(xi))\begin{aligned} \hat{w}_{j+1}&=\hat{w}_j+\alpha \frac{\partial log L(\hat{w})}{\partial \hat{w}}\\ \frac{\partial log L(\hat{w})}{\partial \hat{w}}&=\sum_{i=1}^mx_i(y_i-f(x_i)) \end{aligned} w^j+1​∂w^∂logL(w^)​​=w^j​+α∂w^∂logL(w^)​=i=1∑m​xi​(yi​−f(xi​))​

牛顿法：

w^j+1=w^j−∂2logL(w^)∂w^∂w^T∂logL(w^)∂w^∂logL(w^)∂w^=∑i=1mxi(yi−f(xi))∂2logL(w^)∂w^∂w^T=∑i=1mxixiTf(xi)(f(xi)−1)\begin{aligned} \hat{w}_{j+1}&=\hat{w}_j-\frac{\partial^2 log L(\hat{w})}{\partial \hat{w} \partial \hat{w}^T} \frac{\partial log L(\hat{w})}{\partial \hat{w}}\\ \frac{\partial log L(\hat{w})}{\partial \hat{w}}&=\sum_{i=1}^mx_i(y_i-f(x_i))\\ \frac{\partial^2 log L(\hat{w})}{\partial \hat{w} \partial \hat{w}^T}&=\sum^m_{i=1}x_ix_i^T f(x_i)(f(x_i)-1) \end{aligned} w^j+1​∂w^∂logL(w^)​∂w^∂w^T∂2logL(w^)​​=w^j​−∂w^∂w^T∂2logL(w^)​∂w^∂logL(w^)​=i=1∑m​xi​(yi​−f(xi​))=i=1∑m​xi​xiT​f(xi​)(f(xi​)−1)​

交叉熵损失函数

Loss(f(xi),yi)=−1m∑i=1mlogP(yi∣f(xi))=−1m∑i=1mlog(f(xi)yi(1−f(xi)1−yi)=−1m∑i=1m(yilog(f(xi))+(1−yi)log(1−f(xi)))\begin{aligned} Loss(f(x_i),y_i) &=-\frac{1}{m}\sum_{i=1}^mlogP(y_i|f(x_i))\\ &=-\frac{1}{m}\sum_{i=1}^mlog(f(x_i)^{y_i}(1-f(x_i)^{1-y_i})\\ &=-\frac{1}{m}\sum_{i=1}^m( y_i log(f(x_i))+(1-y_i)log({1-f(x_i)})) \end{aligned} Loss(f(xi​),yi​)​=−m1​i=1∑m​logP(yi​∣f(xi​))=−m1​i=1∑m​log(f(xi​)yi​(1−f(xi​)1−yi​)=−m1​i=1∑m​(yi​log(f(xi​))+(1−yi​)log(1−f(xi​)))​

其他问题

如何让实现多分类

各类别互斥：用softmax输出对各类别的预测概率无明显互斥或交叉：搭建多个二分类器OVR最后结果取max无明显互斥或交叉：搭建多对二分类器OVO最后结果取max

解决类别不平衡问题

人工该表数据分布使其平衡：上采样、下采样阈值移动

Python实现

这里给出一个基于梯度上升法的简单python，暂时先不用sklearn啦~~~

先给出一个思路：先初始化回归系数（一般初始化为1），然后计算数据集的梯度

再使用alpha*gradient更新回归系数，反复循环最终得到回归系数，然后得到可视化结果或分类器。

（一）载入数据：我的这个小数据们放在这里啦链接：data 提取码：033c

filename='data.txt'def loadDataSet(): #读取数据（这里只有两个特征）x_data = []y_data = []fr = open(filename)for line in fr.readlines():lineArr = line.strip().split()x_data.append([1.0, float(lineArr[0]), float(lineArr[1])]) #前面的1，表示方程的常量。比如两个特征X1,X2，共需要三个参数，W1+W2*X1+W3*X2y_data.append(int(lineArr[2]))return x_data,y_data#得到的是list类型

（二）定义sigmoid函数

#定义逻辑函数def Sigmoid(z):return 1.0/(1+np.exp(-z))

（三）定义梯度上升法(或者用随机梯度法都可以啦~)

#定义梯度上升法 x_data=(1;x),y_data是标记的行向量def GradUp(x_data,y_data,iternum):x=np.mat(x_data)y=np.mat(y_data).transpose()m,n=x.shapealpha=0.001weight=np.ones((n,1))#初始化权重wfor k in range(iternum):h=Sigmoid(x*weight)error=y-hweight+=alpha*x.transpose()*errorreturn weight

（四）可视化

这里解释一下为什么画图画的是 y = (-weight[0] - weight[1] * x) / weight[2]：根据sigmoid函数知道，当z=w^Tx>0z=\hat{w}^Tx>0z=w^Tx>0时是正类，所以分界线为超平面w^Tx=0,w^=(b;w)\hat{w}^Tx=0,\hat{w}=(b;w)w^Tx=0,w^=(b;w)。注：纵坐标可不是什么label啊！

#可视化处理画出决策def plotFit(weight):x_data,y_data=loadDataSet()x_array=np.array(x_data)# print(x_array.shape)n = x_array.shape[0]xcord1=[];ycord1=[]xcord2=[];ycord2=[]for i in range(n):if int(y_data[i])==1:xcord1.append(x_array[i,1])ycord1.append(x_array[i,2])else:xcord2.append(x_array[i,1])ycord2.append(x_array[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=np.arange(-3,3,0.1)y = (-weight[0] - weight[1] * x) / weight[2]#超平面w1*x1+w2*x2+b=0分割ax.plot(x,y)plt.xlabel('X1'); plt.ylabel('X2');plt.show()

（五）测试啦：

x_data,y_data=loadDataSet()weight=GradUp(x_data,y_data,200)plotFit(weight)print(weight)

结果：

输出回归系数：

[[ 2.77955136][ 0.36268202][-0.44923697]]

等有以后复习到sklearn的时候再补上这部分代码吧。目前这个代码只写了可视化的一个结果， 如果想得到分类结果，要将数据及回归系数带入sigmoid函数得到预测值y，让y与0.5比较然后得到分类1/0。