机器学习 - 逻辑回归

70

(本系列的 Jupyter Notebook 可以到我的码云下载。)

简介

逻辑回归是用来解决分类任务的,比如判断邮件是否垃圾邮件,肿瘤是否良性等。与线性回归一样,逻辑回归的输入也是一系列特征值,但是输出为 1 和 0 的分类。比如判断邮件是否垃圾邮件时,输入为邮件的关键词,输出为 1 表示垃圾邮件,0 表示非垃圾邮件;判断肿瘤是否良性时,输入的是病人的年龄和肿瘤大小,输出 1 表示肿瘤为良性,输出 0 表示肿瘤为恶性。本篇从零开始讲解逻辑回归。

建模

首先来讨论如何建立模型,模型是用来预测的,而实际在预测一个事件是否发生时,我们都用概率描述它。比如天气预报员说明天有雨,意思是说明天有很大概率有雨,而不是说明天一定有雨。再比如说2018世界杯巴西会夺冠,意思是说巴西会有很大概率夺冠,而不是说巴西一定会夺冠。因此用概率作为模型的输出是比较好的。在实际预测时,可以指定一个阈值(通常为 0.5),如果模型输出的概率大于等于这个阈值,那么就预测为正类(用1表示),如果小于这个阈值,那么就预测为负类(用0表示)。比如天气预报模型输出明天有雨的概率为 0.7,那么就预测为明天有雨,如果输出明天有雨的概率为 0.3,那么就预测为明天没有雨。

综上,我们的模型如下:


\begin{equation} h_\theta(\mathbf{x}) = p(\mathbf{x}, \boldsymbol{\theta}) \tag{1} \end{equation}

其中 $h$ 表示假设函数,$\boldsymbol{\theta}$ 表示模型参数,$\mathbf{x}$ 表示样本的特征向量(即模型输入),$p(\mathbf{x}, \boldsymbol{\theta})$ 表示输出正类的概率。

而我们的预测函数(输出函数)为:


\begin{equation} y = \left\{ \begin{array}{ll} 1, & p \geqslant 0.5 \\ 0, & p \lt 0.5 \end{array} \right. \tag{2} \end{equation}

其中 $y$ 表示预测值,$p$ 表示输出正类的概率。

求解

到这里我们建立了假设函数以及预测函数,但是还有两个问题没有解决:

  • 问题一:$(1)$ 式中 $p$ 是一个关于样本特征 $\mathbf{x}$ 和模型参数 $\boldsymbol{\theta}$ 的函数,这个函数的输出是一个概率,那么如何通过 $\mathbf{x}$ 和 $\boldsymbol{\theta}$ 计算这个概率呢?也就是说如何表示 $p$ 这个函数呢?
  • 问题二:知道了 $p$ 以后,如何计算模型参数 $\boldsymbol{\theta}$ 呢?

首先让我们来讨论问题一。

假设有如下的两类样本:

图1 - 样本

可以很容易用一条直线将它们分开:

图2 - 分隔线

不过请先等一等,我们并不是要把它们分开,而是要求出它们是正类的概率。下面我们标记出一些有代表性的点,然后来看看它们是正类的概率应该是多少:

图3 - 概率分析

先看 $A$ 点,在分割线上的点是比较特殊的,它们既可能是正类,也可能是负类,因此它们的概率应该在 $0.5$ 左右(既可以大于等于 $0.5$,也可以小于 $0.5$),我们可以假设 $A$ 点的概率就是 $0.5$,即 $P\{A\} = 0.5$。 再来看 $B$、 $C$ 两个点,这两个点正在远离分割线,因此它们是正类的概率应该大于 $0.5$,假设它们的概率为 $P\{B\} = 0.6$,$P\{C\} = 0.7$。接着看 $D$、$E$ 两个点,可以看到它们已经离分割线非常远了,所以它们的概率应该都很大,另外它们的间距与 $B$、$C$ 两个点的间距差不多一样,那么是不是它们概率的差值也应该与 $B$、$C$ 概率的差值一样呢?答案是否定的,因为它们离分割线实在太远了,已经完全不受分割线的影响,因此它们的概率应该很接近,举个通俗的例子,比如日本发生了地震,那么北京和天津能感受到房屋晃动的概率应该都差不多一样(都很小) ,尽管天津离日本更近一些。因此离分割线越远,概率的增长或降低速度越小。这里假设 $D$、$E$ 的概率为 $P\{D\} = 0.9$、$P\{E\} = 0.93$。上面的讨论同样适用于 $F$、$G$、$H$、$I$,假设 $P\{F\} = 0.4$、$P\{G\} = 0.3$、$P\{H\} = 0.10$、$P\{I\} = 0.07$。将上面讨论的点绘制在同一坐标系内,并用折线将它们连接起来,得到下面的图:

图4 - 概率分析

这就我们想要的函数!并且这个函数的定义域应该是全体实数 $\mathbb{R}$ (因为样本的特征为实数),并且值域应该是 $(0,1)$(因为计算的是概率)。 那么有没有一个函数的曲线像上图一样,并且定义域为 $\mathbb{R}$ 值域为 $(0,1)$ 呢?有!这就是 $\sigma$ 函数:


$$ \sigma(z) = \dfrac{1}{1+e^{-z}} \tag{3} $$

其中 $z \in \mathbb{R}$ 。它的图形如下图左侧所示(为了比较将分析出的函数图像绘制在右侧):

图5 - sigmoid函数

我们先分析下这个函数的性质:当 $z = 0$ 时,函数的值为 $0.5$。如果 $z$ 非常大,即 $z \to +\infty$,就有 $e^{-z} \to 0$,进而 $\sigma(z) \to 1$。反之,如果 $z$ 非常小,即 $ z \to -\infty $,就有 $e^{-z} \to +\infty$,进而 $\sigma(z) \to 0$。

注意到 $\sigma$ 的参数是 $z$,而我们模型的输入是 $\mathbf{x}$ ,参数是 $\boldsymbol{\theta}$, 那么如何将它们关联起来呢?再来分析图3

图3 - 概率分析

图中的红色直线是分割线,它的方程为:$\boldsymbol{\theta}^T \mathbf{x} = 0$,在这条直线上的点为正类的概率都是 $0.5$(如 $A$ 点)。如果 $\boldsymbol{\theta}^T \mathbf{x} \to +\infty$,那么样本点就朝着正类方向(右上角)无限远离分割线,因此它是正类的概率就无限接近于 $1$(如 $E$ 点)。反之,如果 $ \boldsymbol{\theta}^T \mathbf{x} \to -\infty $,那么样本就朝着负类方向(左下角)无限远离分割线,因此它是正类的概率就无限接近于 $0$(如 $I$ 点)。

现在可以找到 $\mathbf{x}$ 和 $\boldsymbol{\theta}$ 与 $z$ 的对应关系了吧?我们将上面的分析汇总到表格里:

$z$ $\boldsymbol{\theta}^T \mathbf{x}$ 概率值
$0$ $0$ $0.5$
$+\infty$ $+\infty$ $1$
$-\infty$ $-\infty$ $0$

即 $\boldsymbol{\theta}^T \mathbf{x}$ 与 $z$ 成正比例关系!于是有 $z = \alpha \boldsymbol{\theta}^T \mathbf{x}$,其中 $\alpha$ 是任意正的常数。为了简单起见,令 $\alpha = 1$,则有:$ z = \boldsymbol{\theta}^T \mathbf{x} $

终于,到我们找到了想要的函数!我们的假设函数就变成:


$$ h_\theta(\mathbf{x}) = p(\mathbf{x}, \boldsymbol{\theta}) = \sigma(z) = \dfrac{1}{1+e^{-z}} = \dfrac{1}{1+e^{-\boldsymbol{\theta}^T \mathbf{x}}} \tag{4} $$

接着我们讨论第二个问题:知道了 $p$ 以后,如何计算模型参数 $\boldsymbol{\theta}$ 呢?

最大似然法

注意 $(2)$ 式 中 $y$ 的取值只能为 1 或 0,并且取 1 的概率为 $p$,取 0 的概率为 $1-p$,这正好符合 $(0-1)$ 分布。$(0-1)$ 分布的定义如下:

设随机变量 $X$ 只可能取 0 与 1 两个值,它的分布律是


$$ P\{X=k\} = p^k (1-p)^{1-k}, k=0,1\quad (0 < p < 1) \tag{5} $$

,则称 $X$ 服从以 $p$ 为参数的 $(0-1)$ 分布或两点分布。

$(2)$ 式代入 $(5)$ 式得:


$$ P\{Y=y\} = p^y (1-p)^{1-y} \tag{6} $$

现在我们知道了每个样本都服从 $(0-1)$ 分布,并且分布律为 $ p^y (1-p)^{1-y} $ 。在线性回归一篇中我们知道服从某个概率分布的样本可以使用最大似然法求解(如果还不了解最大似然法,请参见线性回归的相关内容)。最大似然法的核心思想是:假设真实值是由服从某个特定概率分布的随机变量产生的,算法的目标则是调整模型参数使模型产生这些真实值的概率达到最大。

首先列出似然函数,也就是所有样本概率的乘积:


$$ \mathcal{L} = \prod\limits_{i=1}^m \left[p^{(i)}\right]^{y^{(i)}} \left[\left(1-p^{(i)}\right)\right]^{1-y^{(i)}} \tag{7} $$

其中 $p^{(i)}$ 表示第 $i$ 个样本的概率,$y^{(i)}$ 表示第 $i$ 个样本的分类。为了方便计算,对上式两边取对数,得到对数似然函数:


$$ \log \mathcal{L} = \sum\limits_{i=1}^m \left[ y^{(i)} \log p^{(i)} + \left( 1 - y^{(i)} \right) \log \left( 1 - p^{(i)} \right) \right] \tag{8} $$

为了求上式的最大值,对其求导,过程如下:


\begin{align*} \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} &= \sum\limits_{i=1}^m y^{(i)} \dfrac{1}{p^{(i)}} \dfrac{\partial p^{(i)}}{\partial \boldsymbol{\theta}} - \left( 1-y^{(i)} \right) \dfrac{1}{1-p^{(i)}} \dfrac{\partial p^{(i)}}{\partial \boldsymbol{\theta}} \tag{9} \end{align*}

$(4)$ 式代入并令 $z = - \boldsymbol{\theta}^T \mathbf{x} $, 得:


\begin{align*} \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} =& \sum\limits_{i=1}^m y^{(i)} \dfrac{1}{p^{(i)}} \dfrac{\partial \sigma^{(i)}}{\partial z^{(i)}} \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} - \left( 1-y^{(i)} \right) \dfrac{1}{1-p^{(i)}} \dfrac{\partial \sigma^{(i)}}{\partial z^{(i)}} \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} \tag{10} \end{align*}

下面的求导过程用到了分式求导公式,这里复习一下:


$$ \left[ \dfrac{u(x)}{v(x)} \right]’ = \dfrac{u’(x)v(x) - u(x)v’(x)}{v^2(x)} \tag{11} $$

下面计算 $\dfrac{\partial \sigma}{\partial z}$:


\begin{align*} \dfrac{\partial \sigma}{\partial z} &= \left(\dfrac{1}{1+e^{-z}}\right)’ \\ &= \dfrac{e^{-z}}{(1+e^{-z})^2} \\ &= \dfrac{1+e^{-z}-1}{(1+e^{-z})^2} \\ &= \dfrac{1}{1+e^{-z}} - \dfrac{1}{(1+e^{-z})^2} \\ &= \dfrac{1}{1+e^{-z}}\left(1-\dfrac{1}{1+e^{-z}}\right) \\ &= \sigma(1-\sigma) = p(1-p) \tag{12} \end{align*}

$(12)$ 式代入 $(10)$ 式,继续整理:


\begin{align*} & \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} \\ &= \sum\limits_{i=1}^m y^{(i)} \dfrac{1}{p^{(i)}} p^{(i)} \left(1 - p^{(i)}\right) \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} - \left( 1 - y^{(i)} \right) \dfrac{1}{1-p^{(i)}} p^{(i)} \left(1 - p^{(i)}\right) \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} \\ &= \sum\limits_{i=1}^m y^{(i)} \left( 1 - p^{(i)} \right) \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} - \left( 1 - y^{(i)} \right) p^{(i)} \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} \\ &= \sum\limits_{i=1}^m \left( y^{(i)} - y^{(i)} p^{(i)} - p^{(i)} + y^{(i)} p^{(i)} \right) \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} \\ &= \sum\limits_{i=1}^m \left( y^{(i)} - p^{(i)} \right) \dfrac{\partial z^{(i)}}{\partial \boldsymbol{\theta}} \\ &= \sum\limits_{i=1}^m \left( y^{(i)} - p^{(i)} \right) \dfrac{\partial \boldsymbol{\theta}^T \mathbf{x}^{(i)}}{\partial \boldsymbol{\theta}} & ① \\ &= \sum\limits_{i=1}^m \left( y^{(i)} - p^{(i)} \right) \mathbf{x}^{(i)} \quad & ② \\ &= \mathbf{X}^T \left( \mathbf{y} - \mathbf{p} \right) & ③ \tag{13} \end{align*}

$① \Rightarrow ②$ 是向量与向量的乘积对向量的导数,推导过程为(注意下面的推导省略了讨厌的角标,即令 $ \mathbf{x} = \mathbf{x}^{(i)}$):


\begin{align*} \dfrac{\partial \boldsymbol{\theta}^T \mathbf{x}}{\partial \boldsymbol{\theta}} &= \dfrac{\partial}{\partial \boldsymbol{\theta}} \begin{bmatrix} \theta_0 & \theta_1 & \cdots & \theta_n \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_n \end{bmatrix} \\ &= \dfrac{\partial}{\partial \boldsymbol{\theta}} \sum\limits_{i=0}^n \theta_i x_i \\ &= \begin{bmatrix} \dfrac{\partial \sum\limits_{i=0}^n \theta_i x_i}{\partial \theta_0} \\ \dfrac{\partial \sum\limits_{i=0}^n \theta_i x_i}{\partial \theta_1} \\ \vdots \\ \dfrac{\partial \sum\limits_{i=0}^n \theta_i x_i}{\partial \theta_n} \end{bmatrix} \\ &= \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_n \end{bmatrix} \\ &= \mathbf{x} \\ &= \mathbf{x}^{(i)} \tag{14} \end{align*}

$ ② \Rightarrow ③ $ 是矩阵与向量乘积的逆推,令 $ a_i = y^{(i)} - p^{(i)}$,$ \mathbf{x}_i = \mathbf{x}^{(i)} = \begin{bmatrix} x_{i0} & x_{i1} & \cdots & x_{in} \end{bmatrix}^T $,推导过程为:


\begin{align*} \sum\limits_{i=1}^m \left( y^{(i)} - p^{(i)} \right) \mathbf{x}^{(i)} &= \sum\limits_{i=1}^m a_i \mathbf{x}_{i} \\ &= a_1 \mathbf{x}_1 + a_2 \mathbf{x}_2 + \cdots + a_m \mathbf{x}_m \\ &= \begin{bmatrix} a_1 x_{10} \\ a_1 x_{11} \\ \vdots \\ a_1 x_{1n} \end{bmatrix} + \begin{bmatrix} a_2 x_{20} \\ a_2 x_{21} \\ \vdots \\ a_2 x_{2n} \end{bmatrix} + \cdots + \begin{bmatrix} a_m x_{m0} \\ a_m x_{m1} \\ \vdots \\ a_m a_{mn} \end{bmatrix} \\ &= \begin{bmatrix} a_1 x_{10} + a_2 x_{20} + \cdots + a_m x_{m0} \\ a_1 x_{11} + a_2 x_{21} + \cdots + a_m x_{m1} \\ \vdots \\ a_1 x_{1n} + a_2 x_{2n} + \cdots + a_m x_{mn} \end{bmatrix} \\ &= \begin{bmatrix} x_{10} & x_{20} & \cdots & x_{m0} \\ x_{11} & x_{21} & \cdots & x_{m1} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1n} & x_{2n} & \cdots & x_{mn} \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_m \end{bmatrix} \\ &= \mathbf{X}^T \mathbf{a} \\ &= \mathbf{X}^T (\mathbf{y} - \mathbf{p}) \tag{15} \end{align*}

$$ \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} = \mathbf{X}^T \left( \mathbf{y} - \mathbf{p} \right) \tag{13} $$

下面我们令 $(13)$ 式为 $\mathbf{0}$ 来求解 $\boldsymbol{\theta}$:


\begin{align*} \begin{array}{rl} \text{令:} & \mathbf{X}^T \left( \mathbf{y} - \mathbf{p} \right) = \mathbf{0} \\ \text{整理,得:} & \mathbf{p} = \mathbf{y} \\ \text{将向量展开,得:} \\ & \begin{bmatrix} \dfrac{1}{1+e^{-\boldsymbol{\theta}^T\mathbf{x}^{(1)}}} \\ \dfrac{1}{1+e^{-\boldsymbol{\theta}^T\mathbf{x}^{(2)}}} \\ \vdots \\ \dfrac{1}{1+e^{-\boldsymbol{\theta}^T\mathbf{x}^{(m)}}} \end{bmatrix} = \begin{bmatrix} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)} \end{bmatrix} \\ \text{整理成关于} \theta_i \text{的方程,得:} \\ & \left\{ \begin{array}{l} x^{(1)}_0 \theta_0 + x^{(1)}_1 \theta_1 + \cdots + x^{(1)}_n \theta_n = \log\dfrac{y^{(1)}}{1-y^{(1)}} \\ x^{(2)}_0 \theta_0 + x^{(2)}_1 \theta_1+ \cdots + x^{(2)}_n \theta_n = \log\dfrac{y^{(2)}}{1-y^{(2)}} \\ \cdots\cdots\cdots\cdots \\ x^{(m)}_0 \theta_0 + x^{(m)}_1 \theta_1 + \cdots + x^{(m)}_n \theta_n = \log\dfrac{y^{(m)}}{1-y^{(m)}} \end{array} \right. \tag{16} \end{array} \end{align*}

将 $\log\dfrac{y^{(i)}}{1-y^{(i)}}$ 记为 $b^{(i)}$,则上面的方程组可简记为:


$$ \mathbf{X} \boldsymbol{\theta} = \mathbf{b} \tag{17} $$

这个方程组有没有解呢?我们先来复习一下线性代数的知识($R(\mathbf{A})$ 表示求矩阵 $\mathbf{A}$ 的秩):

定理 $n$ 元线性方程组 $\mathbf{A}\mathbf{x} = \mathbf{b}$

(ⅰ) 无解的充分必要条件是 $R(\mathbf{A}) \lt R(\mathbf{A}, \mathbf{b}) $
(ⅱ) 有惟一解的充分必要条件是 $R(\mathbf{A}) = R(\mathbf{A}, \mathbf{b}) = n$
(ⅲ) 有无限多解的充分必要条件是 $R(\mathbf{A}) = R(\mathbf{A}, \mathbf{b}) < n$

下面分别求 $(17)$ 式的 $R(\mathbf{X})$ 和 $R(\mathbf{X}, \mathbf{b})$。

假设样本数量为 $m$,样本特征为 $n+1$ (包括 bias),那么 $\mathbf{X}$ 就是一个 $m \times (n+1)$ 的矩阵,而在实际中 $m$ 是远大于 $n+1$ 的,假设没有任何两个样本是完全一样的,那么 $\mathbf{X}$ 的秩就应该为 $n+1$,即 $R(\mathbf{X}) = n + 1$。再来看 $R(\mathbf{X}, \mathbf{b})$,因为 $b^{(i)} = \log\dfrac{y^{(i)}}{1-y^{(i)}}$ 与 $x^{(i)}_{j}$ (样本 $\mathbf{x}^{(i)}$ 的第 $j$ 个特征)并不是线性关系(因为在运算中有个 $\log$),当然也就无法将矩阵 $(\mathbf{X},\mathbf{b})$ 经过有限次的列变换(注意列变换是线性变换)变换成 $(\mathbf{X}, \mathbf{0})$,因此矩阵 $(\mathbf{X}, \mathbf{b})$ 的秩应该为 $n + 2$,即 $R(\mathbf{X}, \mathbf{b}) = n + 2$。

于是有 $ R(\mathbf{X}) \lt R(\mathbf{X}, \mathbf{b}) $,根据上面定理中 (ⅰ) 的结论,方程组 $(17)$ 无解。(请读者自行分析 $m \leqslant n + 1$ 的情况。)


$$ \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} = \mathbf{X}^T \left( \mathbf{y} - \mathbf{p} \right) \tag{13} $$

好了,我们无法令 $(13)$ 式为 $\mathbf{0}$ 来求解 $\boldsymbol{\theta}$,但是我们还有其他方法,比如梯度下降(上升)法,而梯度下降(上升)法在函数为凸(凹)函数时总是可以达到最优值。为了检验这个似然函数是不是凸(凹)函数,我们对它求二次导数。为了书写方便,下面的推导中,将元素为 $a_{ij}$ 的 $m \times n$ 矩阵简记为 $\left( a_{ij} \right)_{m \times n}$,根据 $(13)$ 式,有:


\begin{align*} \dfrac{\partial \log ^2 \mathcal{L}}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} &= \dfrac{\partial}{\partial \boldsymbol{\theta}^T} \mathbf{X}^T (\mathbf{y} - \mathbf{p}) \\ &= - \dfrac{\partial \mathbf{X}^T \mathbf{p}}{\partial \boldsymbol{\theta}^T} \\ &= - \dfrac{\partial}{\partial \boldsymbol{\theta}^T} \begin{bmatrix} x_{10} & x_{20} & \cdots & x_{m0} \\ x_{11} & x_{21} & \cdots & x_{m1} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1n} & x_{2n} & \cdots & x_{mn} \end{bmatrix} \begin{bmatrix} p_1 \\ p_2 \\ \vdots \\ p_m \end{bmatrix} \\ &= - \dfrac{\partial}{\partial \boldsymbol{\theta}^T} \begin{bmatrix} \sum\limits_{k=1}^m x_{k0} p_k \\ \sum\limits_{k=1}^m x_{k1} p_k \\ \vdots \\ \sum\limits_{k=1}^m x_{kn} p_k \end{bmatrix} \\ &= - \begin{bmatrix} \dfrac{\partial}{\partial \theta_0} \sum\limits_{k=1}^m x_{k0} p_k & \dfrac{\partial}{\partial \theta_1} \sum\limits_{k=1}^m x_{k0} p_k & \cdots & \dfrac{\partial}{\partial \theta_n} \sum\limits_{k=1}^m x_{k0} p_k \\ \dfrac{\partial}{\partial \theta_0} \sum\limits_{k=1}^m x_{k1} p_k & \dfrac{\partial}{\partial \theta_1} \sum\limits_{k=1}^m x_{k1} p_k & \cdots & \dfrac{\partial}{\partial \theta_n} \sum\limits_{k=1}^m x_{k1} p_k \\ \vdots & \vdots & \ddots & \vdots \\ \dfrac{\partial}{\partial \theta_0} \sum\limits_{k=1}^m x_{kn} p_k & \dfrac{\partial}{\partial \theta_1} \sum\limits_{k=1}^m x_{kn} p_k & \cdots & \dfrac{\partial}{\partial \theta_n} \sum\limits_{k=1}^m x_{kn} p_k \end{bmatrix} \\ &= - \left( \dfrac{\partial}{\partial \theta_j} \sum\limits_{k=1}^m x_{ki} p_k \right) _{n+1 \times n+1} \\ &= - \left( \sum\limits_{k=1}^m x_{ki} \dfrac{\partial p_k}{\partial \theta_j} \right) _{n+1 \times n+1} \\ &= - \left( \sum\limits_{k=1}^m x_{ki} \dfrac{\partial \sigma_k}{\partial z_k} \dfrac{\partial z_k}{\partial \theta_j} \right) _{n+1 \times n+1} \\ &= - \left( \sum\limits_{k=1}^m x_{ki} p_k (1 - p_k) \dfrac{\partial}{\partial \theta_j} \sum\limits_{r=0}^n \theta_r x_{kr} \right) _{n+1 \times n+1} \\ &= - \left( \sum\limits_{k=1}^m x_{ki} p_k (1 - p_k) x_{kj} \right) _{n+1 \times n+1} \\ &= - \left( \begin{bmatrix} x_{1i} p_1 (1 - p_1) & x_{2i} p_2 (1 - p_2) & \cdots & x_{mi} p_m (1-p_m) \end{bmatrix} \begin{bmatrix} x_{1j} \\ x_{2j} \\ \vdots \\ x_{mj} \end{bmatrix} \right) _{n+1 \times n+1} \\ &= - \begin{bmatrix} x_{10} p_1 (1 - p_1) & x_{20} p_2 (1 - p_2) & \cdots & x_{m0} p_m (1-p_m) \\ x_{11} p_1 (1 - p_1) & x_{21} p_2 (1 - p_2) & \cdots & x_{m1} p_m (1-p_m) \\ \vdots & \vdots & \ddots & \vdots \\ x_{1n} p_1 (1 - p_1) & x_{2n} p_2 (1 - p_2) & \cdots & x_{mn} p_m (1-p_m) \end{bmatrix} \begin{bmatrix} x_{10} & x_{11} & \cdots & x_{1n} \\ x_{20} & x_{21} & \cdots & x_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ x_{m0} & x_{m1} & \cdots & x_{mn} \end{bmatrix} \\ &= - \begin{bmatrix} x_{10} & x_{20} & \cdots & x_{m0} \\ x_{11} & x_{21} & \cdots & x_{m1} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1n} & x_{2n} & \cdots & x_{mn} \end{bmatrix} \begin{bmatrix} p_1 (1 - p_1) & 0 & \cdots & 0 \\ 0 & p_2 (1 - p_2) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & p_m (1-p_m) \end{bmatrix} \mathbf{X} \\ &= - \mathbf{X}^T \begin{bmatrix} p_1 & 0 & \cdots & 0 \\ 0 & p_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & p_m \end{bmatrix} \begin{bmatrix} 1 - p_1 & 0 & \cdots & 0 \\ 0 & 1 - p_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1-p_m \end{bmatrix} \mathbf{X} \\ &= - \mathbf{X}^T \textbf{diag}\left(\mathbf{p}\right) \textbf{diag}\left(\mathbf{1-p}\right) \mathbf{X} \tag{18} \end{align*}

上式即为对数似然函数 $\log \mathcal{L}$ 在任意一点 $\mathbf{X}$ 的海森矩阵,在线性回归中我们已经证明了 $\mathbf{X}^T \mathbf{X}$ 是正定矩阵,而 $\textbf{diag}\left(\mathbf{p}\right) \textbf{diag}\left(\mathbf{1-p}\right)$ 在矩阵相乘时可以看做正的常数,又由于上式前面的负号,因此上式为负定矩阵。因此得出结论:对数似然函数 $\log \mathcal{L}$ 在任意一点 $\mathbf{X}$ 的海森矩阵都为负定矩阵,也就是说对数似然函数 $\log \mathcal{L}$ 为凹函数(注意这里的凸函数和凹函数的定义以基维百科为准)。


$$ \dfrac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}} = \mathbf{X}^T \left( \mathbf{y} - \mathbf{p} \right) \tag{13} $$

到此我们证明了对数似然函数可以使用梯度上升法求得最大值。设学习率为 $\eta$,并以 $(13)$ 式为梯度,则梯度上升可以表示为($:=$ 表示赋值):


$$ \boldsymbol{\theta} := \boldsymbol{\theta} + \eta \mathbf{X}^T \left( \mathbf{y - p} \right) \tag{19} $$

示例

下面我们以鸢尾花数据集训练一个二分类分类器,该样本集包含了三类鸢尾花:Setosa、Versicolour、Virginica。这里仅用 Versicolour 与 Virginica 样本集来训练,判断一个样本是不是 Virginica,数据集中每个样本都有四个特征:花萼长、花萼宽、花瓣长、花瓣宽,为了简单起见,我们仅用花瓣长和花瓣宽作为特征。首先加载数据并按照花瓣长和花瓣宽来绘制样本分布图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from sklearn.datasets import load_iris

iris = load_iris()
print(iris.DESCR)

def plot_iris(iris_data):
plt.figure(figsize=(15, 8))
for index, style in zip((1,2), ('bs', 'g^')):
data = iris_data[iris.target == index]
# 第 2、3 列分别表示花瓣长、宽
plt.plot(data[:,2], data[:,3], style, label=iris.target_names[index])
plt.xticks(fontsize=18); plt.xlabel("petal length", fontsize=18)
plt.yticks(fontsize=18); plt.ylabel("petal width", fontsize=18)
plt.legend(loc='best', fontsize=18)

plot_iris(iris.data)
plt.show()

如下图:

图6 - 鸢尾花样本集

接下来就是套用 $(19)$ 式来进行梯度上升算法了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def sigmoid(z):
return 1 / (1 + np.exp(-z))

x = iris.data[:,(2,3)]
y = (iris.target == 2).astype(np.int).reshape(-1, 1)
m = x.shape[0]
X = np.c_[np.ones([m, 1]), x]
theta = np.random.randn(3, 1)

epoches = 10000
eta = 0.01

for epoch in range(epoches):
theta += eta * X.T.dot(y - sigmoid(X.dot(theta)))

将训练结果绘制出来,如下:

1
2
3
4
5
6
left_right = np.array([2.9, 7])
boundary = -(theta[1] * left_right + theta[0]) / theta[2]

plot_iris(iris.data)
plt.plot(left_right, boundary, 'r-')
plt.show()

结果为:

图7 - 训练结果

sklearn 实现

也可以使用 sklearn 的 LogisticRegression 来达到同样的效果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.linear_model import LogisticRegression

x = iris.data[:,(2,3)]
y = (iris.target == 2).astype(np.int)

log_reg = LogisticRegression(C=10**10)

log_reg.fit(x, y)

left_right = np.array([2.9, 7])
boundary = - (log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]

plot_iris(iris.data)
plt.plot(left_right, boundary, 'r-')
plt.show()

图像为:

图8 - sklearn训练结果

模型的评估

请参见模型评估指标


逻辑回归就讲到这里,如果有任何疑问或建议,请在留言区留言!

版权声明:本文为原创文章,转载请注明出处。http://cynhard.com/2018/07/08/ML-Logistic-Regression/

推荐文章