逻辑回归 二/多分类 +Python代码实现_python实现多分类问题模型 代码实现

前言：

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

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

sigmoid:

$$p=\frac{1}{1+e^{-w^{T}x}}$$

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

损失函数及其梯度计算：

​ 显而易见，我们可以直接将极大似然估计作为损失函数，当燃了，损失函数一般我们是越小越好。而极大似然估计是越大越好，so？我们要对极大似然取反，最终损失函数应当如下：
$$L(w)=\underset{w}{\min }-\frac{1}{m}\prod _{i=1}^n{p}^{y_i}(1-p{)}^{1-y_i}$$
​ 为了便于运算，一般情况下我们会对其进行取log对数，然后对其化简最终得到（此处笔者不作化简过程，因为LaTeX公式太难写，笔者都是从word文档写好再转化过来的，读者可参照其他博主的相关文章）：
$$logL(w)=\underset{w}{\min }-\frac{1}{n}\sum _{i=1}^n[y_i\log (p(x_i))+(1-y_i)log(1-p(x_i))$$
​ 然后有了损失函数，我们就要利用梯度下降算法对其求解，因此，要对其求导求梯度。我们直接对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}{\sum\limits_{i = 1}^{n}\left\lbrack {y_{i}{\log\left( \frac{1}{1 + e^{- w^Tx_i}} \right)} + \left( {1 - y_i} \right){\log\left( {1 - \frac{1}{1+e^{-w^Tx_i}}} \right)}} \right]}}}})\\ &=\frac{1}{n}\sum\limits_{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=\left(\begin{array}{c}{x_1}^T\\ ⋮\\ {x_n}^T\end{array}\right)={\left(\begin{array}{ccc}{x}_{11}& \cdots & x_{1p}\\ ⋮& \ddots & ⋮\\ x_{n1}& \cdots & x_{np}\end{array}\right)}_{n\ast p};y={\left(\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right)}_{n\ast 1};p={\left(\begin{array}{c}{p}_1\\ ⋮\\ p_n\end{array}\right)}_{n\ast 1}$$
​ n*p代表有n个数据，p代表每个数据有p个维度。里面的单个分量(如$x_1,w_1$)都是以列向量的形式存在.

​ 下面我们对梯度简单转化一下（暂时忽略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)​​​
​ 因此
$$\frac{\partial logL(w)}{\partial w}=\frac{1}{n}{x}^T(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=self.x.shape
        #生成一个与（n,1）维的全为1的数组
        self.w=np.zeros_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.0001*self.x.T@(p-self.y)
            #设定停止训练的的loss最小
            if loss<0.01:
                break
    def 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+stats.norm.rvs(1,20,x2.shape)#给第二个维度加上噪声
    y=np.array(x2>norm_x2,dtype=int)#对大于norm的取1，反之取0
    x = np.concatenate([x1, norm_x2], axis=1)#联立x1,x2 (100,2)
    model=Logistic_Regression(x,y) #初始化模型并将数据传入
    model.train(100000) #训练
    print(model.w) #打印模型参数，最后一个维度为偏置项b
    result=model.predict(x) #预测
    plot_figure(x1,x2,norm_x2,result,model.w) #将预测结果画出

if __name__ == '__main__':
    main()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

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}{\sum\limits_{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}$代表第i类的概率。所以某一类的概率就可以写成如下：
$$p(y=y_j|x_i;w_j)=\frac{e^{w^T{x}_i}}{\sum _{i=1}^k{e}^{w^T{x}_i}}$$
​ 以上代表第i个数据，属于第j类的概率。

损失函数：

​ softmax的损失函数如下：
$$logL(w)=-\frac{1}{n}[\sum _{i=1}^n\sum _{j=1}^p{y}_{ij}log(\frac{e^{w_j^T{x}_i}}{\sum _{l=1}^k{w}_l^T{x}_i})]$$
​ 其中 1 { y i = j } 1\left\{y_{i}=j\right\} 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_{jp}$代表的是对第p个分类的w求偏导，下面我们单独提取出一个进行求导（$w_j$表示某一类）
$$\begin{array}{c}\begin{array}{rl}\frac{\partial logL(w)}{\partial w_j}& =\frac{\partial }{\partial w_j}\left(-\frac{1}{n}[\sum _{i=1}^n\sum _{j=1}^p{y}_{ij}log(\frac{e^{w_j^T{x}_i}}{\sum _{l=1}^k{w}_l^T{x}_i})]\right)\\ & =-\frac{1}{n}[\sum _{i=1}^n{x}_i\left(y_{ij}-p(y_i=j|x_i;w_j)\right)]\end{array}\end{array}$$
​ 那么接下来就把它化为矩阵表达就可以完美收工啦(以下转化暂时忽略掉-1/n这个标量)
$$\begin{array}{c}\begin{array}{rl}\frac{\partial logL(w)}{\partial w_j}& =\left(\begin{array}{cccc}{x}_1& x_2& \cdots & x_n\end{array}\right)\left(\begin{array}{c}{y}_{1j}-p_{1j}\\ y_{2j}-p_{2j}\\ ⋮\\ y_{nj}-p_{nj}\end{array}\right)\\ & =x^T\ast (y_j-p_j)\end{array}\end{array}$$
​ 因此：
$$\begin{array}{c}\begin{array}{rl}\frac{\partial logL(w)}{\partial w}& =\left(\begin{array}{ccc}\frac{\partial logL(w)}{\partial w_{j1}}& \frac{\partial logL(w)}{\partial w_{j2}}& \cdots \frac{\partial logL(w)}{\partial w_{jp}}\end{array}\right)\\ & =\left(\begin{array}{ccccc}{x}^T\ast (y_{j1}-p_{j1})& x^T\ast (y_{j2}-p_{j2})& \cdots & x^T\ast (y_{jp}-p_{jp})& \end{array}\right)\\ & =x^T\left(\begin{array}{cccc}{y}_{j1}-p_{j1}& y_{j2}-p_{j2}& \cdots & y_{jp}-p_{jp}\end{array}\right)\\ & =-\frac{1}{n}{x}^T\left(y-p\right)\end{array}\end{array}$$

​ 此处直接将上文忽略的-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=self.x.shape
        #构造one_hot 编码的标签值
        self.onehot_y=np.zeros_like(y,shape=(y.shape[0],self.class_num))
        self.onehot_y[np.arange(self.y.shape[0]),self.y]=1
        #生成一个(维度，类别) 也就是 (n,self.class_num)
        self.w=np.zeros_like(x,shape=(self.n,self.class_num))
    def softmax(self,x):
        '''
        概率生成函数
        x:(400,4)
        '''
        exp=np.exp(x)
        ex_sum=np.sum(exp,axis=1).reshape(-1,1)
        result=exp/ex_sum
        return result
    def 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((p*self.onehot_y))
            print("函数损失:",loss)
            #梯度下降并更新梯度
            GD=(-1/self.m)*self.x.T@(self.onehot_y-p)
            self.w-=0.001*GD
    def predict(self,x):
        #给x增加一个全为1的维度
        x=np.insert(x,x.shape[1],1,1)
        #计算概率
        p=self.softmax(x@self.w)
        #按行求最大值并返回其对应的index
        result=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=np.concatenate([stats.norm.rvs(2,2,(100,1)),stats.norm.rvs(1,2,(100,1))],axis=1)
    y1=np.zeros_like(class1_x1,shape=(class1_x1.shape[0],1))
    # 生成第二个团和对应的标签值
    class1_x2=np.concatenate([stats.norm.rvs(20,2,(100,1)),stats.norm.rvs(1,2,(100,1))],axis=1)
    y2 = np.ones_like(class1_x2, shape=(class1_x2.shape[0], 1))
    # 生成第三个团和对应的标签值
    class1_x3=np.concatenate([stats.norm.rvs(-5,2,(100,1)),stats.norm.rvs(10,2,(100,1))],axis=1)
    y3 = np.ones_like(class1_x3, shape=(class1_x3.shape[0], 1))*2
    # 生成第四个团和对应的标签值
    class1_x4 = np.concatenate([stats.norm.rvs(20, 2, (100, 1)), stats.norm.rvs(10, 2, (100, 1))], axis=1)
    y4 = np.ones_like(class1_x4, shape=(class1_x4.shape[0], 1)) * 3
    #将所有的x合在一起作为训练数据，将所有的y合在一起作为训练数据
    y=np.concatenate([y1,y2,y3,y4],axis=0)
    x=np.concatenate([class1_x1,class1_x2,class1_x3,class1_x4],axis=0)
    model=Logistic_Regression(x,y) #初始化模型并将数据传入
    model.train(100000) #训练
    print(model.w) #打印模型参数，最后一个维度为偏置项b
    result=model.predict(x) #预测
    plot_figure(x,result) #将预测结果画出
if __name__ == '__main__':
    main()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

结束语：

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