逻辑回归 二/多分类 +Python代码实现

逻辑回归 二/多分类 +Python代码实现

前言：

​ 本文主要讲解逻辑回归的二分类和多分类问题。在此，笔者将不做过多引入的讲解，而是直接从理论出发。因为在我看来，网上对这方面的讲解已经很多，写得也非常的好。笔者主要是讲解公式转化为代码这一部分。

关键字：逻辑回归、sigmoid、softmax、Python

sigmoid:

p = 1 1 + e − w T x p=frac{1}{1+e^{-w^Tx}} p=1+e−wTx1​

​ 对于二分类问题，我们使用的是sigmoid函数，采用极大似然估计的方法进行估计。首先二分类的似然函数学过概率论的人都知道，他是这样的
f ( x ) = p y i ( 1 − p ) ( 1 − y i ) f(x)=p^{y_{i}}(1 - p)^{(1 - y_{i})} f(x)=pyi​(1−p)(1−yi​)
​ 为什么是极大似然估计？因为极大似然能够估计出样本的分布空间。如果我们的样本具有代表性，那么理论上此时的样本的模型参数也是接近总体的样本的模型参数。

损失函数及其梯度计算：

​ 显而易见，我们可以直接将极大似然估计作为损失函数，当燃了，损失函数一般我们是越小越好。而极大似然估计是越大越好，so？我们要对极大似然取反，最终损失函数应当如下：
L ( w ) = min ⁡ w − 1 m ∏ i = 1 n p y i ( 1 − p ) 1 − y i L(w)=minlimits_{w}-frac{1}{m}{prodlimits_{i=1}^n{p^{y_i}(1-p)^{1-y_i}}} L(w)=wmin​−m1​i=1∏n​pyi​(1−p)1−yi​
​ 为了便于运算，一般情况下我们会对其进行取log对数，然后对其化简最终得到（此处笔者不作化简过程，因为LaTeX公式太难写，笔者都是从word文档写好再转化过来的，读者可参照其他博主的相关文章）：
l o g L ( w ) = min ⁡ w − 1 n ∑ i = 1 n [ y i log ⁡ ( p ( x i ) ) + ( 1 − y i ) l o g ( 1 − p ( x i ) ) logL(w) = {minlimits_{w}{- frac{1}{n}{sumlimits_{i = 1}^{n}leftlbrack y_{i}{log(p({x_i}))} + (1 - {y_i})log(1 -p({x_i})) right.}}} logL(w)=wmin​−n1​i=1∑n​[yi​log(p(xi​))+(1−yi​)log(1−p(xi​))
​ 然后有了损失函数，我们就要利用梯度下降算法对其求解，因此，要对其求导求梯度。我们直接对w求导(求导过程请参照其他文章)
∂ l o g L ( w ) ∂ w = ∂ ∂ w ( − 1 n ∑ i = 1 n [ y i log ⁡ ( 1 1 + e − w T x i ) + ( 1 − y i ) log ⁡ ( 1 − 1 1 + e − w T x i ) ] ) = 1 n ∑ i = 1 n x i ( 1 1 + e − w T x i − y i ) begin{equation} begin{aligned} frac{partial logL(w)}{partial w} & =frac{partial}{partial w}({{{- frac{1}{n}{sumlimits_{i = 1}^{n}leftlbrack {y_{i}{logleft( frac{1}{1 + e^{- w^Tx_i}} right)} + left( {1 - y_i} right){logleft( {1 - frac{1}{1+e^{-w^Tx_i}}} right)}} right]}}}})\ &=frac{1}{n}sumlimits_{i=1}^nx_i(frac{1}{1+e^{-w^T{x_i}}}-y_i) end{aligned} end{equation} ∂w∂logL(w)​​=∂w∂​(−n1​i=1∑n​[yi​log(1+e−wTxi​1​)+(1−yi​)log(1−1+e−wTxi​1​)])=n1​i=1∑n​xi​(1+e−wTxi​1​−yi​)​​​
​ 因为在程序实现当中，为了程序的简介和提高运算的速度，我们往往是采用矩阵的形式进行预算。因此，要将上面的公式转化为矩阵相乘。当燃了，眼尖一眼就可以看出来。在转化之前，我们先介绍一下各个变量的矩阵表达
x = ( x 1 T ⋮ x n T ) = ( x 11 ⋯ x 1 p ⋮ ⋱ ⋮ x n 1 ⋯ x n p ) n ∗ p ; y = ( y 1 ⋮ y n ) n ∗ 1 ; p = ( p 1 ⋮ p n ) n ∗ 1 x = begin{pmatrix} {x_{1}}^T\ vdots \ {x_{n}}^T\ end{pmatrix} = begin{pmatrix} x_{11} & cdots & x_{1p} \ vdots & ddots & vdots \ x_{n1} & cdots & x_{np} \ end{pmatrix}_{n*p};y = begin{pmatrix} y_{1} \ vdots \ y_{n} \ end{pmatrix}_{n*1};p = begin{pmatrix} p_{1} \ vdots \ p_{n} \ end{pmatrix}_{n*1} x= ​x1​T⋮xn​T​ ​= ​x11​⋮xn1​​⋯⋱⋯​x1p​⋮xnp​​ ​n∗p​;y= ​y1​⋮yn​​ ​n∗1​;p= ​p1​⋮pn​​ ​n∗1​
​ n*p代表有n个数据，p代表每个数据有p个维度。里面的单个分量(如 x 1 , w 1 x_1,w_1 x1​,w1​)都是以列向量的形式存在.

​ 下面我们对梯度简单转化一下（暂时忽略1/n这个标量）
原式 = ( x 1 x 2 ⋯ x n ) ( 1 1 + e − w T x 1 − y 1 1 1 + e − w T x 2 − y 2 ⋮ 1 1 + e − w T x n − y n ) = x T ( p − y ) begin{equation} begin{aligned} 原式 &= left( {begin{matrix} x_{1} & x_{2} & cdots \ end{matrix}~~~~~~x_{n}~} right)begin{pmatrix} begin{matrix} {frac{1}{1 + e^{- w^{T}x_{1}}} - y_{1}} \ {frac{1}{1 + e^{- w^{T}x_{2}}} - y_{2}} \ vdots \ end{matrix} \ {frac{1}{1 + e^{- w^{T}x_{n}}} - y_{n}} \ end{pmatrix}\ &=x^T(p-y) end{aligned} end{equation} 原式​=(x1​​x2​​⋯​      xn​ ) ​1+e−wTx1​1​−y1​1+e−wTx2​1​−y2​⋮​1+e−wTxn​1​−yn​​ ​=xT(p−y)​​​
​ 因此
∂ l o g L ( w ) ∂ w = 1 n x T ( p − y ) frac{partial logL(w)}{partial w}= frac{1}{n}x^T(p-y) ∂w∂logL(w)​=n1​xT(p−y)

代码实现：

​ 注意了，这条直线实际上并不是一条，而是两条，一条是实际上区分数据的，另一条是训练出来的，可以看到基本重合，因此，可以说明训练效果挺好的。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
class Logistic_Regression():def __init__(self,x,y):''':param x: 样本值:param y: 标签值'''#在x的中增加一个全为1的维度。因为后续要将偏置项b写入w之中self.x=np.insert(x,x.shape[1],1,1)self.y=y#m，n为x的维度self.m,self.n&#shape#生成一个与（n,1）维的全为1的数组self.w&#s_like(x,shape=(self.n,1))def sigmoid(self,x):'''概率生成函数'''return 1/(np.exp(-x)+1)def train(self,epoch):''':param epoch: 训练轮次:return:'''for i in range(epoch):#将x*w后的值代入sigmoid函数生成概率p=self.sigmoid(self.x@self.w)#用交叉熵损失计算损失#1e-5为了防止下溢出loss=(-1/self.m*(self.y*np.log(p+1e-5)+(1-self.y)*np.log(1-p+1e-5))).sum()#梯度下降并更新梯度self.w-=0.T@(p-self.y)#设定停止训练的的loss最小if loss<0.01:breakdef predict(self,x):#给x增加一个全为1的维度x=np.insert(x,x.shape[1],1,1)#计算概率p=self.sigmoid(x@self.w)#大于0.5则认为是正类，反之则为负类result=np.array(p>0.5,dtype=int)return result#绘图函数
def plot_figure(x1,x2,norm_x2,y,w):''':param x1: x的第一个维度:param x2: x的第二个维度:param norm_x2: 加入噪声之后的x2:param y: 标签值:return:'''map_color={0:"r",1:"g"}y=[ map_color[i] for i in y.squeeze()]plt.scatter(x1,norm_x2,c=y)w=w.squeeze()x3=-(w[0]*x1.squeeze()+w[2])/w[1]plt.plot(x1.squeeze(),x3)plt.plot(x1,x2)plt.show()def main():x1=np.arange(0,100,1).reshape(-1,1) #生成x的第一个维度x2=np.arange(0,100,1).reshape(-1,1)#生成x的第二个维度norm_x2=x2&#vs(1,20,x2.shape)#给第二个维度加上噪声y=np.array(x2>norm_x2,dtype=int)#对大于norm的取1，反之取0x = np.concatenate([x1, norm_x2], axis=1)#联立x1,x2 (100,2)model=Logistic_Regression(x,y) #初始化模型并将数据传入ain(100000) #训练print(model.w) #打印模型参数，最后一个维度为偏置项bresult=model.predict(x) #预测plot_figure(x1,x2,norm_x2,result,model.w) #将预测结果画出if __name__ == '__main__':main()

softmax:

​ 对于多分类问题，原始的sigmoid函数已经不再适用，因此，要引入softmax函数。在引入之前，我们先对w矩阵做一下定义
w = ( w 1 T ⋮ w n T ) = ( w 11 w 12 ⋯ x 1 p w 21 w 22 ⋯ x 2 p ⋮ ⋮ ⋱ ⋮ w k 1 w k 2 ⋯ w k p ) k ∗ p w = begin{pmatrix} {w_{1}}^T\ vdots \ {w_{n}}^T \ end{pmatrix} = begin{pmatrix} w_{11} & w_{12}&cdots & x_{1p} \ w_{21} & w_{22} & cdots & x_{2p}\ vdots & vdots&ddots & vdots \ w_{k1} & w_{k2}&cdots & w_{kp} \ end{pmatrix}_{k*p} w= ​w1​T⋮wn​T​ ​= ​w11​w21​⋮wk1​​w12​w22​⋮wk2​​⋯⋯⋱⋯​x1p​x2p​⋮wkp​​ ​k∗p​
​ 该矩阵第i列代表的是第i类的参数。k代表是一列有多少个元素，实际上对应的就是x的维度，而p则代表是分类的类别数。

​ 而在此之前，我们也需要对标签值y的矩阵进行重定义，它将不再是一个列向量。而是一个矩阵。
y = ( y 1 T y 2 T ⋮ y n T ) = ( y 11 y 12 ⋯ y 1 p y 21 y 22 ⋯ y 2 p ⋮ ⋮ ⋱ ⋮ y n 1 y n 2 ⋯ y n p ) n ∗ p y=begin{pmatrix} y_1^T \ y_2^T \ vdots \ y_n ^T end{pmatrix} =begin{pmatrix} y_{11} & y_{12} & cdots & y_{1p}\ y_{21} & y_{22} & cdots & y_{2p}\ vdots & vdots & ddots & vdots\ y_{n1} & y_{n2} & cdots & y_{np}\ end{pmatrix}_{n*p} y= ​y1T​y2T​⋮ynT​​ ​= ​y11​y21​⋮yn1​​y12​y22​⋮yn2​​⋯⋯⋱⋯​y1p​y2p​⋮ynp​​ ​n∗p​
​ 其中n代表样本数量，p代表类别数量。

​ 因此，我们定义softmax函数：
g ( x ) = 1 ∑ i = 1 k e w T x ( e p 1 e p 2 ⋮ e p n ) g(x)=frac{1}{sumlimits_{i=1}^ke^{w^Tx}} begin{pmatrix} e_{p_1} \ e_{p_2} \ vdots \ \ e_{p_n} end{pmatrix} g(x)=i=1∑k​ewTx1​ ​ep1​​ep2​​⋮epn​​​ ​
​ 其中， e p i e_{p_i} epi​​代表第i类的概率。所以某一类的概率就可以写成如下：
p ( y = y j ∣ x i ; w j ) = e w T x i ∑ i = 1 k e w T x i p(y=y_j|x_i;w_j)=frac{e^{w^Tx_i}}{sumlimits_{i=1}^ke^{w^Tx_i}} p(y=yj​∣xi​;wj​)=i=1∑k​ewTxi​ewTxi​​
​ 以上代表第i个数据，属于第j类的概率。

损失函数：

​ softmax的损失函数如下：
l o g L ( w ) = − 1 n [ ∑ i = 1 n ∑ j = 1 p y i j l o g ( e w j T x i ∑ l = 1 k w l T x i ) ] logL(w)=-frac{1}{n}[sumlimits_{i=1}^nsumlimits_{j=1}^p{y_{ij} {log(frac{e^{w_j^Tx_i}}{sumlimits_{l=1}^kw_l^Tx_i})}}] logL(w)=−n1​[i=1∑n​j=1∑p​yij​log(l=1∑k​wlT​xi​ewjT​xi​​)]
​ 其中 1 { y i = j } 1left{y_{i}=jright} 1{yi​=j}是示性函数，为真是取1，反之取0。

​ 有了损失函数，我们就可以依照上面的方法对其进行求导然后利用梯度下降求解啦！
∂ l o g L ( w ) ∂ w = ( ∂ l o g L ( w ) ∂ w j 1 ∂ l o g L ( w ) ∂ w j 2 ⋯ ∂ l o g L ( w ) ∂ w j p ) begin{equation} begin{aligned} frac{partial logL(w)}{partial w}= begin{pmatrix} frac{partial logL(w)}{partial w_{j1}} & frac{partial logL(w)}{partial w_{j2}} & cdots frac{partial logL(w)}{partial w_{jp}} end{pmatrix} end{aligned} end{equation} ∂w∂logL(w)​=(∂wj1​∂logL(w)​​∂wj2​∂logL(w)​​⋯∂wjp​∂logL(w)​​)​​​
​ 这里的 w j p w_{jp} wjp​代表的是对第p个分类的w求偏导，下面我们单独提取出一个进行求导（ w j w_{j} wj​表示某一类）
∂ l o g L ( w ) ∂ w j = ∂ ∂ w j ( − 1 n [ ∑ i = 1 n ∑ j = 1 p y i j l o g ( e w j T x i ∑ l = 1 k w l T x i ) ] ) = − 1 n [ ∑ i = 1 n x i ( y i j − p ( y i = j ∣ x i ; w j ) ) ] begin{equation} begin{aligned} frac{partial logL(w)}{partial w_{j}}&= frac{partial }{partial w_{j}}left( -frac{1}{n}[sumlimits_{i=1}^nsumlimits_{j=1}^p{y_{ij} {log(frac{e^{w_j^Tx_i}}{sumlimits_{l=1}^kw_l^Tx_i})}}]right )\ &=-frac{1}{n}[sumlimits_{i=1}^n{x_i}left(y_{ij} -p(y_i=j|x_i;w_j)right)] end{aligned} end{equation} ∂wj​∂logL(w)​​=∂wj​∂​ ​−n1​[i=1∑n​j=1∑p​yij​log(l=1∑k​wlT​xi​ewjT​xi​​)] ​=−n1​[i=1∑n​xi​(yij​−p(yi​=j∣xi​;wj​))]​​​
​ 那么接下来就把它化为矩阵表达就可以完美收工啦(以下转化暂时忽略掉-1/n这个标量)
∂ l o g L ( w ) ∂ w j = ( x 1 x 2 ⋯ x n ) ( y 1 j − p 1 j y 2 j − p 2 j ⋮ y n j − p n j ) = x T ∗ ( y j − p j ) begin{equation} begin{aligned} frac{partial logL(w)}{partial w_{j}}&= begin{pmatrix} x_1 & x_2 & cdots & x_n end{pmatrix} begin{pmatrix} y_{1j}-p_{1j} \ y_{2j}-p_{2j} \ vdots \ y_{nj}-p_{nj} \ end{pmatrix} \ &=x^T*(y_j-p_j) end{aligned} end{equation} ∂wj​∂logL(w)​​=(x1​​x2​​⋯​xn​​) ​y1j​−p1j​y2j​−p2j​⋮ynj​−pnj​​ ​=xT∗(yj​−pj​)​​​
​ 因此：
∂ l o g L ( w ) ∂ w = ( ∂ l o g L ( w ) ∂ w j 1 ∂ l o g L ( w ) ∂ w j 2 ⋯ ∂ l o g L ( w ) ∂ w j p ) = ( x T ∗ ( y j 1 − p j 1 ) x T ∗ ( y j 2 − p j 2 ) ⋯ x T ∗ ( y j p − p j p ) ) = x T ( y j 1 − p j 1 y j 2 − p j 2 ⋯ y j p − p j p ) = − 1 n x T ( y − p ) begin{equation} begin{aligned} frac{partial logL(w)}{partial w}&= begin{pmatrix} frac{partial logL(w)}{partial w_{j1}} & frac{partial logL(w)}{partial w_{j2}} & cdots frac{partial logL(w)}{partial w_{jp}} end{pmatrix}\ &=begin{pmatrix} x^T*(y_{j1}-p_{j1}) & x^T*(y_{j2}-p_{j2}) & cdots & x^T*(y_{jp}-p_{jp}) & end{pmatrix}\ &=x^T begin{pmatrix} y_{j1}-p_{j1} & y_{j2}-p_{j2} & cdots & y_{jp}-p_{jp} end{pmatrix}\ &=-frac{1}{n}x^Tleft( y-pright) end{aligned} end{equation} ∂w∂logL(w)​​=(∂wj1​∂logL(w)​​∂wj2​∂logL(w)​​⋯∂wjp​∂logL(w)​​)=(xT∗(yj1​−pj1​)​xT∗(yj2​−pj2​)​⋯​xT∗(yjp​−pjp​)​​)=xT(yj1​−pj1​​yj2​−pj2​​⋯​yjp​−pjp​​)=−n1​xT(y−p)​​​

​ 此处直接将上文忽略的-1/n写回。

代码实现：

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
class Logistic_Regression():def __init__(self,x,y):''':param x: 样本值:param y: 标签值'''#在x的中增加一个全为1的维度。因为后续要将偏置项b写入w之中self.x=np.insert(x,x.shape[1],1,1)self.y=np.array(y.squeeze(),dtype=int).reshape(-1)#计算y的类别数self.class_num=len(set(self.y))#m，n为x的维度self.m,self.n&#shape#构造one_hot 编码的标签值hot_y&#s_like(y,shape=(y.shape[0],self.class_num))hot_y[np.shape[0]),self.y]=1#生成一个(维度，类别) 也就是 (n,self.class_num)self.w&#s_like(x,shape=(self.n,self.class_num))def softmax(self,x):'''概率生成函数x:(400,4)'''exp&#p(x)ex_sum=np.sum(exp,axis=1).reshape(-1,1)result=exp/ex_sumreturn resultdef train(self,epoch):''':param epoch: 训练轮次:return:'''for i in range(epoch):#将x*w后的值代入sigmoid函数生成概率p=self.softmax(self.x@self.w)#(400,4)#用交叉熵损失计算损失loss=-(1/self.m)*np.sum((hot_y))print("函数损失:",loss)#梯度下降并更新梯度GD=(-1/self.m)*T@(hot_y-p)self.w-=0.001*GDdef predict(self,x):#给x增加一个全为1的维度x=np.insert(x,x.shape[1],1,1)#计算概率p=self.softmax(x@self.w)#按行求最大值并返回其对应的indexresult=np.argmax(p,axis=1)return result#绘图函数
def plot_figure(x,y):map_color={0:"r",1:"g",2:"b",3:"y"}color=[map_color[i] for i in y]plt.scatter(x[:,0],x[:,1],c=color)plt.show()
def main():#生成第一个团和对应的标签值class1_x1&#atenate([vs(2,2,(100,1)),vs(1,2,(100,1))],axis=1)y1&#s_like(class1_x1,shape=(class1_x1.shape[0],1))# 生成第二个团和对应的标签值class1_x2&#atenate([vs(20,2,(100,1)),vs(1,2,(100,1))],axis=1)y2 = np.ones_like(class1_x2, shape=(class1_x2.shape[0], 1))# 生成第三个团和对应的标签值class1_x3&#atenate([vs(-5,2,(100,1)),vs(10,2,(100,1))],axis=1)y3 = np.ones_like(class1_x3, shape=(class1_x3.shape[0], 1))*2# 生成第四个团和对应的标签值class1_x4 = np.concatenate([vs(20, 2, (100, 1)), vs(10, 2, (100, 1))], axis=1)y4 = np.ones_like(class1_x4, shape=(class1_x4.shape[0], 1)) * 3#将所有的x合在一起作为训练数据，将所有的y合在一起作为训练数据y&#atenate([y1,y2,y3,y4],axis=0)x&#atenate([class1_x1,class1_x2,class1_x3,class1_x4],axis=0)model=Logistic_Regression(x,y) #初始化模型并将数据传入ain(100000) #训练print(model.w) #打印模型参数，最后一个维度为偏置项bresult=model.predict(x) #预测plot_figure(x,result) #将预测结果画出
if __name__ == '__main__':main()

结束语：

​ 好啦，以上就是这篇文章的全部内容了。笔者也是在学习中，对于公式如何进行矩阵化的表达是很多大佬都忽略掉的问题。当燃也有可能是我太笨了，才会产生这种问题，亦或者是对线性代数缺乏感性的认识。也许对于很多人而言，可窥一眼而知全貌。但不管怎么说，我都希望能够帮助到与我有着同样的小伙伴。另外，本文公式中有一些不严谨的地方，有任何错误还望能够不吝赐教。就这样把，阿里嘎多（用LateX敲公式真累啊，什么时候CSDN能直接上传word文档）。