统计学习方法（3）逻辑回归（Logistic Regression）

统计学习方法（3）逻辑回归

1、从线性回归到逻辑回归（模型）

1.1 线性模型：

给定数据集 { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x n , y n ) } \{(x_1, y_1), (x_2,y_2),...,(x_n,y_n)\} {(x1​,y1​),(x2​,y2​),...,(xn​,yn​)}，求参数$\omega$满足如下回归模型：
$$\hat{y}={\omega }_0+{\omega }_1{x}_1+{\omega }_2{x}_2+...+{\omega }_n{x}_n$$
线性回归模型的向量形式：
$$\hat{y}={\omega }^T\cdot \mathbf{x}$$
其中：${\theta }^T$表示向量$\theta$的转置（行向量变为了列向量）;
$\mathbf{x}$为每个样本中特征值的向量形式，包括 $x_1$，…， $x_n$，而且 $x_0$ =1。

模型求解：
基于最小二乘思想（方差最小）求解模型参数$\omega$：
$${\hat{\omega }}^{\ast }=argmin\frac{1}{m}\sum _{i=1}^m({\hat{\omega }}^T{x}^{(i)}-y(i){)}^2$$
写成向量形式：
$${\hat{\omega }}^{\ast }=argmin({\hat{\omega }}^T\mathbf{x}-\mathbf{y}{)}^{\mathbf{T}}({\hat{\omega }}^{\mathbf{T}}\mathbf{x}-\mathbf{y})$$
另$E_{\hat{\omega }}=(\hat{\omega }\mathbf{x}-\mathbf{y}{)}^{\mathbf{T}}(\hat{\omega }\mathbf{x}-\mathbf{y})$，对其求导：
$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=2X(X^T\omega -y)$$
令导数为0，求正态方程得到：
$$\hat{\omega }=(X^{T}X{)}^{-1}{X}^{T}y$$

注：当$X^{T}X$列数大于行数（特征数大于样例数），为非满秩矩阵，需引入正则项。

若无法求解正态方程或数据量较大，可采用梯度下降法：
对任意一点$x^{(i)}$的梯度：
$$\frac{\partial E_{\hat{\omega }}}{\partial \hat{\omega }}=2x^{(i)}({\omega }^T{x}^{(i)}-y)$$
随机梯度下降：
$${\omega }_{j+1}:={\omega }_j-\alpha \nabla f(x)={\omega }_j+\alpha (y^{(i)}-{\omega }^T{x}^{(i)})x^{(i)}$$

   \;

1.2 逻辑回归：

考虑二分类问题，由于线性回归输出值是 { − ∞ , + ∞ } \{-\infty,+\infty\} {−∞,+∞}，应选取合适的映射函数$g^{-}(\cdot )$，将模型输出映射为0/1值，这里想到单位阶跃函数可以做到，但是单位阶跃函数不连续，于是想到使用类似单位阶跃函数，又能单调可微的函数，所以可以使用对数几率函数，即sigmoid函数：
$$y=\frac{1}{1+e^{-z}}$$
sigmoid函数有一个非常好的性质，即当$z$趋于正无穷时，$g(z)$趋于1，而当$z$趋于负无穷时，$g(z)$趋于0，这非常适合于我们的分类概率模型。另外，它还有一个很好的导数性质：
$$g\prime (z)=g(z)(1-g(z))$$
这个通过函数对g(z)求导很容易得到，后面我们会用到这个式子。

线性模型带入sigmoid函数，得到：
$$h_{\omega }(x,b)=\frac{1}{1+e^{-({\omega }^T\mathbf{x}+b)}}$$
为了方便，将权重向量和输入向量加以扩充，于是得到逻辑回归模型：
$$h_{\omega }(\mathbf{x})=\frac{1}{1+{\mathbf{e}}^{-{\omega }^{\mathbf{T}}\mathbf{x}}}$$
其中$x$为样本输入，$h_{\omega }(\mathbf{x})$为模型输出，可以理解为某一分类的概率大小。而$\omega$为分类模型要求的模型参数。对于模型输出$h_{\omega }(\mathbf{x})$，我们让它和我们的二元样本输出$y$（假设为0和1）有这样的对应关系：如果$h_{\omega }(\mathbf{x})>0.5$ ，即${\omega }^T\mathbf{x}>0$, 则$y$为1；如果$h_{\omega }(\mathbf{x})<0.5$，即${\omega }^T\mathbf{x}<0$, 则$y$为0；$h_{\omega }(\mathbf{x})=0.5$是临界情况，此时${\omega }^T\mathbf{x}=0$， 从逻辑回归模型本身无法确定分类。

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

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

按照逻辑回归的定义，假设我们的样本输出是0或者1两类，属于伯努利分布（0-1分布），那么就有：
$$p(y=1|x)=h_{\omega }(\mathbf{x})=\frac{{\mathbf{e}}^{{\omega }^{\mathbf{T}}\mathbf{x}}}{1+{\mathbf{e}}^{{\omega }^{\mathbf{T}}\mathbf{x}}}$$$$p(y=0|x)=1-h_{\omega }(\mathbf{x})=\frac{1}{1+{\mathbf{e}}^{{\omega }^{\mathbf{T}}\mathbf{x}}}$$
其中，$\omega =({\omega }^{(1)},{\omega }^{(2)},...,{\omega }^{(n)},b{)}^T，x=(x^{(1)},x^{(2)},...,x^{(n)},1{)}^T$

---

2、策略

得到逻辑回归模型之后，考虑使用什么策略来学习模型。

2.1 极大似然估计

对于二分类，由于模型 y i ∈ { 0 , 1 } y_i\in \{0,1\} yi​∈{0,1}服从0-1分布，设：
$$P(Y=1|x)=\varphi (z)，P(Y=0|x)=1-\varphi (z)$$则概率分布可写为：
$$p(y_i|x_i;\hat{\omega })=\varphi (z_i{)}^{y_i}(1-\varphi (z_i){)}^{1-y_i}$$
由于概率分布已知，根据观测数据，可采用极大似然估计法来估计参数$\omega$和$b$。

通过概率分布写出联合概率密度：
$$L(\omega )=\prod _{i=1}^{n}p(y^{(i)}|\theta ,\omega )=\prod _{i=1}^n\varphi (z^{(i)}{)}^{y^{(i)}}(1-\varphi (z^{(i)}){)}^{(1-y^{(i)})}$$取对数可得：
$$l(\omega )=\sum _{i=1}^{n}lnp(y^{(i)}|\theta ,\omega )=\sum _{i=1}^n{y}^{(i)}ln[\varphi (z^{(i)})]+(1-y^{(i)})ln[(1-\varphi (z^{(i)}))]$$
对$l(\omega )$求极大值，可求出最有可能的$\omega$。
$$-l(\omega )=\sum _{i=1}^n-lnp(y^{(i)}|\theta ,\omega )=-\sum _{i=1}^n{y}^{(i)}ln[\varphi (z^{(i)})]+(1-y^{(i)})ln[(1-\varphi (z^{(i)}))]$$
最大化$l(\omega )$相当于最小化$-l(\omega )$，因此可以将$-l(\omega )$作为损失函数$J(\omega )=-l(\omega )$，求解参数。

2.2 损失函数

同其他统计学习方法一样，逻辑回归也可以使用损失函数来学习模型，逻辑回归使用的损失函数为对数损失函数：
$$L(Y,P(Y|X))=-logP(Y|X)$$
该模型的损失函数为：
$$P(\varphi (z),y|x;\omega )=\{\begin{array}{ll}-ln\varphi (z),& \text{if y=1}\\ -ln(1-\varphi (z)),& \text{if y=0}\end{array}$$
由损失函数得到经验风险函数：
$$J(\hat{\omega })=-\frac{1}{n}\sum _{i=1}^n(y_{i}ln\varphi (z^{(i)})+(1-y_i)ln(1-\varphi (z^{(i)})))$$

同极大似然法得到的对数似然函数比较，可以发现对数损失函数与极大似然估计的对数似然函数本质上是等价的，所以逻辑回归最小化对数损失函数求参数，实际上与采用极大化似然函数来求参数是一致的，实际逻辑回归的对数损失函数也就是用极大化似然估计得到的。

---

3、算法

由于对数似然函数是凸函数，故可采用数值优化算法如梯度下降法、拟牛顿法求其最优解。

以下采用梯度下降法：
令
$$J(\theta )=-l(\omega )=-\sum _{i=1}^n{y}^{(i)}ln[\varphi (z^{(i)})]+(1-y^{(i)})ln[(1-\varphi (z^{(i)}))]$$
$$\frac{\partial J(\theta )}{\partial {\theta }_j}=-\sum _{i=1}^n[y^{(i)}\frac{1}{\varphi (z^{(i)})}-(1-y^{(i)})\frac{1}{(1-\varphi (z^{(i)}))}]\frac{\partial \varphi (z^{(i)})}{\partial {\theta }_j}$$
由于对于sigmoid函数：
$${\varphi }^{\prime }(z)=\varphi (z)(1-\varphi (z))$$
$$\frac{\partial ({\omega }^{T}x+b)}{\partial \omega }=x$$
故：
$$=-\sum _{i=1}^n[y^{(i)}\frac{1}{\varphi (z^{(i)})}-(1-y^{(i)})\frac{1}{(1-\varphi (z^{(i)}))}]\varphi (z^{(i)})(1-\varphi (z^{(i)}))\frac{\partial z^{(i)}}{\partial {\theta }_j}$$
$$=-\sum _{i=1}^n[y^{(i)}(1-\varphi (z^{(i)}))-(1-y^{(i)})\varphi (z^{(i)})]x_j^{(i)}$$
$$=-\sum _{i=1}^n(y^{(i)}-\varphi (z^{(i)}))x_j^{(i)}$$
$$=-\sum _{i=1}^n(y^{(i)}-f(x))x_j^{(i)}$$
从而得到批量梯度下降：
$${\theta }_{j+1}:={\theta }_j+\eta \sum _{i=1}^n(y^{(i)}-f(x))x_j^{(i)}$$
随机梯度下降：
$${\theta }_{j+1}:={\theta }_j+\eta (y^{(i)}-f(x))x_j^{(i)}，foriinrange(n)$$

---

4、实现

批量梯度下降：

    def train(self, X, y):
        """
        批量梯度下降：训练模型

        :param X: 输入集，n*m(n行m列)
        :param y: 标签集，1*m(1行m列)
        :return: 训练好的模型
        """
        n, m = np.shape(X)
        # X = (X, 1), w=(w;b)
        # w^T*x+b => w^T*X
        X = np.mat(np.concatenate((X, np.ones((n, 1))), axis=1))
        y = np.mat(y).T
        w = np.ones((m + 1, 1))

        alpha = self.alpha
        tol = self.tol
        records = []

        for i in range(self.max_iter):
            # 实际值与预测值之差
            error = y - 1.0 / (1 + np.exp(-X * w))
            # 误差变化记录
            records.append((error.T * error)[0, 0])
            if records[-1] <= tol:
                break
            # 负梯度方向更新参数
            w += alpha * X.T * error

        self.w = w[:-1]
        self.b = w[-1]

随机梯度下降：

    def train(self, X, y):
        """
        随机梯度下降

        :param X: 训练集
        :param y: 标签
        :return: w,b
        """
        n, m = np.shape(X)
        w = np.ones((m + 1, 1))

        alpha = self.alpha
        records = []

        for i in range(self.max_iter):
            row = random.randint(0, n - 1)
            x_i = np.mat(np.concatenate((X[row], np.ones(1))))
            y_i = y[row]
            error = y_i - 1.0 / (1 + np.exp(-np.matmul(x_i, w))[0, 0])
            records.append(error * error)

            w += alpha * error * x_i.T

        self.w = w[:-1]
        self.b = w[-1]

博客：
https://www.cnblogs.com/ModifyRong/p/7739955.html