逻辑回归模型及lbfgs算法公式推导

本文最后更新于：2022年10月25日凌晨

机器学习第一次作业记录。小组成员：方桂安，刘玥，周敏。

一、描述逻辑回归模型

1.1数据

\[ D=\{(x_1,y_1),(x_2,y_2),\dots,(x_N,y_N)\},x_i\in \mathbb{R}^n,y_i\in\{0,1\} \]

1.2模型

最初模型：

\[ f(x_i)=\omega^Tx_i+b,使得f(x_i)\simeq g(y_i) \]

我们的标记变量y的范围是0或1，所以我们需要一个函数能够将上述x的线性组合转化为0或1，最理想的是阶跃函数。

\[ 阶跃函数：y=g^{-1}(\omega^Tx+b)= \begin{cases} 0, & \omega^Tx+b<0\\ 0.5, & \omega^Tx+b=0\\ 1, & \omega^Tx+b>0 \end{cases} \]

但由于阶跃函数不连续，不满足单调可微的条件。所以我们希望通过一个一定程度上近似阶跃函数的“替代函数”，并且希望它单调可微。由此，我们想到了逻辑斯蒂函数。

它的图像如下：

因为Logistic 回归主要用于分类问题，以二分类为例，对于所给数据集假设存在这样的一条直线可以将数据完成线性可分。

当我们要找到分类概率 P(Y=1) 与输入向量 x 的直接关系时，我们引入Sigmoid函数，然后通过比较概率值来判断类别。

引入sigmoid函数具体实现如下：

但因为逻辑斯蒂函数的值域在[0,1]之间，无法直接输出0或1。在此基础上，考虑到$\omega^Tx+b$取值是连续的，因此它不能拟合离散变量。可以考虑用它来拟合条件概率$p(y=1|x)$，因为概率的取值也是连续的,我们将逻辑斯蒂函数的输出作为输入x能预测到y为1的概率，并利用对数几率函数，得到下面三个式子。通过此方法，我们将线性模型转换为概率模型。

1.3策略

在策略上，我们采用极大似然法。即选择最优的w，b使得我们输入x得到的正确的y的概率最大，即下式：

\[ (w^*,b^*)=\mathop{argmax}\limits_{(w,b)}\prod_{i=1}^Np(y_i|x_i;\omega,b) \]

我们这里做一点变换：

因为上式是连乘的函数，我们通过对数似然函数将之转化为求和，即下式：

为了方便计算，我们做以下处理

\[ assume \ that\ \hat{\omega}=(\omega;b),\hat{x}=(x;1) \]

则上式可化为

\[ \hat{\omega^*}=\mathop{argmin}\limits_{\hat{\omega}}\sum_{i=1}^N(-y_i\hat{\omega }x_i+ln(1+e^{\hat{\omega}^T\hat{x}_i})) \]

这是一个凸函数，可用经典的数值优化算法，如梯度下降法、牛顿法求解。

最终，我们学得的逻辑斯蒂回归模型为

二、描述训练模型所使用的算法

2.1梯度下降法

2.1.1问题分析

首先，我们的目标是下式

\[ E(\hat{\omega})=\sum_{i=1}^N(-y_i\hat{\omega}^T\hat{x}_i+ln(1+e^{\hat{\omega}^T\hat{x}_i})),\hat{\omega}^*=\mathop{argmin}_{\hat{\omega}}E(\hat{\omega}) \]

梯度下降法是一种迭代算法：我们选取适当的初始值$\hat{\omega}^{(0)}$，不断迭代，更新$\hat{\omega}$的值，进行目标函数的极小化，直到收敛。由于负梯度方向是使得函数值下降最快的方向，在迭代的每一步，以负梯度方向更新$\hat{\omega}$的值，从而达到减小函数值的目的。如下图形象化表示：

2.1.2核心思想：

$E(\hat{\omega})$是具有一阶连续偏导数的凸函数，其极值点在一阶导数为零的地方取得
一阶泰勒展开：$E(\hat{\omega})\thickapprox E(\hat{\omega}^{(k)})+\nabla E(\hat{\omega}^{(k)})(\hat{\omega}-\hat{\omega}^{(k)})$，其中，$\nabla E(\hat{\omega}^{(k)})$是$E(\hat{\omega})$在$\hat{\omega}^{(k)}$的梯度：

\[ \nabla E(\hat{\omega}^{(k)})=\frac{\partial E(\hat{\omega})}{\partial \hat{\omega}}|_{\hat{\omega}=\hat{\omega}^{(k)}} \]

求取第k+1次迭代值：$\hat{\omega}^{k+1}=\hat{\omega}^{(k)}+\eta_k*(-\nabla E(\hat{\omega}^{(k)}))$，其中$\eta_k$是步长，由我们最初指定。梯度推导：

\[ E(\hat{\omega})=\sum_{i=1}^N(-y_i\hat{\omega}^T\hat{x}_i+ln(1+e^{\hat{\omega}^T\hat{x}_i}))\\ \nabla E(\hat{\omega}^{(k)})=\sum_{i=1}^N-y_i\hat{x}_i+\frac{1}{1+e^{\hat{\omega}^T\hat{x}_i}}*e^{\hat{\omega}^T\hat{x}_i}*\hat{x}_i\\ =-\sum_{i=1}^Nx_i(y_i-\frac{e^{\hat{\omega}^T\hat{x}_i}}{1+e^{\hat{\omega}^T\hat{x}_i}}) \]

2.1.3伪代码：

输入：目标函数$E(\hat{\omega})$，梯度函数$\nabla E(\hat{\omega})$，计算精度ε，步长$\eta_k$；

输出： $E(\hat{\omega})$的极小点$\hat{\omega}^*$。

（1）取初始值$\hat{\omega}^{(0)}\in \mathbb{R}^{d+1}$，置k=0；

（2）计算$E(\hat{\omega}^{(k)})$；

（3）计算梯度$\nabla E(\hat{\omega}^{(k)})$，当$||\nabla E(\hat{\omega}^{(k)})||<\varepsilon$时，令$\hat{\omega}^*=\hat{\omega}^{(k)}$，

停止迭代；

（4）置$\hat{\omega}^{(k+1)}=\hat{\omega}^{(k)}+\eta_k(-\nabla E(\hat{\omega}^{(k)}))$，计算$E(\hat{\omega}^{(k+1)})$，

当$||E(\hat{\omega}^{(k+1)})-E(\hat{\omega}^{(k)})||<\varepsilon$或$||\hat{\omega}^{(k+1)}-\hat{\omega}^{(k)}||<\varepsilon$时，

令$\hat{\omega}^*=\hat{\omega}^{(k)}$，停止迭代；

（5）否则，置k=k+1，转步骤（3）。

2.1.4分析

优点：方法简单，易理解

缺点：迭代次数多，下降速度慢，如下图，我们采用梯度下降法，迭代近50000次才收敛

且准确率如下，可以看出准确率不高。

2.2牛顿法

2.2.1核心思想：

$E(\hat{\omega})$是具有二阶连续偏导数的函数

二阶泰勒展开：$E(\hat{\omega})\thickapprox E(\hat{\omega}^{(k)})+\nabla E(\hat{\omega}^{(k)})(\hat{\omega}-\hat{\omega}^{(k)})+\frac{1}{2}(\hat{\omega}-\hat{\omega}^{(k)})^TH(\hat{\omega}^{(k)})(\hat{\omega}-\hat{\omega}^{(k)})$

\[ \nabla E(\hat{\omega})=\frac{\partial E(\hat{\omega})}{\partial \hat{\omega}}|_{(d+1)\times 1},H(\hat{\omega})=\frac{\partial ^2E(\hat{\omega})}{\partial \hat{\omega}_i \partial \hat{\omega}_j}|_{(d+1)\times 1} \]

利用二阶泰勒展开$E(\hat{\omega})$取极小点的必要条件$\nabla E(\hat{\omega})=0$，在第k次迭代$\hat{\omega}^{(k)}$，求$\nabla E(\hat{\omega}^{(k)})+H(\hat{\omega}^{(k)})(\hat{\omega}-\hat{\omega}^{(k)})=0$的点，作为第k+1次迭代值$\hat{\omega}^{(k+1)}$

2.2.2伪代码

输入：目标函数$E(\hat{\omega})$，梯度函数$\nabla E(\hat{\omega})$，海森矩阵$H(\hat{\omega})$，精度ε；

输出：$E(\hat{\omega})$的极小点$\hat{\omega}^*$。

（1）取初始值$\hat{\omega}^{(0)}\in \mathbb{R}^{n+1}$，置k=0；

（2）计算梯度$\nabla E(\hat{\omega}^{(k)})$；

（3）当$||E(\hat{\omega}^{(k)})||<\varepsilon$时，令$\hat{\omega}^*=\hat{\omega}^{(k)}$，停止迭代；

否则，计算海森矩阵$H(\hat{\omega}^{(k)})$ ；

（4）置$\hat{\omega}^{(k+1)}=\hat{\omega}^{(k)}-(H(\hat{\omega}))^{(-1)}\nabla E(\hat{\omega}^{(k)})$；

（5）置k=k+1，转步骤（2）。

2.2.3分析

牛顿法优点：下降速度快，属于二次收敛

缺点：海森矩阵计算复杂度高，且要求可逆才能计算，所以我们查阅资料，将采用拟牛顿法。

2.3 BFGS算法:

由于上述牛顿公式中可以看出，我们的海森矩阵不易得到，因此我们有以下迭代公式来逼近海森矩阵：

\[ H_{k+1}=H_k+\frac{y_ky_k^T}{y_k^Ts_k}-\frac{H_ks_ks_k^TH_k^T}{s_k^TH_k^Ts_k} \]

但计算量还是很大，矩阵相乘太多。所以我们最终采取$Sherman-Morrison$公式进行变换可得：

\[ H_{k+1}=\left(I-\frac{s_{k} y_{k}^{T}}{y_{k}^{T} s_{k}}\right) H_{k}\left(I-\frac{y_{k} s_{k}^{T}}{y_{k}^{T} s_{k}}\right)+\frac{s_{k} s_{k}^{T}}{y_{k}^{T} s_{k}} \quad(1) \]

公式推导如下：

\[ \begin{array}{l} \text { Sherman Morrison 公式: }\\ \left(\mathrm{A}+\frac{u u^{T}}{t}\right)^{-1}=A^{-1}-\frac{A^{-1} u u^{T} A^{-1}}{t+u^{T} A^{-1} u}\\ \left(\mathrm{H}+\frac{y y^{T}}{y^{T} \mathrm{~s}}-\frac{H s s^{T} \mathrm{H}}{s^{T} H s}\right)^{-1}\\ =\left(\mathrm{H}+\frac{y y^{T}}{y^{T} \mathrm{~s}}\right)^{-1}+\left(\mathrm{H}+\frac{y y^{T}}{y^{T} \mathrm{~s}}\right)^{-1} \frac{H s s^{T} H}{s^{T} H^{T} s-s^{T} H\left(\mathrm{H}+\frac{y y^{T}}{y^{T} \mathrm{~S}}\right)^{-1} \mathrm{Hs}}\left(\mathrm{H}+\frac{y y^{T}}{y^{T} \mathrm{~S}}\right)^{-1}\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} S+y^{T} H^{-1} y}\right)+\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} S+y^{T} H^{-1} y}\right) \frac{H s s^{T} H}{s^{T} H S-s^{T} H\left(H^{-1}-\frac{H^{-1} y y^{T}-1}{y^{T} s+y^{T} H^{-1} y}\right) H s}\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} S+y^{T} H^{-1} y}\right)\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right)+\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right) \frac{H s s^{T} H}{\frac{s^{T} y y^{T} s}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}}\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~S}+y^{T} H^{-1} y}\right)\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right)+\frac{H^{-1} H s s^{T} H H^{-1}}{\frac{s^{T} y y^{T} S}{y^{T} s+y^{T} H^{-1} y}}-\frac{H^{-1} H s s^{T} H}{\frac{s^{T} y y^{T} s}{}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y} H^{-1} \frac{y y^{T}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y} H^{-1}\\ -\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y} \frac{H s s^{T} H}{\frac{s^{T} y y^{T} S}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}} H^{-1}\\ +H^{-1} \frac{y y^{T}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y} H^{-1} \frac{H s s^{T} H}{\frac{s^{T} y y^{T} S}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}} H^{-1} \frac{y y^{T}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y} H^{-1}\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right)+\frac{s s^{T}\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right)}{s^{T} y y^{T} s}-\frac{s s^{T} y y^{T} H^{-1}}{s^{T} y y^{T} S}-\frac{H^{-1} y y^{T} S S^{T}}{s^{T} y y^{T} S}\\ +\frac{H^{-1} y y^{T} S s^{T} y y^{T} H^{-1}}{\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right) s^{T} y y^{T} S}\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right)+\frac{s s^{T}\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right)}{\left(s^{T} y\right)^{2}}-\frac{s\left(s^{T} y\right) y^{T} H^{-1}}{\left(s^{T} y\right)^{2}}-\frac{H^{-1} y\left(y^{T} s\right) s^{T}}{\left(s^{T} y\right)^{2}}\\ +\frac{H^{-1} y\left(y^{T} S s^{T} y\right) y^{T} H^{-1}}{\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right) s^{T} y y^{T} S}\\ =\left(H^{-1}-\frac{H^{-1} y y^{T} H^{-1}}{y^{T} \mathrm{~s}+y^{T} H^{-1} y}\right)+\frac{s s^{T}\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right)}{\left(s^{T} y\right)^{2}}-\frac{s y^{T} H^{-1}}{s^{T} y}-\frac{H^{-1} y s^{T}}{s^{T} y}+\frac{H^{-1} y y^{T} H^{-1}}{\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right)}\\ =H^{-1}+\frac{s s^{T}\left(y^{T} \mathrm{~s}+y^{T} H^{-1} y\right)}{\left(s^{T} y\right)^{2}}-\frac{s y^{T} H^{-1}}{s^{T} y}-\frac{H^{-1} y s^{T}}{s^{T} y}\\ =H^{-1}+\frac{s s^{T} y^{T} \mathrm{~s}}{\left(s^{T} y\right)^{2}}+\frac{s s^{T} y^{T} H^{-1} y}{\left(s^{T} y\right)^{2}}-\frac{s y^{T} H^{-1}}{s^{T} y}-\frac{H^{-1} y s^{T}}{s^{T} y}\\ =H^{-1}\left(I-\frac{y s^{T}}{s^{T} y}\right)-\frac{s y^{T} H^{-1}}{s^{T} y}\left(I-\frac{y s^{T}}{s^{T} y}\right)+\frac{s s^{T}}{s^{T} y}\\ =\left(I-\frac{s y^{T}}{s^{T} y}\right) H^{-1}\left(I-\frac{y s^{T}}{s^{T} y}\right)+\frac{s s^{T}}{s^{T} y} \end{array} \]

2.4 L-BFGS算法

$$ \[\begin{array}{c} H_{k+1}=V_{k}^{T} H_{k} V_{k}+\rho_{k} s_{k} s_{k}^{T} \\ \end{array}\]

给定初始矩阵$H_0=I$，利用上式，可得：

\[ \begin{aligned} H_{1}&=V_{0}^{T} H_{0} V_{0}+\rho_{0} s_{0} s_{0}^{T}\\ H_{2} &=V_{1}^{T} H_{1} V_{1}+\rho_{1} s_{1} s_{1}^{T} \\ &=V_{1}^{T}\left(V_{0}^{T} H_{0} V_{0}+\rho_{0} s_{0} s_{0}^{T}\right) V_{1}+\rho_{1} s_{1} s_{1}^{T} \\ &\left.=V_{1}^{T} V_{0}^{T} H_{0} V_{0} V_{1}+V_{1}^{T} \rho_{0} s_{0} s_{0}^{T}\right) V_{1}+\rho_{1} s_{1} s_{1}^{T} \\ & \\ \quad H_{k+1} &=\left(V_{k}^{T} V_{k-1}^{T} \ldots V_{1}^{T} V_{0}^{T}\right) H_{0}\left(V_{0} V_{1} \ldots V_{k-1} V_{k}\right) \\ &+\left(V_{k}^{T} V_{k-1}^{T} \ldots V_{1}^{T}\right) \rho_{1} s_{1} s_{1}^{T}\left(V_{1} \ldots V_{k-1} V_{k}\right) \\ &+\ldots \\ &+\left(V_{k}^{T}\right) \rho_{k-1} s_{k-1} s_{k-1}^{T}\left(V_{k}\right) \\ &+\rho_{k} s_{k} s_{k}^{T} \end{aligned} \]

只保留最近的m步后，上式的迭代公式变为：

\[ \begin{aligned} H_{k+1} &=\left(V_{k}^{T} V_{k-1}^{T} \ldots V_{k-m}^{T}\right) H_{0}\left(V_{k-m} \ldots V_{k-1} V_{k}\right) \\ &+\left(V_{k}^{T} V_{k-1}^{T} \ldots V_{k-m+1}^{T}\right) \rho_{k-m} s_{k-m} s_{k-m}^{T}\left(V_{k-m+1} \ldots V_{k-1} V_{k}\right) \\ +& \ldots \\ &+\left(V_{k}^{T}\right) \rho_{k-1} s_{k-1} s_{k-1}^{T}\left(V_{k}\right) \\ &+\rho_{k} s_{k} s_{k}^{T} \end{aligned} \]

所求方向为：

\[ \begin{aligned} H_{k} \nabla f &=\left(V_{K-1}^{T} V_{K-2}^{T} \ldots V_{K-m}^{T}\right) H_{0}\left(V_{K-m} V_{K-m+1} \ldots V_{K-1}\right) \nabla f \\ &+\left(V_{K-1}^{T} \ldots V_{K-m+1}^{T}\right)\rho_{k-m} s_{k-m} s_{k-m}^T(V_{k-m+1}\dots V_{k-1}V_{k}) \nabla f\\ &+\ldots \\ &+V_{k-1} \rho_{k-1} s_{k-1}s_{k-1}^TV_k\nabla f \\ &+\rho_{k} s_{k}s_{k}^T\nabla f \end{aligned} \]

Two-Loop 算法：

\[ \begin{array}{l} q_{k} \leftarrow \nabla f_{k} \\ \text { for } i=k-1 \text { to } k-m \text { do } \\ \quad \alpha_{i}=\rho_{i} s_{i}^{T} q_{i+1} \\ q_{i}=q_{i+1}-\alpha_{i} y_{i} \\ \text { end for } \\ r_{k-m-1}=H_{0} q_{k-m} \\ \text { for } i=k-m, k-m+1 \text { to } k-1 \text { do } \\ \quad \beta_{i}=\rho_{i} y_{i}^{T} r_{i-1} \\ r_{i}=r_{i-1}+s_{i} \alpha_{i}-\beta_{i} \\ \text { end for } \\ \text { End, The result is } H_{k+1} \nabla f=r \end{array} \]

Two-Loop算法解析---第一个循环：

重写公式：

\[ \begin{aligned} H_{k} \nabla f &=\left(V_{K-1}^{T} V_{K-2}^{T} \ldots V_{K-m}^{T}\right) H_{0}\left(V_{K-m} V_{K-m+1} \ldots V_{K-1}\right) \nabla f \\ &+\left(V_{K-1}^{T} \ldots V_{K-m+1}^{T}\right) s_{k-m} \alpha_{k-m} \\ &+\ldots \\ &+V_{k-1} s_{k-1} \alpha_{k-2} \\ &+s_{k-1} \alpha_{k-1} \end{aligned} \]

Two-Loop算法解析---第二个循环：

初始:

\[ r_{k-\mathrm{m}}=V_{k-\mathrm{m}} H_{0} V_{k-\mathrm{m}} V_{k-\mathrm{m}+1} \ldots V_{k-1} \nabla \mathrm{f}+\mathrm{s}_{k-\mathrm{m}} \alpha_{k-m} \]

得：

\[ \begin{aligned} r_{k-\mathrm{m}+\mathrm{i}} &=V_{k-\mathrm{m}+\mathrm{i}} \ldots V_{k-\mathrm{m}} H_{0} V_{k-\mathrm{m}} \ldots V_{k-\mathrm{m}+\mathrm{i}} \nabla \mathrm{f} \\ &+\left(V_{k-\mathrm{m}+\mathrm{i}} \ldots V_{k-\mathrm{m}+1}\right) \mathrm{s}_{k-\mathrm{m}} \alpha_{k-m} \\ &+\left(V_{k-\mathrm{m}+\mathrm{i}} \ldots V_{k-\mathrm{m}+2}\right) \mathrm{s}_{k-\mathrm{m}+1} \alpha_{k-m+1}\\ &+\dots\\ &s_{k-m+1}\alpha_{k-m+i} \end{aligned} \]

$r_{k-1}$即是所求的搜索方向d。

使用LBFGS求解逻辑回归模型代码如下：

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn.preprocessing import StandardScaler
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# 进一步处理数据集和测试集，将输入和输出分割
train.columns=list(['x1','x2','y'])
test.columns=list(['x1','x2','y'])
X_train = np.asarray(train.get(['x1', 'x2']))
y_train = np.asarray(train.get('y'))
X_test = np.asarray(test.get(['x1', 'x2']))
y_test = np.asarray(test.get('y'))
# 使用 sklearn 的 LogisticRegression 作为模型
# 其中有 penalty，solver，multi_class 几个比较重要的参数，不同的参数有不同的准确率
model = LogisticRegression(solver='newton-cg')
# newton-cg sag lbfgs liblinear


# 对数据进行标准化
ss = StandardScaler()
X_train = ss.fit_transform(X_train) 
X_test = ss.fit_transform(X_test)
# 拟合
model.fit(X_train, y_train)

# 预测测试集
predictions = model.predict(X_test)

# 打印准确率
print('测试集准确率：', accuracy_score(y_test, predictions))

weights = np.column_stack((model.intercept_, model.coef_)).transpose()
#print(weights)

三、绘制ROC曲线和PR曲线

该部分出现的英语缩写：
TP: True Positive
FP: False Positive
FN: False Negative
TN: True Negative
P: Precision
R: Recall
TPR: True Positive Rate
FPR: False Positive Rate

3.1 ROC曲线

3.1.1介绍

ROC全称是“受试者工作特征”(Receiver Operating Characteristic)曲线，它源于“二战”中用于敌机检测的雷达信号分析技术，二十世纪六七十年代开始被用于一些心理学、医学检测应用中，此后被引入机器学习领域，用来评判分类、检测结果的好坏。因此，ROC曲线是非常重要和常见的统计分析方法。

为了绘制ROC曲线，我们需要计算出两个重要量的值（

TPR、FPR），分别以它们为横、纵坐标作图。其中的TP、FP、TN、FN来自于混淆矩阵，且TP+FP+TN+FN=样本总数。

\[ TPR=\frac{TP}{TP+FN} \]

\[ FPR=\frac{FP}{FP + TN} \]

3.1.2画图流程

给定m+个正例和m-个负例，根据学习器预测结果对样例进行排序
然后把分类阈值设为最大，即把所有样例均预测为反例，此时真正例率和假正例率均为0，在坐标(0,0)处标记一个点
将分类阈值依次设为每个样例的预测值，即依次将每个样例划分为正例，设前一个标记点坐标为(x,y)，当前若为真正例，则对应标记点的坐标为$\left ( x,y+\frac{1}{m^{+}} \right )$；当前若为假正例，则对应标记点的坐标为$\left ( x+\frac{1}{m^{-}},y \right )$
最后用线段连接相邻点

3.1.3 AUC分析

ROC曲线下方的面积也有着重要意义（英语：Area under the Curve of ROC (AUC ROC)），其意义是：

因为是在1x1的方格里求面积，AUC必在0~1之间。
假设阈值以上是正例，以下是反例；
简单说：AUC值越大的分类器，正确率越高。

从AUC判断分类器（预测模型）优劣的标准：

AUC = 1，是完美分类器，采用这个预测模型时，存在至少一个阈值能得出完美预测。绝大多数预测的场合，不存在完美分类器。
0.5 < AUC < 1，优于随机猜测。这个分类器（模型）妥善设置阈值的话，能有预测价值。
AUC = 0.5，跟随机猜测一样（例：丢铜板），模型没有预测价值。
AUC < 0.5，比随机猜测还差；但只要总是反预测而行，就优于随机猜测。

假设ROC曲线由为{ ( x1,y1 ),⋯,( xN′,yN′ ) }的点按需连接而成且有x₁=0,x_N'=1，则AUC可估算为：

\[ AUC=\frac{1}{2} \sum_{j=1}^{N{}'-1} \left ( x_{j+1}-x_{j} \right ) \left ( y_{j+1}+y_{j} \right ) \]

如图即为使用本次作业所提供数据绘制的ROC曲线。由于测试样例有限，所以仅能获得有限个（真正例率，假正例率）坐标对，无法产生光滑的ROC曲线；由此计算得到的AUC的值为0.9648，可以得知该模型的性能较优。

完整代码如下：

def draw_roc(confidence_scores, data_labels):
    #真正率，假正率
    fpr, tpr, thresholds = roc_curve(data_labels, confidence_scores)
    plt.figure()
    plt.grid()
    plt.title('ROC Curve')
    plt.xlabel('FPR')
    plt.ylabel('TPR')
 
    from sklearn.metrics import auc
    auc=auc(fpr, tpr) #AUC计算
    plt.plot(fpr,tpr,'k--', label = 'roc_curve(AUC=%0.4f)' % auc)
    plt.legend()
    plt.show()

3.2 PR曲线

3.2.1介绍

PR曲线全称为查准率-查全率曲线，查准率P与查全率R分别定义为：

\[ P=\frac{TP}{TP+FP}，R=\frac{TP}{TP+FN} \]

查准率和查全率是一对矛盾的度量。一般来说，查准率高时，查全率往往偏低；而查全率高时，查准率往往偏低。

3.2.2画图流程

绘制PR曲线的流程与ROC曲线类似，我们需要根据学习器的预测结果按正例可能性大小对样例进行排序，再逐个样本的选择阈值，在该样本之前的都属于正例，该样本之后的都属于负例。每一个样本作为划分阈值时，都可以计算对应的precision和recall，那么就可以以此绘制曲线。

3.2.3 AP分析

其中平衡点是曲线上“查准率=查全率”时的取值，可用于度量PR曲线有交叉的分类器性能高低。与AUC类似，PR曲线下方面积也有重要意义。PR曲线下的面积称之为AP(Average Precision)，通常来说一个越好的分类器，AP值越高。

对于连续的PR曲线，有：

\[ AP=\int_{0}^{1} p\left ( r \right ) \mathrm{d}r \]

但由于曲线可能出现不可导的部分，故我们常常求其近似值：

\[ p_{\text {interp }}(r)=\max _{\tilde{r} \geq r} p(\tilde{r}) \]

对于离散的PR曲线，有：

\[ \mathrm{AP}=\sum_{k=1}^{n} p(k) \Delta r(k) \]

另外PR曲线平衡点更用常用的是F1度量：

\[ F 1=\frac{2 \times P \times R}{P+R}=\frac{2 \times T P}{\text { 样例总数 }+T P-T N} \]

比F1度量更一般的形式是F_β：

\[ F_{\beta}=\frac{\left(1+\beta^{2}\right) \times P \times R}{\left(\beta^{2} \times P\right)+R} \]

β=1：标准F1
β>1：偏重查全率（逃犯信息检索）
β<1：偏重查准率（商品推荐系统）

如图即为使用本次作业所提供数据绘制的PR曲线。在现实任务中，PR曲线是非单调、不平滑的，在很多局部有上下波动；由此计算得到的AP的值为0.9751，可以得知该模型的性能较优。

完整代码如下：

def draw_pr(confidence_scores, data_labels):
    plt.figure()
    plt.title('PR Curve')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.grid()
 
    #精确率，召回率，阈值
    precision,recall,thresholds = precision_recall_curve(data_labels,confidence_scores)
 
    from sklearn.metrics import average_precision_score
    AP = average_precision_score(data_labels, confidence_scores) # 计算AP
    plt.plot(recall, precision,'k--', label = 'pr_curve(AP=%0.4f)' % AP)
    plt.legend()
    plt.show()

四、总结模型训练过程中的收获

4.1加深了对逻辑斯蒂回归的理解

4.1.1简述对模型的理解：

因为线性回归模型产生的预测值是一系列实值。为了使得输出的预测结果变成分类所需的0和1，我们需要在线性回归的基础式子外再套一个函数将其输出变成0和1，又要求该函数单调可微，所以我们引入logistic函数，将输出的预测结果成功转为概率值。这样，逻辑斯蒂回归模型被成功应用于解决分类模型。

4.1.2关于算法的择优：

在代码实现过程中，我们最开始使用的是梯度下降法，但是迭代速度较慢，拟合效果不是很好；之后我们选择了牛顿法，但是因为计算海森矩阵的复杂度太高，我们选择用一种拟牛顿法——‘L-BFGS’来逼近海森矩阵，最终达到了我们理想的效果。

梯度下降法和牛顿法/拟牛顿法相比，两者都是迭代求解，不过梯度下降法是梯度求解，而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。相对而言，使用牛顿法/拟牛顿法收敛更快。

4.2实现了代码技能的提升

在代码实现过程中，我们调用了机器学习工具包sklearn中的重要函数——LogisticRegression函数，熟悉了它的常用参数及意义，下面以表格形式列出我们在此次模型训练中使用到的参数。

参数	意义	备注
penalty	str类型，可选项有{‘L1’,‘L2’}，用来确定惩罚项的规范。‘newton-cg’，‘sag’和‘lbfgs’仅支持‘L2’惩罚项。	该参数是为了添加惩罚项，避免过拟合，用以提高函数的泛化能力。我们在本次模型训练中使用的是‘L2’。
solver	可选的优化算法有{‘newton-cg’，‘lbfgs’,‘liblinear’,‘sag’}	小数据集中，liblinear是一个好选择，sag和saga对大数据更快；多分类问题中，除了liblinear其它四种算法都可以使用；newton-cg，lbfgs和sag仅能使用L2惩罚项；我们经过对比，选择的算法是lbfgs。
multi_class	str类型，可选参数有{‘ovr’，‘multinomial’} 如果是二元分类问题则两个选项一样，如果是多元分类则ovr将进行多次二分类，分别为一类别和剩余其它所有类别; multinomial则分别进行两两分类，需要T(T-1)/2次分类。	在多分类中，ovr快，精度低; multinomial慢，精度高。