机器学习入门 07 逻辑回归

1 什么是逻辑回归（Logistic Regression）

1.1 介绍

​ 一般用来解决分类问题，只能解决二分类问题。

​ 将样本的特征和样本的概率联系起来，概率是一个数，所以可以叫做回归问题。

​ 在多项式回归中，$\hat{y}=f(x)={\theta }^{⊺}\cdot x_b$，${\theta }^{⊺}$是系数，$x_b$是添加了$x_b\equiv 1$的向量（矩阵）。$\hat{y}$的值域为$(-\infty ,\infty )$，而概率$\hat{p}$的值域为$[0,1]$，这里使用Sigmoid函数将y的值域范围转换为p的值域范围。

​ Sigmoid函数：$\sigma (t)=\frac{1}{1+e^{-t}}$，$\{\begin{array}{c}t>0,p>0.5\\ t<0,p<0.5\end{array}$

​ 函数图像：

​ 令 $\hat{p}=\sigma ({\theta }^{⊺}\cdot x_b)=\frac{1}{1+e^{-{\theta }^{⊺}\cdot x_b}}$，最终分类 $\hat{y}=\{\begin{array}{c}1,\hat{p}\ge 0.5\\ 0,\hat{p}\le 0.5\end{array}$

​ 问题：

​ 对于给定的样本数据集X和y，如何找到参数$\theta$，使得用这样的方式，可以最大程度的获得样本数据集X对应的分类输出y。

2 逻辑回归的损失函数

$$\begin{gathered} cost=\{\begin{array}{c}-log(\hat{p})ify=1\\ -log(1-\hat{p})ify=0\end{array}\Rightarrow cost=-ylog(\hat{p})-(1-y)log(1-\hat{p}) \\ J(\theta )=-\frac{1}{m}{\sum }_{i=1}^m[y^{(i)}log({\hat{p}}^{(i)})+(1-y^{(i)})log(1-{\hat{p}}^{(i)})],{\hat{p}}^{(i)}=\sigma (X_b^{(i)}\theta )=\frac{1}{1+e^{-X_b^{(i)}\theta }} \end{gathered}$$

3 使用梯度下降法求$\theta$ 使得$J(\theta )$最小

$$\begin{gathered} J(\theta )=-\frac{1}{m}{\sum }_{i=1}^m[y^{(i)}log({\hat{p}}^{(i)})+(1-y^{(i)})log(1-{\hat{p}}^{(i)})],{\hat{p}}^{(i)}=\sigma (X_b^{(i)}\theta )=\frac{1}{1+e^{-X_b^{(i)}\theta }} \\ J(\theta )=-\frac{1}{m}{\sum }_{i=1}^m[y^{(i)}log(\sigma (X_b^{(i)}\theta ))+(1-y^{(i)})log(1-\sigma (X_b^{(i)}\theta ))] \\ \{\begin{array}{c}\sigma (t)=\frac{1}{1+e^{-t}}={(1+e^{-t})}^{-1},\sigma (t)\prime ={(1+e^{-t})}^{-2}{e}^{-t}\\ [log\sigma (t)]\prime =\frac{1}{\sigma (t)}\sigma (t)\prime =\frac{1}{{(1+e^{-t})}^{-1}}{(1+e^{-t})}^{-2}{e}^{-t}=\frac{e^{-t}}{1+e^{-t}}=\frac{1+e^{-t}-1}{1+e^{-t}}=1-\sigma (t)\\ [log(1-\sigma (t))]\prime =\frac{1}{1-\sigma (t)}(-1)\sigma (t)\prime =-\frac{1}{1-\sigma (t)}{(1+e^{-t})}^{-2}{e}^{-t}=-\frac{1}{\frac{1+e^{-t}}{1+e^{-t}}-\frac{1}{1+e^{-t}}}{(1+e^{-t})}^{-2}{e}^{-t}\\ =-\frac{1+e^{-t}}{e^{-t}}{(1+e^{-t})}^{-2}{e}^{-t}=-{(1+e^{-t})}^{-1}=-\sigma (t)\end{array} \\ \{\begin{array}{c}\frac{d(y^{(i)}log(\sigma (X_b^{(i)}\theta )))}{d{\theta }_j}=y^{(i)}(1-\sigma (X_b^{(i)}\theta ))X_j^{(i)}\\ \frac{d((1-y^{(i)})log(1-\sigma (X_b^{(i)}\theta )))}{d{\theta }_j}=(1-y^{(i)})(-\sigma (X_b^{(i)}\theta ))X_j^{(i)}\end{array}\begin{array}{c}+\\ ⟹\end{array}[y^{(i)}-\sigma (X_b^{(i)}\theta )]X_j^{(i)} \\ \Rightarrow \frac{dJ(\theta )}{d{\theta }_j}=\frac{1}{m}{\sum }_{i=1}^m(\sigma (X_b^{(i)}\theta )-y^{(i)})X_j^{(i)} \\ \Rightarrow \nabla J(\theta )=\left(\begin{array}{c}\partial J/\partial {\theta }_0\\ \partial J/\partial {\theta }_1\\ ⋮\\ \partial J/\partial {\theta }_n\end{array}\right)=\frac{1}{m}\left(\begin{array}{c}{\sum }_{i=1}^m(\sigma (X_b^{(i)}\theta )-y^{(i)})X_0^{(i)})\\ {\sum }_{i=1}^m(\sigma (X_b^{(i)}\theta )-y^{(i)})X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m(\sigma (X_b^{(i)}\theta )-y^{(i)})X_n^{(i)})\end{array}\right)=\frac{1}{m}\left(\begin{array}{c}{\sum }_{i=1}^m({\hat{p}}^{(i)}-y^{(i)})X_0^{(i)})\\ {\sum }_{i=1}^m({\hat{p}}^{(i)}-y^{(i)})X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m({\hat{p}}^{(i)}-y^{(i)})X_n^{(i)})\end{array}\right)(\Rightarrow \frac{1}{m}\left(\begin{array}{c}{\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_0^{(i)})\\ {\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_1^{(i)})\\ ⋮\\ {\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)})X_n^{(i)})\end{array}\right)) \\ =\frac{1}{m}{X}_b^{⊺}[\sigma (X_b\theta )-y] \end{gathered}$$

4 Python实现逻辑回归算法

逻辑回归模块：

import numpy as np
from sklearn.metrics import accuracy_score

class LogisticRegression:

    def __init__(self):
        """初始化Logistic Regression模型"""
        self.coef_ = None
        self.intercept_ = None
        self._theta = None

    def _sigmoid(self, t):
        return 1./(1.+np.exp(-t))

    def fit(self, X_train, y_train, eta=0.01, n_iters=1e4):
        """根据训练数据集X_train, y_train, 使用梯度下降法训练Logistic Regression模型"""
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of X_train must be equal to the size of y_train"

        def J(theta, X_b, y):
            """目标函数"""
            y_hat = self._sigmoid(X_b.dot(theta))
            try:
                return -np.sum(y*np.log(y_hat)+(1-y)*np.log(1-y_hat)) / len(y)
            except:
                return float('inf')

        def dJ(theta, X_b, y):
            """梯度"""
            return X_b.T.dot(self._sigmoid(X_b.dot(theta)) - y) / len(X_b)

        def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):

            theta = initial_theta
            cur_iter = 0

            while cur_iter < n_iters:
                gradient = dJ(theta, X_b, y)
                last_theta = theta
                theta = theta - eta * gradient
                if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
                    break

                cur_iter += 1

            return theta

        X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
        initial_theta = np.zeros(X_b.shape[1])
        self._theta = gradient_descent(
            X_b, y_train, initial_theta, eta, n_iters)

        self.intercept_ = self._theta[0]
        self.coef_ = self._theta[1:]

        return self

    def predict_proba(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果概率向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
        return self._sigmoid(X_b.dot(self._theta))

    def predict(self, X_predict):
        """给定待预测数据集X_predict，返回表示X_predict的结果向量"""
        assert self.intercept_ is not None and self.coef_ is not None, \
            "must fit before predict!"
        assert X_predict.shape[1] == len(self.coef_), \
            "the feature number of X_predict must be equal to X_train"

        proba=self.predict_proba(X_predict)
        return np.array(proba>=0.5,dtype='int')

    def score(self, X_test, y_test):
        """根据测试数据集 X_test 和 y_test 确定当前模型的准确度"""

        y_predict = self.predict(X_test)
        return accuracy_score(y_test, y_predict)

    def __repr__(self):
        return "LogisticRegression()"

使用上述模块：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
# 准备数据
iris=datasets.load_iris()

X=iris.data
y=iris.target

# 只选取其中两个分类
X=X[y<2,:2]
y=y[y<2]

X.shape # (100, 2)
# 绘图
plt.scatter(X[y==0,0],X[y==0,1],color='r')
plt.scatter(X[y==1,0],X[y==1,1],color='b')
plt.show()

from sklearn.model_selection import train_test_split
import LogisticRegression
# 分割数据集
X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=666)

log_reg=LogisticRegression.LogisticRegression()
log_reg.fit(X_train,y_train)
# 模型准确率
log_reg.score(X_test,y_test) # 1.0

5 决策边界

$$\begin{gathered} \{\begin{array}{c}\hat{p}=\sigma ({\theta }^{⊺}\cdot x_b)=\frac{1}{1+e^{-{\theta }^{⊺}\cdot x_b}},\hat{y}=\{\begin{array}{c}1,\hat{p}\ge 0.5\\ 0,\hat{p}<0.5\end{array}\\ p=\sigma (t)=\frac{1}{1+e^{-t}}\Rightarrow \begin{array}{c}t>0,p>0.5\\ t<0,p<0,5\end{array}\end{array} \\ \Rightarrow \hat{y}=\{\begin{array}{c}1,\hat{p}\ge 0.5,{\theta }^{⊺}\cdot x_b\ge 0\\ 0,\hat{p}<0.5,{\theta }^{⊺}\cdot x_b<0\end{array} \\ \Rightarrow {\theta }^{⊺}\cdot x_b=0 \end{gathered}$$

决策边界为: ${\theta }^{⊺}\cdot x_b=0$

若X有两个特征：${\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2=0\Rightarrow x_2=\frac{-{\theta }_0-{\theta }_1{x}_1}{{\theta }_2}$

绘图演示：

def x2(x1):
    return (-log_reg.coef_[0]*x1-log_reg.intercept_)/log_reg.coef_[1]
    
x1_plot=np.linspace(4,8,1000)
x2_plot=x2(x1_plot)

plt.scatter(X[y==0,0],X[y==0,1],color='r')
plt.scatter(X[y==1,0],X[y==1,1],color='b')
plt.plot(x1_plot,x2_plot)
plt.show()