二项逻辑斯蒂回归（逻辑回归）

二项逻辑斯蒂回归

1 定义

二项逻辑斯蒂回归模型是如下的条件概率分布：
$$\begin{gathered} P(Y=1|x)=\frac{\exp (\omega \cdot x+b)}{1+\exp (\omega \cdot x+b)} \\ P(Y=0|x)=\frac{1}{1+\exp (\omega \cdot x+b)} \end{gathered}$$
这里， $x\in R^n$ 是输入， Υ ∈ { 0 , 1 } \Upsilon\in\left\{0,1\right\} Υ∈{0,1} 是输出，$\omega \in R^n$ 和 $b\in R$ 是参数，$\omega$ 称为权值向量，$b$ 称为偏置，$\omega \cdot x$ 为 $\omega$ 和 $x$ 的内积。

可以将权值向量和输入向量加以扩充，仍记作 $\omega ,x$ ，即
$$\begin{array}{l}\omega =\left(\begin{array}{cccccc}{\omega }^{(1)}& {\omega }^{(2)}& {\omega }^{(3)}& \cdots & {\omega }^{(n)}& b\end{array}\right)\\ x=\left(\begin{array}{cccccc}{x}^{(1)}& x^{(2)}& x^{(3)}& \cdots & x^{(n)}& 1\end{array}\right)\end{array}$$
这时，逻辑斯蒂回归模型如下：
$$\begin{gathered} P(Y=1|x)=\frac{\exp (\omega \cdot x)}{1+\exp (\omega \cdot x)} \\ P(Y=0|x)=\frac{1}{1+\exp (\omega \cdot x)} \end{gathered}$$

2 模型参数估计

应用极大似然估计法估计模型参数，从而得到逻辑斯蒂回归模型。

对数似然函数为
$$L(\omega )=\sum _{i=1}^N\left[y_i(\omega \cdot x_i)-\log \left(1+\exp (\omega \cdot x_i)\right)\right]$$
从理论上来讲，直接求 $L(\omega )$ 的极大值，就能得到 $\omega$ 的估计值。但在实际的数据场景中，我们常采用梯度下降法进行求解。

$L(\omega )$ 对 $\omega$ 的每一个元素求偏导，
$$\frac{\partial L(\omega )}{\partial {\omega }^{(j)}}=\sum _{i=1}^N\left[x_i^{(j)}{y}_i-\frac{\exp (\omega \cdot x_i)x_i^{(j)}}{1+\exp (\omega \cdot x_i)}\right]$$
再由所有偏导组成向量，得到的就是梯度。

3 学习算法¹

输入：训练数据集 T = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ⋯   , ( x N , y N ) } T=\left\{\left(x_{1}, y_{1}\right),\left(x_{2}, y_{2}\right), \cdots,\left(x_{N}, y_{N}\right)\right\} T={(x1​,y1​),(x2​,y2​),⋯,(xN​,yN​)}，其中 $x_i\in \chi =R^n$， y i ∈ Υ = { − 1 , + 1 } y_{i}\in\Upsilon=\left\{-1,+1\right\} yi​∈Υ={−1,+1}，$i=1,2,\cdots ,N$；学习率 $\eta \left(0<\eta \le 1\right)$，又称为步长；梯度阈值 $\delta$ ；

输出：$\omega$；

1. 选取初始值 ${\omega }_0$；

2. 在训练集中选取数据 $\left(x_i,y_i\right)$；

3. 如果 $\exists j\in [1,n]$ ，使得：
   $$\left[x_i^{(j)}{y}_i-\frac{\exp (\omega \cdot x_i)x_i^{(j)}}{1+\exp (\omega \cdot x_i)}\right]>\delta$$
   则：
   $$\omega \leftarrow \omega +\eta \left[x_i^{(j)}{y}_i-\frac{\exp (\omega \cdot x_i)x_i^{(j)}}{1+\exp (\omega \cdot x_i)}\right]$$

4. 转至2，直到对于 $\forall j\in [1,n]$ ，都有
   $$\left[x_i^{(j)}{y}_i-\frac{\exp (\omega \cdot x_i)x_i^{(j)}}{1+\exp (\omega \cdot x_i)}\right]\le \delta$$

4 代码实现

1. 了解数据

   • 输入向量 $\chi$ 和 ${\rm Y}$
     $$\begin{array}{l}\chi =\left(\begin{array}{ccccc}{x}_1^{(1)}& x_1^{(2)}& x_1^{(3)}& \cdots & x_1^{(n)}\\ x_2^{(1)}& x_2^{(2)}& x_2^{(3)}& \cdots & x_2^{(n)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_N^{(1)}& x_N^{(2)}& x_N^{(3)}& \cdots & x_N^{(n)}\end{array}\right){\rm Y}=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_N\end{array}\right)\end{array}$$
   • 扩充后的输入向量
     $$\chi =\left(\begin{array}{cccccc}{x}_1^{(1)}& x_1^{(2)}& x_1^{(3)}& \cdots & x_1^{(n)}& 1\\ x_2^{(1)}& x_2^{(2)}& x_2^{(3)}& \cdots & x_2^{(n)}& 1\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ x_N^{(1)}& x_N^{(2)}& x_N^{(3)}& \cdots & x_N^{(n)}& 1\end{array}\right)$$
     对应的代码为：

   def preprocessing(X):
       X_plus = np.ones(X.shape[0]).reshape(-1, 1)
       X_new = np.hstack([X, X_plus])
       return X_new
   X = preprocessing(X)
   

2. 初始化参数。

   • 权值向量 $\omega$ 是用来和样本 $x_i$ 求内积的一个向量，对应于扩充后的输入向量：

     $$\omega =\left(\begin{array}{cccccc}{\omega }^{(1)}& {\omega }^{(2)}& {\omega }^{(3)}& \cdots & {\omega }^{(n)}& b\end{array}\right)=\left(\begin{array}{ccccc}0& 0& 0& \cdots & 0\end{array}\right)$$
     对应的代码为：

     w = np.zeros(X.shape[1]).reshape(1, -1)
     

   • 超参数 $\eta$ ，赋予默认值 1
     $$\eta =1$$

     eta = 1
     

   • 梯度更新的阈值 $\delta$ ，赋予默认值 0.01

     delta = 1e-2
     

3. 计算梯度

   • 定义
     $$\nabla L(\omega )=\left(\frac{\partial L(\omega )}{\partial {\omega }^{(1)}}\cdots \frac{\partial L(\omega )}{\partial {\omega }^{(n+1)}}\right)$$

   • 公式解析
     $$\frac{\partial L(\omega )}{\partial {\omega }^{(j)}}=\sum _{i=1}^N\left[x_i^{(j)}{y}_i-\frac{\exp (\omega \cdot x_i)x_i^{(j)}}{1+\exp (\omega \cdot x_i)}\right]=\sum _{i=1}^N\left[x_i^{(j)}\left(y_i-\frac{\exp (\omega \cdot x_i)}{1+\exp (\omega \cdot x_i)}\right)\right]$$
     等式左边，是对 $\omega$ 的第 $j$ 个元素求偏导，得到的就是梯度第 $j$ 个元素的值。

     等式右边是一个求和公式。先看求和符号，$i\in [1,N]$ 表示所有的样本；再来看求和的具体内容，$x_i$ 是第 $i$ 个样本的特征向量，$x_i^{(j)}$ 第 $i$ 个样本的第 $j$ 个特征，即特征向量的第 $j$ 个元素，二者有如下关系：
     $$x_i=\left(\begin{array}{ccccccc}{x}_i^{(1)}& x_i^{(2)}& \cdots & x_i^{(j)}& \cdots & x_i^{(n)}& 1\end{array}\right)$$
     所以等式右边表示的是：所有样本的、第 $j$ 个特征的、计算值的求和。

   • 公式计算过程

     由内向外看，对于样本 $x_i$ ，最内部的内积，在这里是一个数值：
     $$z_i=\omega \cdot x_i=\left(\begin{array}{cccccc}{\omega }^{(1)}& {\omega }^{(2)}& {\omega }^{(3)}& \cdots & {\omega }^{(n)}& b\end{array}\right)\left(\begin{array}{c}{x}_i^{(1)}\\ x_i^{(2)}\\ \cdots \\ x_i^{(j)}\\ \cdots \\ x_i^{(n)}\\ 1\end{array}\right)={\omega }^{(1)}{x}_i^{(1)}+{\omega }^{(2)}{x}_i^{(2)}+\cdots +{\omega }^{(n)}{x}_i^{(n)}+b$$
     紧接着点乘的，是一个Sigmoid函数，变换后，仍旧是一个数值：
     $${\hat{z}}_i=\frac{\exp (z_i)}{1+\exp (z_i)}$$

     紧接着，是与 $y_i$ 做差，再与 $x_i^{(j)}$ 相差，这两个变量都是数值，所以第 $i$ 个样本、第 $j$ 个特征最终的计算值为：
     $$x_i^{(j)}\left(y_i-{\hat{z}}_i\right)$$
     最后，把所有样本的上述计算值求和，就是对 $\omega$ 的第 $j$ 个元素求偏导的结果。

   • 公式变形

     教材中给出的上述公式，每次只能计算梯度中的一个值，如果想一次性计算出整个梯度的值，要如何做呢？

     将上述计算过程中的 $x_i$ 变成 $\chi$ ，则 $z$ 的值为：
     $$z=\chi \cdot \omega =\left(\begin{array}{cccccc}{x}_1^{(1)}& x_1^{(2)}& x_1^{(3)}& \cdots & x_1^{(n)}& 1\\ x_2^{(1)}& x_2^{(2)}& x_2^{(3)}& \cdots & x_2^{(n)}& 1\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ x_N^{(1)}& x_N^{(2)}& x_N^{(3)}& \cdots & x_N^{(n)}& 1\end{array}\right)\left(\begin{array}{c}{\omega }^{(1)}\\ {\omega }^{(2)}\\ ⋮\\ {\omega }^{(n)}\\ b\end{array}\right)=\left(\begin{array}{c}{\omega }^{(1)}{x}_1^{(1)}+{\omega }^{(2)}{x}_1^{(2)}+\cdots +{\omega }^{(n)}{x}_1^{(n)}+b\\ {\omega }^{(1)}{x}_2^{(1)}+{\omega }^{(2)}{x}_2^{(2)}+\cdots +{\omega }^{(n)}{x}_2^{(n)}+b\\ ⋮\\ {\omega }^{(1)}{x}_N^{(1)}+{\omega }^{(2)}{x}_N^{(2)}+\cdots +{\omega }^{(n)}{x}_N^{(n)}+b\end{array}\right)=\left(\begin{array}{c}{x}_1\cdot \omega \\ x_2\cdot \omega \\ ⋮\\ x_i\cdot \omega \\ ⋮\\ x_N\cdot \omega \end{array}\right)=\left(\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_i\\ ⋮\\ z_N\end{array}\right)$$

     z = np.dot(X, w.T)  # 将点乘转换为矩阵乘法
     

     使用Sigmoid函数对 $z$ 进行变换：

     $$\hat{z}=\frac{\exp (z)}{1+\exp (z)}$$

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

     和 ${\rm Y}$ 做差：

     $${\rm Y}-\hat{z}=\left(\begin{array}{c}{y}_1-{\hat{z}}_1\\ y_2-{\hat{z}}_2\\ ⋮\\ y_i-{\hat{z}}_i\\ ⋮\\ y_N-{\hat{z}}_N\end{array}\right)$$

     y - z
     

     再和 $\chi$ 相乘：
     $$\left(\begin{array}{cccccc}{x}_1^{(1)}& \cdots & x_1^{(j)}& \cdots & x_1^{(n)}& 1\\ x_2^{(1)}& \cdots & x_2^{(j)}& \cdots & x_2^{(n)}& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ x_i^{(1)}& \cdots & x_i^{(j)}& \cdots & x_i^{(n)}& 1\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ x_N^{(1)}& \cdots & x_N^{(j)}& \cdots & x_N^{(n)}& 1\end{array}\right)\times \left(\begin{array}{c}{y}_1-{\hat{z}}_1\\ y_2-{\hat{z}}_2\\ ⋮\\ y_i-{\hat{z}}_i\\ ⋮\\ y_N-{\hat{z}}_N\end{array}\right)=\left(\begin{array}{cccccc}{x}_1^{(1)}(y_1-{\hat{z}}_1)& \cdots & x_1^{(j)}(y_1-{\hat{z}}_1)& \cdots & x_1^{(n)}(y_1-{\hat{z}}_1)& (y_1-{\hat{z}}_1)\\ x_2^{(1)}(y_2-{\hat{z}}_2)& \cdots & x_2^{(j)}(y_2-{\hat{z}}_2)& \cdots & x_2^{(n)}(y_2-{\hat{z}}_2)& (y_2-{\hat{z}}_2)\\ ⋮& \cdots & ⋮& \ddots & ⋮& ⋮\\ x_i^{(1)}(y_i-{\hat{z}}_i)& \cdots & x_i^{(j)}(y_i-{\hat{z}}_i)& \cdots & x_i^{(n)}(y_i-{\hat{z}}_i)& (y_i-{\hat{z}}_i)\\ ⋮& \cdots & ⋮& \ddots & ⋮& ⋮\\ x_N^{(1)}(y_N-{\hat{z}}_N)& \cdots & x_N^{(j)}(y_N-{\hat{z}}_N)& \cdots & x_N^{(n)}(y_N-{\hat{z}}_N)& (y_N-{\hat{z}}_N)\end{array}\right)$$
     需要注意的是，这里既不是点乘，也不是矩阵乘法，而是将 $ \Upsilon - \hat z$ 中的元素，与 $\chi$ 每一列元素对应相乘。之所以要这样做，是因为我们最终想得到的就是：
     $$x_i^{(j)}\left(y_i-{\hat{z}}_i\right)$$
     得益于Python，我们可以直接写作：

     X * (y - z)
     

     最后，我们需要对所有样本的、第 $j$ 个特征的计算值进行求和，得到的就是一个由偏导组成的向量，即梯度：
     $$\left(\frac{\partial L(\omega )}{\partial {\omega }^{(1)}}\cdots \frac{\partial L(\omega )}{\partial {\omega }^{(n+1)}}\right)=\left(\begin{array}{ccccc}{\sum }_{i=1}^N\left[x_i^{(1)}(y_i-{\hat{z}}_i)\right]& {\sum }_{i=1}^N\left[x_i^{(2)}(y_i-{\hat{z}}_i)\right]& \cdots & {\sum }_{i=1}^N\left[x_i^{(n)}(y_i-{\hat{z}}_i)\right]& {\sum }_{i=1}^N\left[(y_i-{\hat{z}}_i)\right]\end{array}\right)$$
     这里的代码表示为：

     (X * (y - z)).sum(axis=0)
     

     最核心的梯度部分梳理完了，将代码整理如下：

     def sigmoid(x):
         return 1/(1 + np.exp(-x))
     z = np.dot(X, w.T)
     # 求梯度（方向导数）
     grad = (X * (y - sigmoid(z))).sum(axis=0)
     

4. 参数更新

   根据已经计算出来的梯度，更新参数 $\omega$ ：
   $$\hat{\omega }=\left(\begin{array}{ccccc}{\omega }^{(1)}& {\omega }^{(2)}& \cdots & {\omega }^{(n)}& b\end{array}\right)+\eta \left(\frac{\partial L(\omega )}{\partial {\omega }^{(1)}}\cdots \frac{\partial L(\omega )}{\partial {\omega }^{(n+1)}}\right)=\left({\omega }^{(1)}+\eta \frac{\partial L(\omega )}{\partial {\omega }^{(1)}}\cdots {\omega }^{(n)}+\eta \frac{\partial L(\omega )}{\partial {\omega }^{(n)}}b+\eta \frac{\partial L(\omega )}{\partial {\omega }^{(n+1)}}\right)$$

   w += eta * grad
   

5. 停止条件

   模型中设置了两个停止条件，只要满足任一个，模型都会停止迭代。

   • 当梯度小于给定的阈值时，意味着接下来的参数更新只能带来很少的收益，也就意味着参数达到了我们认可的一种最优状态。在这里，只有当梯度中的每一个元素都小于给定的阈值时，我们才停止更新：

     if (np.abs(grad) <= delta).all():
         print('停止迭代')
     

   • 设定最大迭代次数，是为了防止模型的过拟合。所以当迭代次数达到提前设定的最大值时，也要停止更新：

     loop = 1
     while loop <= max_iter:
         
         pass  # 模型训练过程
     
         loop += 1
     

6. 模型预测

   求解出参数 $\omega$ 的值，根据模型的定义，我们就能预测新样本所属的标签：
   $$\begin{gathered} P(Y=1|x)=\frac{\exp (\omega \cdot x)}{1+\exp (\omega \cdot x)} \\ P(Y=0|x)=\frac{1}{1+\exp (\omega \cdot x)} \end{gathered}$$

   X_test = preprocessing(X_test)
   p = sigmoid(np.dot(X_test, w.T))
   p[np.where(p>=0.5)] = 1
   p[np.where(p<0.5)] = 0
   

5 代码封装

class LogisticRegression(object):
    """
    Binomial Logistic Regression classifier.
    
    Parameters
    ----------
    eta : float, default=0.1
        Step size for each iteration.
    
    max_iter : int, default=10000
        Maximum number of iterations taken for the solvers to converge.
        
    delta : float, default=1e-2
        deltaerance for stopping criteria.
    
    method : {'BGD', 'SGD'}, default='BGD'
        The way of Gradient Descent. 
        The 'BGD' is Batch Gradient Descent. The 'SGD' is Stochastic Gradient Descent.
        
    Attributes
    ----------
    w : ndarray of shape (1, n_features)
        Coefficient of the features in the model.
        
    """
    def __init__(self, 
                 eta=0.1, 
                 max_iter=10000, 
                 delta=1e-2, 
                 method='BGD'):
        self.eta = eta
        self.max_iter = max_iter
        self.delta = delta
        self.method = method
        self.w = None
        
    # Y = 1 时的模型
    def sigmoid(self, x):
        """
        Sigmoid function.
        
        Parameters
        ----------
        x : float or ndarray of shape (n_samples, )
            The independent variable of the Sigmoid function.
            
        Returns
        -------
        y : float or ndarray of shape (n_samples, )
            Function value
        """
        return 1/(1 + np.exp(-x))
    
    def preprocessing(self, X):
        """
        Extend input vector.
        
        Parameters
        ----------
        X : ndarray of shape (n_samples, n_features)
            Input vector.
            
        Returns
        -------
        X_new : ndarray of shape (n_samples, n_features + 1)
            Extend input vector.    
        """
        X_plus = np.ones(X.shape[0]).reshape(-1, 1)
        X_new = np.hstack([X, X_plus])
        return X_new
    
    def fit(self, X, y):
        """
        Fit the model according to the given training data.

        Parameters
        ----------
        X : ndarray of shape (n_samples, n_features)
            Training vector, where n_samples is the number of samples and
            n_features is the number of features.

        y : ndarray of shape (n_samples, )
            Target vector relative to X.

        Returns
        -------
        self
            Fitted estimator.
        """
        # 初始化 w
        X = self.preprocessing(X)
        self.w = np.zeros(X.shape[1]).reshape(1, -1)
        
        if self.method == 'BGD':
            loop = 1
            while loop <= self.max_iter:
                # w * x
                z = np.dot(X, self.w.T)
                # 求梯度（方向导数）
                grad = (X * (y - self.sigmoid(z))).sum(axis=0)
                # 下降（上升）的幅度
                if (np.abs(grad) <= self.delta).all():
                    break
                else:
                    self.w += self.eta * grad
                loop += 1
        elif self.method == 'SGD':
            # 随机获取样本点
            index = random.randint(0, X.shape[0] - 1)
            xi = X[index]
            yi = y[index]
            
            loop = 1
            while loop <= self.max_iter:
                z = np.dot(xi, self.w.T)
                grad = xi * (yi - self.sigmoid(z))
                if (np.abs(grad) <= self.delta).all():
                    break
                else:
                    self.w += self.eta * grad
                loop += 1
        else:
            pass
    
    def predict(self, X):
        """
        Predict class labels for samples in X.

        Parameters
        ----------
        X : ndarray of shape (n_samples, n_features)
            Samples.

        Returns
        -------
        p : ndarray of shape (n_samples,)
            Predicted class label per sample.
        """
        X = self.preprocessing(X)
        p = self.sigmoid(np.dot(X, self.w.T))
        p[np.where(p>=0.5)] = 1
        p[np.where(p<0.5)] = 0
        
        return p

---

1. 文章作者整理的内容，原书中并没有这一部分的表述 ↩︎