chapter5-逻辑斯谛回归-LogisticRegression

Posted on | In | | Visitors: ℃
Words count in article: 3.2k | Reading time ≈ 13

chapter5-逻辑斯谛回归-LogisticRegression

statistics learning

1.浅谈线性回归

由于前面最小二乘法和梯度下降算法，已经大力讨论了回归模型，因此，本章只进行简单的回顾回归模型。分析线性回归原理和与最小二乘法之间的区别。

线性回归与最小二乘法的最大区别，就在于损失函数的迭代。也就是如何优化损失函数。最小二乘法顾名思义，就是采用最小二乘法进行迭代，损失函数如下：
$$
J(\theta_0,\theta_1,\cdots,\theta_j) =\sum_{i=1}^m(h_θ(x_0^{(i)},x_1^{(i)},\cdots,x_j^{(i)}) - y^{(j)})^2=\sum_{i=1}^n(\sum_{j=1}^m\theta_jx_j^{(i)}-y^{(i)})^2
$$

对损失函数求导为，更新参数$\theta$：
$$
\theta = (X^TX)^-1X^TY
$$
对于线性回归算法来说，就是利用梯度下降算法，损失函数如下：
$$
J(\theta_0,\theta_1,\cdots,\theta_j) =\frac{1}{2m}\sum_{i=1}^m(h_θ(x_0^{(i)},x_1^{(i)},\cdots,x_j^{(i)}) - y^{(j)})^2=\frac{1}{2m}\sum_{i=1}^n(\sum_{j=1}^m\theta_jx_j^{(i)}-y^{(i)})^2
$$

对损失函数求偏导，更新参数$\theta$：
$$
\theta = \theta-\alpha X^T(\theta X-Y)
$$

综合上述，对于线性回归来说，利用最小二乘法与梯度下降算法去优化模型，区别就在于求全导(每一个错误样本都会更新)，还是偏导上(逐渐逼近最佳拟合函数)。相对于不同的数据量来说，小样本，最小二乘法会更快的收敛，时间会更少。当样本量较大时，最小二乘法的计算量会逐渐增大，这时，没有梯度下降算法收敛得快，但不一定是最优解。

2. LogisticRegression theory
1) 线性回归和逻辑回归的关系

​ 逻辑斯谛回归是在线性回归的基础之上，加入了非线性函数，就是sigmod函数，进行分类或者多分类。 因此logisticRegression是分类模型，那为什么带有回归两字呢？因为线性回归的y是连续的，而分类问题，y是离散的。那怎么做才能让连续分类问题转化为离散分类问题呢？利用sigmod函数，可以在闭区间$[0,1]$中的连续值，转化为0和1标签，多分类同理。sigmod函数如下：
$$
g(z) = \frac{1}{1+e^{-z}}
$$
可视化如图所示：

它有一个非常好的性质，即当z趋于正无穷时，g(z)趋于1，而当z趋于负无穷时，g(z)趋于0，这非常适合于我们的分类概率模型。另外，它还有一个很好的导数性质：　
$$
g^{‘}(z) = g(z)(1-g(z))
$$
这个通过函数对g(z)求导很容易得到，后面我们会用到这个式子。

如果我们令g(z)中的z为：${z = \theta x}$，这样就得到了二元逻辑回归模型的一般形式：
$$
h_{\theta}(x) = \frac{1}{1+e^{-\theta x}}
$$
　　其中x为样本输入，$h_{\theta}(x)$为模型输出，可以理解为某一分类的概率大小。而$\theta$为分类模型的要求出的模型参数。对于模型输出$h_{\theta}(x)$，让它和二元样本输出y（假设为0和1）有这样的对应关系，如果$h_{\theta}(x)$ >0.5 ，即$\theta x > 0$, 则y为1。如果$h_{\theta}(x) < 0.5$，即$x\theta < 0$, 则y为0。$y=0.5$是临界情况，此时$\theta x = 0$为， 从逻辑回归模型本身无法确定分类。

　 $h_{\theta}(x)$的值越小，而分类为0的的概率越高，反之，值越大的话分类为1的的概率越高。如果靠近临界点，则分类准确率会下降。

　此处我们也可以将模型写成矩阵模式：
$$
h_{\theta}(X) = \frac{1}{1+e^{-X\theta}}
$$
　其中$h_{\theta}(X)$为模型输出，为$ m \times 1$的维度。$X$为样本特征矩阵，为$m\times n$的维度。$\theta$为分类的模型系数，为$n\times1$的向量。

　理解了二元分类回归的模型，接着我们就要看模型的损失函数了，我们的目标是极小化损失函数来得到对应的模型系数$\theta$。

2) 二元逻辑回归的损失函数

回顾下线性回归的损失函数，由于线性回归是连续的，所以可以使用模型误差的的平方和来定义损失函数。但是逻辑回归不是连续的，自然线性回归损失函数定义的经验就用不上了。不过我们可以用最大似然法来推导出我们的损失函数。

我们知道，按照第二节二元逻辑回归的定义，假设我们的样本输出是0或者1两类。那么我们有
$$
P(y=1|x,\theta ) = h_{\theta}(x)\
P(y=0|x,\theta ) = 1- h_{\theta}(x)
$$

把这两个式子写成一个式子(服从伯努利分布)，就是：
$$
P(y|x,\theta ) = h_{\theta}(x)^y(1-h_{\theta}(x))^{1-y}
$$
其中y的取值只能是0或者1。

得到了y的概率分布函数表达式，我们就可以用似然函数最大化来求解我们需要的模型系数$\theta$。

为了方便求解，这里我们用对数似然函数最大化，对数似然函数取反即为我们的损失函数$J\theta$。其中：

似然函数的代数表达式为：
$$
L(\theta) = \prod\limits_{i=1}^{m}(h_{\theta}(x^{(i)}))^{y^{(i)}}(1-h_{\theta}(x^{(i)}))^{1-y^{(i)}}
$$
其中m为样本的个数。

对似然函数对数化取反的表达式，即损失函数表达式为：
$$
J(\theta) = -lnL(\theta) = -\sum\limits_{i=1}^{m}(y^{(i)}log(h_{\theta}(x^{(i)}))+ (1-y^{(i)})log(1-h_{\theta}(x^{(i)})))
$$
损失函数用矩阵法表达更加简洁：
$$
J(\theta) = -Y^Tlogh_{\theta}(X) - (E-Y)^T log(E-h_{\theta}(X))
$$
其中E为全1向量。

3) 二元逻辑回归损失函数的优化

对于二元逻辑回归的损失函数极小化，有比较多的方法，最常见的有梯度下降法，坐标轴下降法，牛顿法等。这里推导出梯度下降法中$\theta$每次迭代的公式。由于代数法推导比较的繁琐，我习惯于用矩阵法来做损失函数的优化过程，这里给出矩阵法推导二元逻辑回归梯度的过程。

对于$J(\theta) = -Y^T logh_{\theta}(X) - (E-Y)^T log(E-h_{\theta}(X))$，我们用$J(\theta)$对$\theta$向量求导可得：
$$
\frac{\partial}{\partial\theta}J(\theta) = X^T[\frac{1}{h_{\theta}(X)}\odot h_{\theta}(X)\odot (E-h_{\theta}(X))\odot (-Y)] \ \qquad + X^T[\frac{1}{E-h_{\theta}(X)}\odot h_{\theta}(X)\odot (E-h_{\theta}(X))\odot (E-Y)]
$$
这一步我们用到了向量求导的链式法则，和下面三个基础求导公式的矩阵形式：
$$
\frac{\partial}{\partial x}logx = 1/x \ \
\frac{\partial}{\partial z}g(z) = g(z)(1-g(z)) \quad (g(z)为sigmoid函数) \
\frac{\partial x\theta}{\partial \theta} = x
$$
对于刚才的求导公式我们进行化简可得：
$$
\frac{\partial}{\partial\theta}J(\theta) = X^T(h_{\theta}(X) - Y )
$$
从而在梯度下降法中每一步向量$\theta$的迭代公式如下：
$$
\theta = \theta - \alpha X^T(h_{\theta}(X) - Y )
$$
其中，$\alpha$为梯度下降法的步长。

实践中，我们一般不用操心优化方法，大部分机器学习库都内置了各种逻辑回归的优化方法，不过了解至少一种优化方法还是有必要的。

4) 二元逻辑回归的正则化

　逻辑回归也会面临过拟合问题，所以我们也要考虑正则化。常见的有$L1$正则化和$L2$正则化。

　逻辑回归的$L1$正则化的损失函数表达式如下，相比普通的逻辑回归损失函数，增加了$L1$的范数做作为惩罚，超参数$\alpha$作为惩罚系数，调节惩罚项的大小。

二元逻辑回归的$L1$正则化损失函数表达式如下：
$$
J(\theta) = -Y^T\bullet logh_{\theta}(X) - (E-Y)^T\bullet log(E-h_{\theta}(X)) +\alpha ||\theta||_1
$$
其中$||\theta||_1$为$\theta$的$L1$范数。

逻辑回归的$L1$正则化损失函数的优化方法常用的有坐标轴下降法和最小角回归法。

二元逻辑回归的L2正则化损失函数表达式如下：　　
$$
　J(\theta) = -Y^T\bullet logh_{\theta}(X) - (E-Y)^T\bullet log(E-h_{\theta}(X)) + \frac{1}{2}\alpha||\theta||_2^2
$$
其中$||\theta||_2$为$\theta$的$L2$范数。

逻辑回归的$L2$正则化损失函数的优化方法和普通的逻辑回归类似。

5） 多元逻辑回归

前面几节我们的逻辑回归的模型和损失函数都局限于二元逻辑回归，实际上二元逻辑回归的模型和损失函数很容易推广到多元逻辑回归。比如总是认为某种类型为正值，其余为0值，这种方法为最常用的$one-vs-rest$，简称$OvR$.

另一种多元逻辑回归的方法是$Many-vs-Many(MvM)$，它会选择一部分类别的样本和另一部分类别的样本来做逻辑回归二分类。最常用的是$One-Vs-One（OvO）$。$OvO$是$MvM$的特例。每次我们选择两类样本来做二元逻辑回归。

这里只介绍多元逻辑回归的softmax回归的一种特例推导：

首先回顾下二元逻辑回归

$$
P(y=1|x,\theta ) = h_{\theta}(x) = \frac{1}{1+e^{-x\theta}} = \frac{e^{x\theta}}{1+e^{x\theta}}\

　　　　P(y=0|x,\theta ) = 1- h_{\theta}(x) = \frac{1}{1+e^{x\theta}}
$$

其中y只能取到0和1。则有：
$$
ln\frac{P(y=1|x,\theta )}{P(y=0|x,\theta)} = \theta x
$$

如果我们要推广到多元逻辑回归，则模型要稍微做下扩展。

我们假设是K元分类模型,即样本输出y的取值为$1，2，\dots，K$。

根据二元逻辑回归的经验，我们有：
$$
ln\frac{P(y=1|x,\theta )}{P(y=K|x,\theta)} = \theta_1 x \

ln\frac{P(y=2|x,\theta )}{P(y=K|x,\theta)} = \theta_2 x　\

…\

ln\frac{P(y=K-1|x,\theta )}{P(y=K|x,\theta)} = \theta_{K-1} x\
$$

上面有$K-1$个方程。

加上概率之和为1的方程如下：
$$
\sum\limits_{i=1}^{K}P(y=i|x,\theta ) = 1
$$
从而得到K个方程，里面有K个逻辑回归的概率分布。

解出这个K元一次方程组，得到K元逻辑回归的概率分布如下：
$$
P(y=k|x,\theta ) = e^{x\theta_k} \bigg/ 1+\sum\limits_{t=1}^{K-1}e^{x\theta_t}　 k = 1,2,…K-1 \
P(y=K|x,\theta ) = 1 \bigg/ 1+\sum\limits_{t=1}^{K-1}e^{x\theta_t}
$$
多元逻辑回归的损失函数推导以及优化方法和二元逻辑回归类似，这里就不累述。

meaching learning

a. 手动创建模型
  1. 导入相关工具包
    1
    2
    3
    4
    5
    6
    7
    8
    from math import exp
    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    %matplotlib inline

    from sklearn.datasets import load_iris
    from sklearn.model_selection import train_test_split
  1. 导入数据
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    # data
    def create_data():
        iris = load_iris()
        df = pd.DataFrame(iris.data, columns=iris.feature_names)
        df['label'] = iris.target
        df.columns = ['sepal length', 'sepal width', 'petal length', 'petal width', 'label']
        data = np.array(df.iloc[:100, [0,1,-1]])
        # print(data)
        return data[:,:2], data[:,-1]
    #split train and test set data    
    X, y = create_data()
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
  2. 创建模型
    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
    class LogisticReressionClassifier:
        def __init__(self, max_iter=200, learning_rate=0.01):
            self.max_iter = max_iter
            self.learning_rate = learning_rate

        def sigmoid(self, x):
            return 1 / (1 + exp(-x))

        def data_matrix(self, X):
            data_mat = []
            for d in X:
                data_mat.append([1.0, *d])
            return data_mat

        def fit(self, X, y):
            # label = np.mat(y)
            data_mat = self.data_matrix(X)  # m*n
            self.weights = np.zeros((len(data_mat[0]), 1), dtype=np.float32)

            for iter_ in range(self.max_iter):
                for i in range(len(X)):
                    result = self.sigmoid(np.dot(data_mat[i], self.weights))
                    error = y[i] - result
                    self.weights += self.learning_rate * error * np.transpose(
                        [data_mat[i]])
            print('LogisticRegression Model(learning_rate={},max_iter={})'.format(
                self.learning_rate, self.max_iter))

        # def f(self, x):
        #     return -(self.weights[0] + self.weights[1] * x) / self.weights[2]

        def score(self, X_test, y_test):
            right = 0
            X_test = self.data_matrix(X_test)
            for x, y in zip(X_test, y_test):
                result = np.dot(x, self.weights)
                if (result > 0 and y == 1) or (result < 0 and y == 0):
                    right += 1
            return right / len(X_test)
  3. 训练模型与可视化
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    lr_clf = LogisticReressionClassifier()
    lr_clf.fit(X_train, y_train)

    lr_clf.score(X_test, y_test)

    x_ponits = np.arange(4, 8)
    y_ = -(lr_clf.weights[1]*x_ponits + lr_clf.weights[0])/lr_clf.weights[2]
    plt.plot(x_ponits, y_)

    #lr_clf.show_graph()
    plt.scatter(X[:50,0],X[:50,1], label='0')
    plt.scatter(X[50:,0],X[50:,1], label='1')
    plt.legend()

    结果如图：

b. sklearn 工具包
  1. 创建与训练
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    from sklearn.linear_model import LogisticRegression 
    #创建模型
    clf = LogisticRegression(max_iter=200)
    #模型训练
    clf.fit(X_train, y_train)
    #得到分数
    clf.score(X_test, y_test) #1.0
    #输出参数
    print(clf.coef_, clf.intercept_) #[[ 1.92439982 -3.1874392 ]] [-0.57523308]

    x_ponits = np.arange(4, 8)
    y_ = -(clf.coef_[0][0]*x_ponits + clf.intercept_)/clf.coef_[0][1]
    plt.plot(x_ponits, y_)

    plt.plot(X[:50, 0], X[:50, 1], 'bo', color='blue', label='0')
    plt.plot(X[50:, 0], X[50:, 1], 'bo', color='orange', label='1')
    plt.xlabel('sepal length')
    plt.ylabel('sepal width')
    plt.legend()

    结果如下：

参考链接

逻辑回归原理

什么是概率分布函数

概率分布函数与似然函数的关系