逻辑回归算法实现

---

---

算法原理

模型

逻辑回归模型是一个二分类的对数线性模型。令样本数据集为
{ x i ; y i } i = 1 n , x i ∈ R t , y i ∈ { 0 , 1 } \{x_i;y_i\}^n_{i=1},x_i\in\mathbb{R}^t,y_i\in\{0,1\} {xi​;yi​}i=1n​,xi​∈Rt,yi​∈{0,1}
其中，$x_i$是样本的特征，$y_i$是样本的类别，限定为0或1。逻辑回归模型的数学表达式为
$$f(x_i)=\frac{1}{1+e^{-({\omega }^T{x}_i+b)}},i=1,2,...,n$$
单单看这个表达式可以发现，它跟上节讲到的感知机模型是很类似的，感知机使用sign函数对线性表达式进行转化，而逻辑回归是用了sigmoid函数，这也是为什么感知机中的类别标签是1和-1，而逻辑回归的类别标签是0和1。但是在算法思想上两者却完全不同，感知机利用了错误驱动，通过错误分类集来构造损失函数，然后使用梯度下降法进行求解。而接下来我们会讲到，逻辑回归是通过构造概率模型并最大化概率的方式来进行分类。

原理

首先，模型的前提是所有样本数据都是独立同分布的。

单个样本$x_i$

现在已知有一个样本数据($x_i$，$y_i$)，我们的目的是希望将$x_i$代入模型
$$f(x_i)=\frac{1}{1+e^{-({\omega }^T{x}_i+b)}}$$
得到的$f(x_i)$刚好等于$y_i$。所以在这里要先定义一个条件概率，表示在已知样本是$x_i$的条件下，类别$y_i$是1的概率
$$p_1=P(y_i=1|x_i)=\frac{1}{1+e^{-({\omega }^T{x}_i+b)}}$$
当类别$y_i$等于1时，$p_1$要尽可能的大。由$p_1$的公式可以得到在已知样本是$x_i$的条件下，类别$y_i$是0的概率为
$$p_0=P(y_i=0|x_i)=1-P(y_i=1|x_i)=\frac{e^{-({\omega }^T{x}_i+b)}}{1+e^{-({\omega }^T{x}_i+b)}}$$
当类别$y_i$等于0时，$p_0$要尽可能的大。

可以将$p_1$和$p_0$整合起来，令
$$p=P(y_i|x_i)=y_i{p}_1+(1-y_i)p_0=\{\begin{array}{ll}{p}_1& y_i=1\\ p_0& y_i=0\end{array}$$
这样，无论$y_i$是0还是1，都是要让$p$尽可能的大，即最大化。

整个数据集

对于单个样本数据($x_i$，$y_i$)，需要让$p$尽可能的大，那推广到整个数据集上，就是要求解
$$\underset{\omega ,b}{\max }\prod _{i=1}^{n}P(y_i|x_i)$$
这里之所以将每个样本数据对应的$p$相乘，是因为每个样本数据都是独立同分布的。另外，由于每个$p$的数值都在0和1之间，因此它们的积也是在0和1之间。为了将乘法拆开，需要做一个负对数的转换（一般采用自然对数）
$$\underset{\omega ,b}{\max }\prod _{i=1}^{n}P(y_i|x_i)\to \underset{\omega ,b}{\max }[log(\prod _{i=1}^{n}P(y_i|x_i))]\to \underset{\omega ,b}{\min }[-log(\prod _{i=1}^{n}P(y_i|x_i))]\to \underset{\omega ,b}{\min }[-\sum _{i=1}^n(logP(y_i|x_i))]$$
为了避免误会，箭头代表问题的等价转换。第一个箭头转换是因为对数函数是增函数，第二个箭头是因为添加负号，同时变成最小化问题，第三个箭头是根据对数的性质。进一步，我们就得到了逻辑回归的损失函数，也就是第三个箭头右边的公式的均值化，详细展开如下。
$$L(\omega ,b)=-\frac{1}{n}\sum _{i=1}^n(logP(y_i|x_i))=-\frac{1}{n}\sum _{i=1}^n[y_{i}log(p_1)+(1-y_i)log(p_0)]=-\frac{1}{n}\sum _{i=1}^n[y_{i}log(\frac{1}{1+e^{-({\omega }^T{x}_i+b)}})+(1-y_i)log(\frac{e^{-({\omega }^T{x}_i+b)}}{1+e^{-({\omega }^T{x}_i+b)}})]$$

模型求解

采用梯度下降法求解，更新$\omega$和$b$的值，不断迭代，从而最小化损失函数。
$$\omega \to \omega -\lambda \frac{\partial L}{\partial \omega }$$
$$b\to b-\lambda \frac{\partial L}{\partial b}$$
其中，$\lambda$是学习步长
$$\frac{\partial L}{\partial \omega }=-\frac{1}{n}\sum _{i=1}^n{x}_i(y_i-f(x_i))$$
$$\frac{\partial L}{\partial b}=-\frac{1}{n}\sum _{i=1}^n(y_i-f(x_i))$$

程序实现

sigmoid函数

def sigmoid(X):
    y = 1 / (1+np.exp(-X))
    return y

初始化函数

def __init__(self, alpha=0.01, iteration=1000):
    """
    alpha     学习步长
    w         权重向量
    iteration 最大迭代次数
    error     记录损失函数值
    """
    self.alpha = alpha
    self.iteration = iteration

训练函数

def fit(self, X, y):
    [data_num, fea_num] = X.shape
    X_ = np.vstack((np.ones(data_num).reshape(1,-1), X.T))
    y_ = y.reshape(-1,1)
    self.w = np.zeros(fea_num+1).reshape(-1,1)
    self.error = []
    i = 1

    while i<self.iteration:
        # 迭代w值
        dif_y = y_ - sigmoid(np.dot(self.w.T,X_)).reshape(-1,1)
        self.w = self.w + self.alpha * 1/data_num * np.dot(X_, dif_y)

        # 记录损失
        p1 = sigmoid(np.dot(self.w.T,X_).T)
        temp = y_*np.log(p1) + (1-y_)*np.log(1-p1)
        self.error.append(-temp.mean())

        # 迭代次数
        if (i+1)%100==0:
            print('epoch--', i+1)
        i += 1

训练函数的核心代码在第十一和第十二行，是参数$\omega$的迭代更新公式，参考模型求解部分的公式，代码中将权重参数$\omega$和偏置参数$b$合并在一起进行计算。可转化成矩阵计算的形式如下图（*代表矩阵乘法）。

图中显示的各矩阵的维数分别为：
$$\omega :(t+1)\ast 1$$
$$X_{:}(t+1)\ast n$$
$$y和f(X):n\ast 1$$
其中，$\omega$包含了偏置$b$，所以需要再加上一维变成$t+1$维。另外，逻辑回归不一定能将样本全部正确分类，所以需要设置迭代次数。

预测函数

def predict(self, X):
    X_ = np.vstack((np.ones(X.shape[0]).reshape(1,-1), X.T))
    sig = sigmoid(np.dot(self.w.T, X_))[0]
    return (sig>=0.5)*1

整合全部代码

class LogisticRegression:

    def __init__(self, alpha=0.01, iteration=1000):
        """
        alpha     学习步长
        w         权重向量
        iteration 最大迭代次数
        error     记录损失函数值
        """
        self.alpha = alpha
        self.iteration = iteration
    
    def fit(self, X, y):
        [data_num, fea_num] = X.shape
        X_ = np.vstack((np.ones(data_num).reshape(1,-1), X.T))
        y_ = y.reshape(-1,1)
        self.w = np.zeros(fea_num+1).reshape(-1,1)
        self.error = []
        i = 1
        
        while i<self.iteration:
            # 迭代w值
            dif_y = y_ - sigmoid(np.dot(self.w.T,X_)).reshape(-1,1)
            self.w = self.w + self.alpha * 1/data_num * np.dot(X_, dif_y)
                
            # 记录损失
            p1 = sigmoid(np.dot(self.w.T,X_).T)
            temp = y_*np.log(p1) + (1-y_)*np.log(1-p1)
            self.error.append(-temp.mean())
            
            # 迭代次数
            if (i+1)%100==0:
                print('epoch--', i+1)
            i += 1
    
    def predict(self, X):
        X_ = np.vstack((np.ones(X.shape[0]).reshape(1,-1), X.T))
        sig = sigmoid(np.dot(self.w.T, X_))[0]
        return (sig>=0.5)*1

实例化演示

导入相关库和数据并查看数据分布

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

iris = load_iris()
X = iris.data[:100,:2]
y = (iris.target[:100]==1)*1

plt.scatter(X[y==0,0], X[y==0,1], color='red')
plt.scatter(X[y==1,0], X[y==1,1], color='green')

实例化函数并画出划分平面

# 训练
clf=LogisticRegression()
clf.fit(X,y)
line1 = X[:,0]
line2 = -(clf.w[0]+clf.w[1]*line1)/clf.w[2]

# 画图
plt.scatter(X[y==0,0],X[y==0,1],color='red')
plt.scatter(X[y==1,0],X[y==1,1],color='green')
plt.plot(line1,line2)

至此，逻辑回归算法讲解完毕，以上是个人对逻辑回归算法的一些理解，如有错误欢迎指出。