（Death）PRML学习笔记 4

第四章 线性分类（逻辑回归）

一、概念

1. 分类目标

• 分类的目标是将输入变量 x 分到 K 个离散的类别 C_k 中的某一类
• 输入空间被划分为不同的决策区域，决策区域的边界称为决策边界或决策面
• 线性分类模型指的是决策面为输入向量 x 的线性函数，因此若输入的 x 为 D 维，决策面就是 (D-1) 维的超平面
• 若数据集可以被线性决策面分类，则称该数据集是线性可分的

2. 问题表示

• 在分类问题中，目标值（target）表示类别标签。如二分类中，常用二元表示法，即 t 可以取 0 或 1，来表示 C₁ 或 C₂ 类。这样表示可以将 t 理解为分到 C₁ 的概率，而概率只有 0 或 1；对于多分类问题，常用 1-of-k 编码方式
• 上一章线性回归的模型的形式为：
  $$y\left(x\right)=w^{T}x+w_0$$
  但到了线性分类问题中，我们想得到的是离散的类别标签或是属于某类的概率值，而非连续的 y 值，因此我们可以推广此模型：
  $$y\left(x\right)=f\left(w^{T}x+w_0\right)$$
  其中，f(·) 被称为激活函数，且是一个非线性函数（f 的反函数称为链接函数）；因此决策面依然是 x 的线性函数，但 y 不再是参数 w 的线性函数

二、判别函数

1. 二分类

• 线性判别函数最简单的形式是输入向量的线性函数，即：
  $$y\left(x\right)=w^{T}x+w_0$$
  其中 w 被称为权向量，w₀ 称为偏置，偏置的相反数有时被称为阈值
• 对于一个 x，若 y(x) ≥ 0，那么它被分到 C₁ 中，否则被分到 C₂ 中；因此决策面方程为 y(x) = 0
• 考虑两个点 x_A、x_B 都位于决策面上，则
  $$w^T\left(x_A-x_B\right)=0$$
  因此 w 与决策面正交，即 w 确定了决策面的方向
  若 x 为决策面内一点，则原点到决策面的垂直距离为：
  $$\frac{w^{T}x}{||w||}=-\frac{w_0}{||w||}$$
  可以看到，w₀ 确定了决策面的位置
• 此外，若设某点 x 到决策面的距离为 r，那么通过下列推导可求出 r 的表达式：（x_⊥ 是决策面上的垂足）

$$\begin{gathered} x=x_{\perp }+r\frac{w}{||w||} \\ {w}^{T}x+w_0=w^T{x}_{\perp }+r||w||+w_0 \\ r=\frac{y\left(x\right)}{||w||} \end{gathered}$$

• 与上一章相同，我们可以引入 x₀ = 1，使得 y(x) = w^Tx
  此时决策面为 D 维，且这个超平面会穿过 D+1 维的扩展输入空间的原点

2. 多分类

• 多分类的判别函数如果简单的使用二分类的思想，即 “一对多” 的二分类或 “两两之间” 分别二分类，都会导致输入空间中有无法分类的区域，因此可以引入 K 类判别函数
• K 类判别函数由 K 个线性函数组成：
  $$y_k\left(x\right)=w_k^{T}x+w_{k0}$$
  若对于所有的 j ≠ k，都有 y_k(x)＞y_j(x)，则把 x 分到 k 类；类别 C_k 与 C_j 的决策面为 y_k(x) = y_j(x)，且对应于一个 D - 1 维超平面：
  $${\left(w_k-w_j\right)}^{T}x+\left(w_{k0}-w_{j0}\right)=0$$
  这样的判别函数的决策区域总是单连通的，且是凸的，因此不会出现无法分类的区域

3. 最小平方方法

• 由于最小平方方法是在给定输入向量的情况下，近似了目标值的条件期望 E[t|x]，条件期望通常由后验概率给出，但通常这些概率很难近似；另外，最小平方方法对应于高斯条件分布假设下的最大似然法，但数据按不同类划分情况通常并不是高斯分布

• 最小平方方法对离群点（异常点）比较敏感：（由于异常点离决策面过远，导致平方和误差函数惩罚了正确的划分）
  

  绿色线为 logistic 回归模型的决策面，红色线为最小平方方法的决策面

• 虽然最小平方方法对于任意的 x，y(x) 的和为 1，但每个 y(x) 并未被限制在 (0,1) 之间，因此无法表示成概率的形式

• 由于决策面必须交于一点，因此如下图所示的数据分类就无法使用最小平方方法：
  

  左图是最小平方方法，右图是 logistic 回归方法

4. Fisher线性判别函数（二分类）

• Fisher 线性判别的核心思想是降低维度，假设我们有一个 D 维输入向量 x，然后使用下式投影到一维（一条线上）：
  $$y=w^{T}x$$
  若在 y 上设置一个阈值 - w₀，把 y ≥ - w₀ 的样本分到 C₁ 类，其余的分到 C₂ 类，即可得到一个标准的线性分类器
  由于投影可能造成数据重叠，因此通过调整 w 的值，可以选择一个让类别之间分开最大的投影方式
• 二分类情况下，若 C₁ 有 N₁ 个点，C₂ 有 N₂ 个点，则两类的均值向量为：
  $${\mathbf{m}}_1=\frac{1}{N_1}\sum _{n\in C_1}{x}_n,{\mathbf{m}}_2=\frac{1}{N_2}\sum _{n\in C_2}{x}_n$$
  则最简单的度量类别之间分开程度的方式就是投影后类别均值间的距离：
  $$\begin{gathered} m_2-m_1={\mathbf{w}}^T\left({\mathbf{m}}_2-{\mathbf{m}}_1\right) \\ {m}_k={\mathbf{w}}^{\mathbf{T}}{\mathbf{m}}_{\mathbf{k}} \end{gathered}$$
  由于 w 本身有大小，因此可以将 w 限制为单位长度，则可以使用拉格朗日乘数法求解有限制条件的最大值问题
  由于投影后为常量，因此：
  $$\mathbf{w}\propto \left({\mathbf{m}}_2-{\mathbf{m}}_1\right),\text{由于}{\mathbf{w}}^{\mathbf{T}}\mathbf{w}=||\mathbf{w}|{|}^2\left(\text{常数}\right)$$
• 为了减少重叠，引入类内方差，若使得每个类别的类内方差都较小，就会减少重叠：
  $$s_k^2=\sum _{n\in C_k}{\left(y_n-m_k\right)}^2$$
  其中 y_n = w^Tx_n，总的类内方差就可以把各个类的类内方差加和
• Fisher准则定义为：类间方差与类内方差的比值，即：
  $$J\left(w\right)=\frac{{\left(m_2-m_1\right)}^2}{s_1^2+s_2^2}$$
  用 w 将各值替换可看出决策函数对 w 的依赖：
  $$\begin{gathered} J\left(w\right)=\frac{w^T{S}_{B}w}{w^T{S}_{W}w} \\ {S}_B=\left(m_2-m_1\right){\left(m_2-m_1\right)}^T \\ {S}_W=\sum _{n\in C_k}\left(x_n-m_k\right){\left(x_n-m_k\right)}^T,k=1,2 \end{gathered}$$
  S_B 称为类间协方差矩阵，S_W 称为类内协方差矩阵
• 若将上式关于 w 求导，则 J(w) 取得最大值的条件为：
  $$\left({\mathbf{w}}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{B}}\mathbf{w}\right){\mathbf{S}}_{\mathbf{W}}\mathbf{w}=\left({\mathbf{w}}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{W}}\mathbf{w}\right){\mathbf{S}}_{\mathbf{B}}\mathbf{w}$$
  由于我们只关心 w 的方向，因此可以将标量忽略；式中括号内的两项可以看作是平方项标量，S_B 本身也可看作平分项标量，因此可忽略。又因为 w 与 (m₂ - m₁) 成正比，则可将右侧 w 替换；最后两边再同乘 S_W^-1，得：
  $$\mathbf{w}\propto {\mathbf{S}}_{\mathbf{W}}^{-1}\left({\mathbf{m}}_2-{\mathbf{m}}_1\right)$$
• 对于二分类问题，Fisher 准则可以看作是最小平方方法的一个特例

5. Fisher线性判别函数（多分类）

• 多分类的 Fisher 可以由二分类进行推广到 K 分类，类内协方差矩阵 S_W 容易推广：
  $$\begin{gathered} S_W=\sum _{k=1}^K{S}_k \\ {S}_k=\sum _{n\in C_k}\left(x_n-m_k\right){\left(x_n-m_k\right)}^T \\ {m}_k=\frac{1}{N_k}\sum _{n\in C_k}{x}_n \end{gathered}$$
• 想得到类间协方差矩阵 S_B 需要引入整体的协方差矩阵 S_T：
  $$\begin{gathered} S_T=\sum _{n=1}^N\left(x_n-m\right){\left(x_n-m\right)}^T \\ m=\frac{1}{N}\sum _{n=1}^N{x}_n=\frac{1}{N}\sum _{k=1}^K{N}_k{m}_k \end{gathered}$$
  m 为全体数据的均值，且：
  $$S_T=S_W+S_B$$
  这样就可以推出类间协方差矩阵 S_B：
  $$S_B=\sum _{k=1}^K{N}_k\left(m_k-m\right){\left(m_k-m\right)}^T$$
• 当投影到 y 空间后，依然可以定义类似的类内和类间协方差矩阵，准则函数也可以与之改变，其中一种选择是：
  $$J\left(W\right)=Tr\left\{s_W^{-1}{s}_B\right\}=Tr\left\{{\left(W^T{S}_{W}W\right)}^{-1}\left(W^T{S}_{B}W\right)\right\}$$
• S_B 由 K 个矩阵的和组成，而每个矩阵都是两个向量的外积，因此秩都为 1；又可以证明这些矩阵中只有 (K - 1) 个是相互独立的，因此加和以后得到的 S_B 的秩最大就是 (K - 1)，即最多有 (K - 1) 个非零特征值

6. 感知器算法

• 感知器算法对应于二分类的模型，首先输入向量 x 进行一个固定的非线性变换得到一个特征向量 Φ(x)，然后构造一个线性模型：
  $$\begin{gathered} y\left(x\right)=f\left(w^T\varphi \left(x\right)\right) \\ f\left(a\right)=\{\begin{array}{l}+1,a⩾0\\ -1,a<0\end{array} \end{gathered}$$
  其中 f(·) 为非线性激活函数

• 由于感知器算法模型是分段常函数，因此用一般的误差函数会有不连续和梯度为 0 等问题，于是引入另外的一个误差函数，称为感知器准则
  我们希望找到一个权向量 w 使得对于 C₁ 中的模式 x_n 都有 w^TΦ(x_n)＞0，而对于 C₂ 中的刚好相反；因为目标变量 t 是 ±1，所以对于所有的 x_n，都希望有 w^TΦ(x_n)·t_n＞0
  于是，对于正确分类的模式，感知器准则为零误差；而对于误分类的模式，因为 w^TΦ(x_n)·t_n＜0，因此我们可以最小化 -w^TΦ(x_n)·t_n 来让误分类的情况最小
  因此感知器准则为：
  $$E_P\left(w\right)=-\sum _{n\in M}{w}^T\varphi \left(x_n\right)t_n$$
  其中 M 表示所有误分类模式的集合。由于误差函数对于正确分类的 x_n 值为0，对于错误分类的为 w 的线性函数，因此误差函数是分段线性的

• 对这个误差函数使用随机梯度下降，w 的变化为：
  $$w^{\left(\tau +1\right)}=w^{\left(\tau \right)}-\eta \nabla E_P\left(w\right)=w^{\left(\tau \right)}+\eta \varphi \left(x_n\right)t_n$$
  η 为学习率，由于将 w 乘一个常数后，y(x,w) 分类效果不变，因此可以令 η 为 1，而式子依然具有一般性

• 感知器学习算法可以简单的表示为：对于训练样本进行循环处理，对于每个 x_n，我们都计算感知器函数。若模式正确分类，则 w 不变；若分类错误，则对于类别 C₁，我们将 w 加上 Φ(x_n)，对于类别 C₂，我们将 w 减去 Φ(x_n)，如下图：
  

  黑色线为决策面，黑色箭头为每次的 w，且指向被分为红色的区域；每次寻找误分类数据，将 w 与误分类数据向量叠加，不断更新，最终得到合理的 w
  不断更新过程中，我们可以证明，每一次权值更新过程中，当前错误分类的 x_n 对误差函数的贡献会逐渐减小；但同时其他样本也会改变，可能会使正确分类的变成误分类，因此感知器学习规则并不保证每次都能减小整体的误差函数

三、概率生成式模型

1. 概念

• 这种方法中，我们首先对类条件概率 p(x|C_k) 和先验概率 p(C_k) 建模，之后通过贝叶斯定理计算后验概率 p(C_k|x)
• 首先考虑二分类情况，C₁ 的后验概率可以写成：
  $$\begin{gathered} p\left(C_1|x\right)=\frac{p\left(x|C_1\right)p\left(C_1\right)}{p\left(x|C_1\right)p\left(C_1\right)+p\left(x|C_2\right)p\left(C_2\right)}=\frac{1}{1+\exp \left(-a\right)}=\sigma \left(a\right) \\ a=\ln \frac{p\left(x|C_1\right)p\left(C_1\right)}{p\left(x|C_2\right)p\left(C_2\right)} \end{gathered}$$
  其中 σ(a) 为 logistic sigmoid 函数，满足对称性：
  $$\sigma \left(-a\right)=1-\sigma \left(a\right)$$
  它的反函数为：（logit 函数 / log odds 函数）
  $$a=\ln \left(\frac{\sigma }{1-\sigma }\right)$$
  它表示两类概率比值的对数 $ \ln \left[ \frac{p\left( C_1|x \right)}{p\left( C_2|x \right)} \right] $
• 当 k＞2时，多分类情况下的后验概率为：
  $$\begin{gathered} p\left(C_k|x\right)=\frac{p\left(x|C_k\right)p\left(C_k\right)}{{\sum }_j^{}p\left(x|C_j\right)p\left(C_j\right)}=\frac{\exp \left(a_k\right)}{{\sum }_j^{}\exp \left(a_j\right)} \\ {a}_k=\ln \left(p\left(x|C_k\right)p\left(C_k\right)\right) \end{gathered}$$
  它被称为归一化指数（softmax 函数）
  若对于所有的 j ≠ k，都有 a_k≫a_j，那么 p(C_k|x) 约等于 1，而 p(C_j|x) 约等于 0

2. 输入变量连续情况

• 假设类条件概率密度为高斯分布，求解后验概率
  假定所有类别的协方差相同，这样类别 C_k 的类条件概率为：
  $$p\left(x|C_k\right)=\frac{1}{{\left(2\pi \right)}^{\frac{D}{2}}}\frac{1}{|\Sigma {|}^{\frac{1}{2}}}\exp \left\{-\frac{1}{2}{\left(x-{\mu }_k\right)}^T{\Sigma }^{-1}\left(x-{\mu }_k\right)\right\}$$

• 若考虑两类情况：
  $$\begin{gathered} p\left(C_1|x\right)=\sigma \left(w^{T}x+w_0\right) \\ w={\Sigma }^{-1}\left({\mu }_1-{\mu }_2\right) \\ {w}_0=-\frac{1}{2}{\mu }_1^T{\Sigma }^{-1}{\mu }_1+\frac{1}{2}{\mu }_2^T{\Sigma }^{-1}{\mu }_2+\ln \frac{p\left(C_1\right)}{p\left(C_2\right)} \end{gathered}$$
  可以看到高斯概率密度中 x 的二次型消失了，这是由于我们假设协方差矩阵相同，最终得到了 x 的线性函数的 sigmoid 函数；若不假设协方差矩阵相同，则判别函数会变成二次，如下图：
  

  左图为三个类别的类条件概率密度，右图为各自的后验概率密度和决策边界，可见左侧为线性决策面，右侧两个为二次决策面
  因此决策面为线性的，又可以看到先验概率密度 p(C_k) 仅出现在 w₀ 中，因此先验改变会使决策面平移

• 若考虑 K 个分类的情况：
  $$\begin{gathered} a_k\left(x\right)=w_k^{T}x+w_{k0} \\ {w}_k={\Sigma }^{-1}{\mu }_k \\ {w}_{k0}=-\frac{1}{2}{\mu }_k^T{\Sigma }^{-1}{\mu }_k+\ln p\left(C_k\right) \end{gathered}$$

3. 最大似然解

• 首先考虑两类的情况，假设每个类别都有一个高斯类条件概率，且协方差矩阵相同
  设数据集 {x_n,t_n}，其中 n = 1,2,…,N；t_n = 1 表示 C₁，t_n = 0 表示 C₂
  最后把先验概率记作 p(C₁) = π，则 p(C₂) = 1 - π，则：
  $$\begin{gathered} p\left(x_n,C_1\right)=p\left(C_1\right)p\left(x_n|C_1\right)=\pi \mathcal{N}\left(x_n|{\mu }_1,\Sigma \right) \\ p\left(x_n,C_2\right)=p\left(C_2\right)p\left(x_n|C_2\right)=\left(1-\pi \right)\mathcal{N}\left(x_n|{\mu }_2,\Sigma \right) \end{gathered}$$
• 因此似然函数为：
  $$p\left(t,X|\pi ,{\mu }_1,{\mu }_2,\Sigma \right)=\prod _{n=1}^N{\left[\pi \mathcal{N}\left(x_n|{\mu }_1,\Sigma \right)\right]}^{t_n}{\left[\left(1-\pi \right)\mathcal{N}\left(x_n|{\mu }_2,\Sigma \right)\right]}^{1-t_n}$$
  其中 t = (t₁,…,t_N)^T，取对数为：
  ln ⁡ p = ∑ n = 1 N { t n ln ⁡ [ π N ( x n ∣ μ 1 , Σ ) ] + ( 1 − t n ) ln ⁡ [ ( 1 − π ) N ( x n ∣ μ 2 , Σ ) ] } \ln p=\sum_{n=1}^N{\left\{ t_n\ln \left[ \pi \mathcal{N}\left( x_n|\mu _1,\Sigma \right) \right] +\left( 1-t_n \right) \ln \left[ \left( 1-\pi \right) \mathcal{N}\left( x_n|\mu _2,\Sigma \right) \right] \right\}} lnp=n=1∑N​{tn​ln[πN(xn​∣μ1​,Σ)]+(1−tn​)ln[(1−π)N(xn​∣μ2​,Σ)]}
• 首先考虑关于 π 的最大化，上述对数似然函数中与 π 相关的项为：
  ∑ n = 1 N { t n ln ⁡ π + ( 1 − t n ) ln ⁡ ( 1 − π ) } \sum_{n=1}^N{\left\{ t_n\ln \pi +\left( 1-t_n \right) \ln \left( 1-\pi \right) \right\}} n=1∑N​{tn​lnπ+(1−tn​)ln(1−π)}
  对 π 求导，令导数为 0，可得：
  $$\pi =\frac{1}{N}\sum _{n=1}^N{t}_n=\frac{N_1}{N}$$
  其中 N₁ 表示 C₁ 中数据点的个数，因此此式容易推广到多类情况
• 现在考虑关于 μ₁ 的最大化，上述对数似然函数与 μ₁ 相关的项为：
  ∑ n = 1 N { t n ln ⁡ N ( x n ∣ μ 1 , Σ ) } = − 1 2 ∑ n = 1 N t n ( x n − μ 1 ) T Σ − 1 ( x n − μ 1 ) + C \sum_{n=1}^N{\left\{ t_n\ln \mathcal{N}\left( x_n|\mu _1,\Sigma \right) \right\}}=-\frac{1}{2}\sum_{n=1}^N{t_n\left( x_n-\mu _1 \right) ^T\Sigma ^{-1}}\left( x_n-\mu _1 \right) +C n=1∑N​{tn​lnN(xn​∣μ1​,Σ)}=−21​n=1∑N​tn​(xn​−μ1​)TΣ−1(xn​−μ1​)+C
  对 μ₁ 求导，令导数为 0，可得：
  $${\mu }_1=\frac{1}{N_1}\sum _{n=1}^N{t}_n{x}_n$$
  同理对于 μ₂：
  $${\mu }_2=\frac{1}{N_2}\sum _{n=1}^N\left(1-t_n\right)x_n$$
  μ₁ 和 μ₂ 都是属于对应类别的输入向量 x_n 的均值
• 最后考虑协方差矩阵 Σ 的最大化，上述对数似然函数与之相关的项为：

$$-\frac{1}{2}\sum _{n=1}^N{t}_n\ln |\Sigma |-\frac{1}{2}\sum _{n=1}^N{t}_n{\left(x_n-{\mu }_1\right)}^T{\Sigma }^{-1}\left(x_n-{\mu }_1\right)-\frac{1}{2}\sum _{n=1}^N\left(1-t_n\right)\ln |\Sigma |-\frac{1}{2}\sum _{n=1}^N\left(1-t_n\right){\left(x_n-{\mu }_2\right)}^T{\Sigma }^{-1}\left(x_n-{\mu }_2\right)$$
可以简写为：
$$-\frac{N}{2}\ln |\Sigma |-\frac{N}{2}Tr\left\{{\Sigma }^{-1}S\right\}$$
其中，
$$\begin{gathered} S=\frac{N_1}{N}{S}_1+\frac{N_2}{N}{S}_2 \\ {S}_1=\frac{1}{N_1}\sum _{n\in C_1}{\left(x_n-{\mu }_1\right)\left(x_n-{\mu }_1\right)}^T \\ {S}_2=\frac{1}{N_2}\sum _{n\in C_2}{\left(x_n-{\mu }_2\right)\left(x_n-{\mu }_2\right)}^T \end{gathered}$$

• 求最大似然解，可以得到 Σ = S，它表示对一个与两类都有关系的协方差矩阵求加权平均值：
  $$\Sigma =\frac{1}{N}\sum _{k=1,2}\sum _{n\in C_k}{\left(x_n-{\mu }_k\right)\left(x_n-{\mu }_k\right)}^T$$
  这个结果容易推广到 K 类问题，其中我们假定每个类条件概率密度都是高斯分布，且协方差矩阵相同
  注意，此方法对于离群点依然没有鲁棒性，因为高斯的最大似然估计是没有鲁棒性的

4. 输入变量离散情况

• 当特征值离散时，先考虑二元特征值 x_i 为 0 或 1 的情况，若有 D 个输入，则概率分布会对应于一个大小为 2^D 的矩阵，包含 2^D - 1 个独立变量（要满足加和为1），因此这会随着特征的数量指数增长
• 因此这里我们做出朴素贝叶斯的假设，假设特征值相互独立，条件为 C_k，则类条件分布为：
  $$p\left(x|C_k\right)=\prod _{i=1}^D{\mu }_{ki}^{x_i}{\left(1-{\mu }_{ki}\right)}^{1-x_i}$$
  其中对于每个类别，都有 D 个独立的参数；用 sigmoid 函数表示的话，则有：
  a k ( x ) = ∑ i = 1 D { x i ln ⁡ μ k i + ( 1 − x i ) ln ⁡ ( 1 − μ k i ) } + ln ⁡ p ( C k ) a_k\left( x \right) =\sum_{i=1}^D{\left\{ x_i\ln \mu _{ki}+\left( 1-x_i \right) \ln \left( 1-\mu _{ki} \right) \right\}}+\ln p\left( C_k \right) ak​(x)=i=1∑D​{xi​lnμki​+(1−xi​)ln(1−μki​)}+lnp(Ck​)
  这是输入变量 x_i 的线性函数

5. 指数族分布

• 无论是服从高斯分布还是离散的输入，后验概率都是由一般的线性模型或 sigmoid(K=2) / softmax(K>2) 激活函数所给出的
• 而通过假定类条件概率密度 p(x|C_k) 是指数族分布的成员，则可以将上述结果写成一个统一的一般式：
  $$p\left(x|{\lambda }_k\right)=h\left(x\right)g\left({\lambda }_k\right)\exp \left\{{\lambda }_k^{T}u\left(x\right)\right\}$$
  我们考虑 u(x) = x 的情况，并引入一个缩放参数 s：（每个类别都有自己的参数向量 λ_k，但各类的 s 相同）
  $$p\left(x|{\lambda }_k,s\right)=\frac{1}{s}h\left(\frac{1}{s}x\right)g\left({\lambda }_k\right)\exp \left\{\frac{1}{s}{\lambda }_k^{T}x\right\}$$
• 对于二分类问题，我们将类条件概率密度代入之前的 sigmoid 表达式，可得：
  $$a\left(x\right)=\frac{1}{s}{\left({\lambda }_1-{\lambda }_2\right)}^{T}x+\ln g\left({\lambda }_1\right)-\ln g\left({\lambda }_2\right)+\ln p\left(C_1\right)-\ln p\left(C_2\right)$$
  类似地，对于 K 类问题，可得：
  $$a_k\left(x\right)=\frac{1}{s}{\lambda }_k^{T}x+\ln g\left({\lambda }_k\right)+\ln p\left(C_k\right)$$
  均为 x 的线性函数

四、概率判别式模型

1. 判别式与生成式建模的区别

• 生成式建模的一个例子是寻找类条件概率和先验概率，使用贝叶斯定理得到后验概率；而判别式模型是显式的使用一般的线性模型的函数形式，然后使用最大似然法直接确定它的参数
• 寻找这样的解有一个高效的算法：迭代重加权最小平方法（IRLS）

2. 固定基函数

• 若选择一个基函数向量 Φ(x) 对输入变量 x 进行一个非线性变换，在 Φ 空间中决策面是线性的，因此对应在原始 x 空间中决策面为非线性，如下图：
  

  左图为 x 空间中的分布，右图为 Φ 空间的分布，黑色线为决策面

• 对于许多实际问题，类条件概率密度之间有着很大重叠，因此后验概率不能等于 0 或 1；非线性变换 Φ(x) 也不能消除这些重叠，甚至有时会增加重叠，但恰当的非线性变换可以使后验概率建模过程更简单，从而便于决策

3. logistic 回归

• 首先还是讨论二分类问题，之前我们将后验概率写成了 sigmoid 的形式，在一些一般性的假设条件下，C₁ 的后验概率可以写成作用在特征向量 Φ 上的线性函数的 logistic sigmoid 函数的形式：
  $$\begin{gathered} p\left(C_1|\varphi \right)=y\left(\varphi \right)=\sigma \left(w^T\varphi \right) \\ p\left(C_2|\varphi \right)=1-p\left(C_1|\varphi \right) \end{gathered}$$
• 对于 M 维的特征空间 Φ，这个模型有 M 个参数，因此是线性关系增长；若使用最大似然法调节高斯类条件概率，有 2M 个参数描述均值，M(M+1)/2 个参数描述协方差矩阵，1 个参数描述先验，因此共有 M(M+5)/2+1 个参数量，与 M 呈二次关系增长
• 现在使用最大似然方法确定 logistic 回归模型的参数，似然函数可以写成：
  $$\begin{gathered} p\left(t|w\right)=\prod _{n=1}^N{y}_n^{t_n}{\left(1-y_n\right)}^{1-t_n} \\ {y}_n=p\left(C_1|{\varphi }_n\right)=\sigma \left(a_n\right) \\ {a}_n=w^T\varphi \left(x_n\right) \end{gathered}$$
  通过取似然函数的负对数，定义一个误差函数：
  E ( w ) = − ln ⁡ p ( t ∣ w ) = − ∑ n = 1 N { t n ln ⁡ y n + ( 1 − t n ) ln ⁡ ( 1 − y n ) } E\left( w \right) =-\ln p\left( t|w \right) =-\sum_{n=1}^N{\left\{ t_n\ln y_n+\left( 1-t_n \right) \ln \left( 1-y_n \right) \right\}} E(w)=−lnp(t∣w)=−n=1∑N​{tn​lnyn​+(1−tn​)ln(1−yn​)}
  这个误差函数称为交叉熵误差函数（cross-entropy）
  两侧均对 w 求导，得到误差函数的梯度：（sigmoid 的导数为 σ(1-σ)）
  $$\nabla E\left(w\right)=\sum _{n=1}^N\left(y_n-t_n\right)\varphi \left(x_n\right)$$
  它的形式与线性回归模型中的平方和误差函数的梯度形式完全相同
• 同样的，我们可以提出一种顺序算法，每次只出现一个模式 x_n 使权向量进行更新
• 值得注意的是，最大似然方法对于线性可分的数据集会产生过拟合现象，最大似然解出现在 σ = 0 即 w^TΦ = 0 的情况下，通过引入先验概率或加入一个正则化项可以避免这种过拟合问题

4. 迭代重加权最小平方（IRLS）

• 上一章的线性回归模型中，我们假设的高斯模型的最大似然解有解析解，而 logistic 回归模型不再有解析解，因为它是非线性函数
• 误差函数可以用一种迭代方法求出最小值，这种方法基于牛顿-拉弗森迭代最优化框架（Newton-Raphson），使用了对数似然函数的局部二次近似，权值的更新表达式为：
  $$w^{new}=w^{old}-H^{-1}\nabla E\left(w\right)$$
  其中 H 为海瑟矩阵（Hessian Matrix），它的元素由 E(w) 关于 w 的二阶导数组成
• 若将此方法用于线性回归模型上，误差函数为平方和误差函数，则：
  $$\begin{gathered} \nabla E\left(w\right)=\sum _{n=1}^N\left(w^T\varphi \left(x_n\right)-t_n\right)\varphi \left(x_n\right)={\Phi }^T\Phi w-{\Phi }^{T}t \\ H=\nabla \nabla E\left(w\right)={\Phi }^T\Phi \end{gathered}$$
  其中等式右侧的 Φ 是 N×M 设计矩阵，第 n 行为 Φ^T(x_n)，因此：
  $$w^{new}=w^{old}-{\left({\Phi }^T\Phi \right)}^{-1}\left({\Phi }^T\Phi w^{old}-{\Phi }^{T}t\right)={\left({\Phi }^T\Phi \right)}^{-1}{\Phi }^{T}t$$
  可以看到这就是标准的最小平方解，由于此处误差函数为二次，因此只用一次迭代就得到了精确解
• 将此方法用于 logistic 回归模型上，误差函数为交叉熵误差函数，则：
  $$\begin{gathered} \nabla E\left(w\right)=\sum _{n=1}^N\left(y_n-t_n\right)\varphi \left(x_n\right)={\Phi }^T\left(y-t\right) \\ H=\nabla \nabla E\left(w\right)={\Phi }^{T}R\Phi \end{gathered}$$
  其中，R 为 N×N 的对角矩阵，元素为：
  $$R_{nn}=y_n\left(1-y_n\right)$$
  可以看到 Hessian 矩阵不再是常量，而是通过权矩阵 R 依赖于 w；由于 Hessian 矩阵是正定的，因此误差函数是 w 的一个凸函数，有唯一的最小值，因此：
  $$\begin{gathered} w^{new}=w^{old}-{\left({\Phi }^{T}R\Phi \right)}^{-1}{\Phi }^T\left(y-t\right)={\left({\Phi }^{T}R\Phi \right)}^{-1}{\Phi }^{T}Rz \\ z={\Phi }^{old}-R^{-1}\left(y-t\right) \end{gathered}$$
  其中 z 为一个 N 维向量

5. 多类别 logistic 回归

• 上一节我们对于多分类的生成式模型给出了后验概率，由特征变量的线性函数的 softmax 变换给出：
  $$\begin{gathered} p\left(C_k|\varphi \right)=y_k\left(\varphi \right)=\frac{\exp \left(a_k\right)}{{\sum }_j^{}\exp \left(a_j\right)} \\ {a}_k=w_k^T\varphi \end{gathered}$$
  之前我们使用最大似然估计类条件概率和先验概率，用贝叶斯定理找到对应的后验概率；这里我们直接使用最大似然方法直接确定参数 w
• 我们首先计算 y_k 关于所以激活 a_j 的导数：
  $$\frac{\partial y_k}{\partial a_j}=y_k\left(I_{kj}-y_j\right)$$
  其中 I_kj 为单位矩阵
  接下来写出似然函数，使用 1-of-K 表达方式：
  $$p\left(T|w_1,...,w_K\right)=\prod _{n=1}^N\prod _{k=1}^{K}p{\left(C_k|{\varphi }_n\right)}^{t_{nk}}=\prod _{n=1}^N\prod _{k=1}^K{y}_{nk}^{t_{nk}}$$
  其中，T 是目标变量的一个 N×K 的矩阵，元素为 t_nk，且：
  $${\varphi }_n=\varphi \left(x_n\right),y_{nk}=y_k\left({\varphi }_n\right)$$
  取负对数，可得：
  $$E\left(w_1,...,w_K\right)=-\ln p\left(T|w_1,...,w_K\right)=-\sum _{n=1}^N\sum _{k=1}^K{t}_{nk}\ln y_{nk}$$
  它被称为多分类问题的交叉熵误差函数
• 现在求误差函数关于 w_j 的梯度：
  $$\nabla w_{j}E\left(w_1,...,w_K\right)=\sum _{n=1}^N\left(y_{nj}-t_{nj}\right){\varphi }_n$$
  可以看出，梯度的函数形式是 (y_nj - t_nj) 与基函数 Φ_n 的乘积，这种形式在之前也出现过，因此这是一个更一般的结果，且可以将此公式用于顺序算法
  而多分类的 IRLS 算法，要求出 M×M 的块组成的 Hessian 矩阵，其中 块 i,j 为：
  $$\nabla w_k\nabla w_{j}E\left(w_1,...,w_K\right)=\sum _{n=1}^N{y}_{nk}\left(I_{kj}-y_{nj}\right){\varphi }_n{\varphi }_n^T$$
  该 Hessian 矩阵是正定的，因此误差函数有唯一的最小值

6. probit 回归

• 我们还是考虑二分类的情况，使用一般的线性模型的框架，即：
  $$\begin{gathered} p\left(t=1|a\right)=f\left(a\right) \\ a=w^T\varphi \end{gathered}$$
  f(·) 为激活函数，对每个输入 Φ_n，我们计算 a_n = w^TΦ_n，然后按照此方式设置 t_n ：若 a_n ≥ θ，则设 t_n = 1；否则设 t_n = 0
  若 θ 的值服从概率分布 p(θ)，则对应的激活函数为：
  $$f\left(a\right)={\int }_{-\infty }^{a}p\left(\theta \right)d\theta$$
  若 p(θ) 是标准高斯概率密度，则累计分布函数（激活函数）为：
  $$\Phi \left(a\right)={\int }_{-\infty }^a\mathcal{N}\left(\theta |0,1\right)d\theta$$
  这被称为逆 probit 函数（inverse probit），形状为 sigmoid 形状；且使用更一般的高斯分布也不会改变模型，因为这样只是等价于对线性系数 w 的缩放
• 逆 probit 函数中出现较多的一个函数为 erf 函数（也称 error 函数（不是误差函数））：
  $$\begin{gathered} erf\left(a\right)=\frac{2}{\sqrt{\pi }}{\int }_0^a\begin{array}{c}\exp \left(-{\theta }^2\right)d\theta \end{array} \\ \Phi \left(a\right)=\frac{1}{2}\left\{1+erf\left(\frac{a}{\sqrt{2}}\right)\right\} \end{gathered}$$
  基于 probit 激活函数的一般线性模型被称为 probit 回归
• 使用最大似然法来确定模型参数，结果与 logistic 回归的结果类似；然而对于离群点，当 x→∞ 时，probit 激活函数为 exp(-x²) 速度衰减，而 logistic sigmoid 函数为 exp(-x) 速度衰减，因此 probit 模型对离群点更加敏感
• logistic 模型和 probit 模型都假设数据点被正确标记，因此错误标记的影响会很容易影响模型；我们可以引入一个概率 ε，它是目标值 t 被翻转到错误值的概率，则数据点 x 的目标值 t 的分布为：
  $$p\left(t|x\right)=\left(1-ϵ\right)\sigma \left(x\right)+ϵ\left(1-\sigma \left(x\right)\right)=ϵ+\left(1-2ϵ\right)\sigma \left(x\right)$$
  其中 σ(x) 是 x 的激活函数，ε 可以提前设定，也可以当成参数进行推断

7. 正则链接函数（Canonical link function）（证明略）

• 假设目标变量的条件分布为指数族分布，对应的激活函数选择正则连接函数，则对 “数据点 n 对误差函数的贡献” 关于参数 w 求导，得到的导数形式为误差 y_n - t_n 与特征向量 Φ_n 的乘积，其中 y_n = w^TΦ_n，这个结果是一般性的结果
• 广义线性模型：y 是输入变量（特征变量）的线性组合的非线性函数：
  $$y=f\left(w^T\varphi \right)$$
  其中 f 为激活函数，f^-1 为链接函数

五、拉普拉斯近似

1. 概念

• 拉普拉斯近似的目标是找到定义在一组连续变量上的概率密度的高斯近似
  首先考虑单一连续变量 z 的情形，假设分布 p(z) 的定义为：
  $$\begin{gathered} p\left(z\right)=\frac{1}{Z}f\left(z\right) \\ Z=\int f\left(z\right)dz \end{gathered}$$
  其中 Z 是归一化系数，假设 Z 的值是未知的
  在拉普拉斯方法中，我们目标是寻找一个高斯近似 q(z)，它的中心位于 p(z) 的众数的位置
• 第一步是寻找 p(z) 的众数，即寻找一个点 z₀ 使得 p’(z₀) = 0，即：
  $$\frac{df\left(z\right)}{dz}{\mid }_{z=z_0}^{}=0$$
  由于高斯分布的对数是变量的二次函数，考虑 lnf(z) 以众数 z₀ 为中心的泰勒展开式：（一次项系数为 0）
  $$\begin{gathered} \ln f\left(z\right)\simeq \ln f\left(z_0\right)-\frac{1}{2}A{\left(z-z_0\right)}^2 \\ A=-\frac{d^2}{d{z}^2}\ln f\left(z\right){\mid }_{z=z_0}^{} \end{gathered}$$
  两侧同时取指数，再使用归一化的高斯分布的标准形式，就可以得到归一化的概率分布 q(z)：
  $$q\left(z\right)={\left(\frac{A}{2\pi }\right)}^{\frac{1}{2}}\exp \left\{-\frac{A}{2}{\left(z-z_0\right)}^2\right\}$$
  这就是拉普拉斯近似，但只有 A＞0 时近似效果较好，即 f’'(z₀)＜0
• 将拉普拉斯近似推广，近似 M 维空间 z 上的概率分布 p(z) = f(z)/Z；同理可得：
  $$q\left(z\right)=\frac{|A{|}^{\frac{1}{2}}}{{\left(2\pi \right)}^{\frac{M}{2}}}\exp \left\{-\frac{1}{2}{\left(z-z_0\right)}^{T}A\left(z-z_0\right)\right\}=\mathcal{N}\left(z|z_0,A^{-1}\right)$$
  当 A 是正定矩阵时，此高斯分布有良好的定义，即 z₀ 一定是一个局部极大值点（不是极小值或鞍点）

2. 模型比较和贝叶斯信息准则（BIC）

• 除了近似概率分布 p(z)，我们还可以近似归一化常数 Z：
  $$Z=\int f\left(z\right)dz\simeq f\left(z_0\right)\frac{{\left(2\pi \right)}^{\frac{M}{2}}}{|A{|}^{\frac{1}{2}}}$$
  而这个形式可以获得模型证据的一个近似
• 考虑一个数据集 D 以及一组模型 {M_i}，模型参数为 {θ_i}，对于每个模型都定义一个似然函数 p(D|θ_i,M_i)，再引入一个参数的先验概率 p(θ_i|M_i)，那么我们想计算的是不同模型的模型证据 P(D|M_i)，则：（为了简化记号，省略 M_i）
  $$p\left(D\right)=\int p\left(D|\theta \right)p\left(\theta \right)d\theta$$
  令 f(θ) = p(D|θ)p(θ)，Z = p(D)，则：
  $$\begin{gathered} \ln p\left(D\right)\simeq \ln p\left(D|{\theta }_{MAP}\right)+\ln p\left({\theta }_{MAP}\right)+\frac{M}{2}\ln \left(2\pi \right)-\frac{1}{2}\ln |A| \\ A=-\nabla \nabla \ln p\left({\theta }_{MAP}|D\right) \end{gathered}$$
  θ_MAP 为后验概率分布中众数位置的 θ 值，A 是负对数后验概率的二阶导数组成的 Hessian 矩阵
  第一项为对数似然，后面三项为 Occam 因子，它对模型的复杂度进行惩罚
• 若假设参数的高斯先验分布比较宽，且 Hessian 矩阵满秩，则可以将上式再近似为：
  $$\ln p\left(D\right)\simeq \ln p\left(D|{\theta }_{MAP}\right)-\frac{1}{2}M\ln N$$
  其中 N 是数据点的总是，M 是 θ 中参数的数量，且省略一些常数项；这被称为贝叶斯信息准则（BIC）或者 Schwarz 准则

六、贝叶斯 logistic 回归

1. 后验概率的高斯近似

• 由于贝叶斯推断无法处理 logistic 回归，因此采用拉普拉斯近似来处理贝叶斯 logistic 回归问题，我们的主要目标就是寻找后验概率分布的一个高斯表示
• 首先选择高斯先验：
  $$\begin{gathered} p\left(w\right)=\mathcal{N}\left(w|m_0,S_0\right) \\ p\left(w|t\right)\propto p\left(w\right)p\left(t|w\right) \end{gathered}$$
  m₀、S₀ 为固定的超参数，t = (t₁,…,t_N)^T；两侧取对数，然后代入先验：
  ln ⁡ p ( w ∣ t ) = − 1 2 ( w − m 0 ) T S 0 − 1 ( w − m 0 ) + ∑ n = 1 N { t n ln ⁡ y n + ( 1 − t n ) ln ⁡ ( 1 − y n ) } + C \ln p\left( w|t \right) =-\frac{1}{2}\left( w-m_0 \right) ^TS_0^{-1}\left( w-m_0 \right) +\sum_{n=1}^N{\left\{ t_n\ln y_n+\left( 1-t_n \right) \ln \left( 1-y_n \right) \right\}}+C lnp(w∣t)=−21​(w−m0​)TS0−1​(w−m0​)+n=1∑N​{tn​lnyn​+(1−tn​)ln(1−yn​)}+C
  其中，$ y_n=\sigma \left( w^T\phi _n \right) $
• 最大化后验分布，得到 w_MAP，它定义了高斯分布的均值；这样协方差矩阵就是负对数似然函数的二阶导数矩阵的逆矩阵：
  $$S_N=S_0^{-1}+\sum _{n=1}^N{y}_n\left(1-y_n\right){\varphi }_n{\varphi }_n^T$$
  则后验分布的高斯近似为：
  $$q\left(w\right)=\mathcal{N}\left(w|w_{MAP},S_N\right)$$

2. 预测分布（证明略）

• 给定一个新的特征向量 Φ(x)，类别 C₁ 的预测分布可以通过对后验概率 p(w|t) 积分得到
  通过推导，我们可以得到预测分布的近似：
  $$\begin{gathered} p\left(C_1|t\right)=\int \sigma \left(a\right)\mathcal{N}\left(a|{\mu }_a,{\sigma }_a^2\right)da \\ a=w^T\varphi \\ {\mu }_a=w_{MAP}^T\varphi \\ {\sigma }_a^2={\varphi }^T{S}_N\varphi \end{gathered}$$
  然而这个积分表示一个高斯分布和一个 logistic sigmoid 函数的卷积，不能解析地求值
• 但通过用 logistic sigmoid 函数和逆 probit 函数可以求一个较好的近似：
  $$\begin{gathered} p\left(C_1|\varphi ,t\right)=\sigma \left[\kappa \left({\sigma }_a^2\right){\mu }_a\right] \\ \kappa \left({\sigma }_a^2\right)={\left(1+\frac{\pi {\sigma }_a^2}{8}\right)}^{-\frac{1}{2}} \end{gathered}$$

七、代码实例

1. fisher 线性判别（以 Iris 数据集为例）

from sklearn import datasets
import matplotlib.pyplot as plt
import numpy as np

# 参考 https://blog.csdn.net/mengjizhiyou/article/details/103309372

# 准备数据
iris = datasets.load_iris()
X_origin = iris.data       # (150,4)的矩阵
# 随机分组验证
ac1 = 0     # 一维平均正确率
ac2 = 0     # 二维平均正确率
for times in range(10):		# 随机取k个作为测试集,多次求平均
    k_size = 30    # 分组大小 设为30
    randomindex = np.random.choice(X_origin.shape[0], size = k_size, replace=False)
    X_test = X_origin[randomindex, :]   # 测试集为(30,4)的矩阵
    X = np.delete(X_origin, randomindex, axis=0)    # 训练集为(120,4)的矩阵
    y_origin = iris.target     # (150,_)的向量
    y_test = y_origin[randomindex]   # 测试集的标签为(30,)的矩阵
    y = np.delete(y_origin, randomindex, axis=0)    # 标签为(120,_)的向量
    names = iris.feature_names
    labels = iris.target_names
    y_c = np.unique(y)  # 数组去重后排序，为了查看数据中共几类数据
    np.set_printoptions(precision=4)    # 设置计算精度，精确到小数点后四位

    '''各类样本均值向量'''
    mean_vector = []  # 每类的样本均值向量
    for i in y_c:
        mean_vector.append(np.mean(X[y == i], axis=0))  #求每个类别的样本的均值，axis=0是按列，axis=1是按行求
        #print('类别 %s 样本中心为: %s\n' % (i, mean_vector[i]))

    """类内离散度矩阵"""
    S_W = np.zeros((X.shape[1], X.shape[1]))
    for i in y_c:
        Xi = X[y == i] - mean_vector[i]     #此处将mean广播成(xx,4)
        S_W += np.mat(Xi).T * np.mat(Xi)
    #print('类内离散度矩阵:\n', S_W)

    """类间离散度矩阵"""
    S_B = np.zeros((X.shape[1], X.shape[1]))
    mu = np.mean(X, axis=0)  # 所有样本平均值
    for i in y_c:
        Ni = len(X[y == i])
        S_B += Ni * np.mat(mean_vector[i] - mu).T * np.mat(mean_vector[i] - mu)
    #print('类间离散度矩阵:\n', S_B)

    eigvalues, eigvectors = np.linalg.eig(np.linalg.inv(S_W) * S_B)  # 求特征值，特征向量
    #print('特征值:\n', eigvalues)
    #print('特征向量:\n', eigvectors)
    np.testing.assert_array_almost_equal(np.mat(np.linalg.inv(S_W) * S_B) * np.mat(eigvectors[:, 0].reshape(4, 1)), eigvalues[0] * np.mat(eigvectors[:, 0].reshape(4, 1)), decimal=6, err_msg='', verbose=True)      # 测试特征值和特征向量是否正确

    # 按特征值从大到小排序对应特征向量

    eig_pairs = [(np.abs(eigvalues[i]), eigvectors[:, i]) for i in range(len(eigvalues))]
    eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
    #print('pairs:\n', eig_pairs)
    W = np.hstack((eig_pairs[0][1].reshape(4, 1)))  # 降到一维
    W = W.reshape(4,1)
    #print('W*:\n', W)

    X_trans = X.dot(W)      # 降维后的X (120,1)
    assert X_trans.shape == (150-k_size, 1)
    #print('降维后的X:\n', X_trans)

    Xtest_trans = X_test.dot(W)
    #print('降维后的Xtest:\n', Xtest_trans)

    meanvector_trans = np.matrix(mean_vector).reshape(3,4).dot(W)
    #print('降维后的均值向量组:\n', meanvector_trans)

    accuracy = 0
    yindex = 0
    for vec in Xtest_trans:		# 分类时取离样本点最近的类均值中心的标签
        dist = [0,0,0]
        dist[0] = np.linalg.norm(vec - meanvector_trans[0])
        dist[1] = np.linalg.norm(vec - meanvector_trans[1])
        dist[2] = np.linalg.norm(vec - meanvector_trans[2])
        min_dist = min(dist)
        for i in range(3):
            if min_dist == dist[i]:
                if y_test[yindex] == i:
                    accuracy += 1
        yindex += 1
    ac1 += accuracy / k_size
    #print('准确率为:\n', accuracy/k_size)

    # 画图 x为训练集，o为测试集
    plt.figure(figsize=(5, 4))
    plt.scatter(X_trans[y == 0, 0], X_trans[y == 0, 0], c='r', marker='x')
    plt.scatter(X_trans[y == 1, 0], X_trans[y == 1, 0], c='g', marker='x')
    plt.scatter(X_trans[y == 2, 0], X_trans[y == 2, 0], c='b', marker='x')
    plt.title('Fisher')
    plt.xlabel('LD0')
    plt.ylabel('LD1')
    plt.legend(labels, loc='best', fancybox=True)
    plt.scatter(Xtest_trans[y_test == 0, 0], Xtest_trans[y_test == 0, 0], c='r', marker='.')
    plt.scatter(Xtest_trans[y_test == 1, 0], Xtest_trans[y_test == 1, 0], c='g', marker='.')
    plt.scatter(Xtest_trans[y_test == 2, 0], Xtest_trans[y_test == 2, 0], c='b', marker='.')
    plt.show()
    
print('一维平均正确率为:\n%.3f' % (ac1/10))

2. logistic 回归模型

• 首先是线性 logistic 回归，我们考虑一个二分类问题（不设置正则化）

# 导入需要用到的库
import numpy as np
import scipy.optimize as op
import matplotlib.pyplot as plt

def plot_data(X, y):		# 绘制所给数据的散点图
    plt.figure()
    pos = y==1		# 将绘制点分为两类
    neg = y==0
    plt.plot(X[pos,0],X[pos,1],'r+',markersize=5, markeredgewidth=1,label='Admitted')
    plt.plot(X[neg,0],X[neg,1],'bo',markersize=3, markeredgewidth=1,label='Not admitted')
    plt.xlabel("Exam 1 Score")
    plt.ylabel("Exam 2 Score")
    plt.legend(loc='upper right')
    
def plot_decision_boundary(theta, X, y):
    """绘制分类面"""
    plot_data(X[:, 1:], y)
    d = X.shape[1]
    if d <= 3:
        plot_x = np.array([np.min(X[:, 1])-2, np.max(X[:, 1])+2])
        plot_y = -1.0 / theta[2]*(theta[1]*plot_x + theta[0])
        plt.plot(plot_x, plot_y, 'm-', label="Decision Boundary")
        plt.xlim([30, 100])
        plt.ylim([30, 100])
    else:
        n_grid = 50
        u = np.linspace(-1, 1.5, n_grid)
        v = np.linspace(-1, 1.5, n_grid)
        z = np.zeros((n_grid, n_grid))
        for i in range(n_grid):
            for j in range(n_grid):
                uu, vv = np.array([u[i]]), np.array([v[j]])
                z[i, j] = np.dot(map_feature(uu, vv), theta)
        z = z.T
        CS = plt.contour(u, v, z, linewidths=2, levels=[0.0], colors=['m'])
        CS.collections[0].set_label('Decision boundary')
    plt.legend()
    
def sigmoid(z):
    """计算 sigmoid 函数"""

    z = np.asarray(z)
    g = np.zeros_like(z)
    g = 1/(1 + np.exp(-z))

    return g

def cost_function(theta, X, y):
    """无正则项的代价函数"""

    J = 0.0
    m = 1.0*len(y)
    one1 = np.ones_like(y)
    J = np.sum(-y.transpose() @ np.log(sigmoid(X @ theta)) - (one1 - y).transpose() @ np.log(one1 - sigmoid(X @ theta)))/m		# 交叉熵损失函数
    
    return J

def cost_gradient(theta, X, y):
    """无正则项的代价函数的梯度"""
    
    m = 1.0*len(y)
    grad = np.zeros_like(theta)
    grad = (sigmoid(X @ theta) - y).transpose() @ X/m
    
    return grad

def predict(theta, X):
    """预测函数，标签是 0 或 1"""
    m = X.shape[0]
    pred = np.zeros((m, 1), dtype=np.bool)
    ypre = X @ theta
    for i in range(len(sigmoid(ypre))):
        if sigmoid(ypre)[i] >= 0.5:
                pred[i][0] = True

    return pred

def logistic_regression():
    """建立逻辑回归模型"""

    # 加载数据
    data = np.loadtxt("LR_data1.txt", delimiter=",")
    X, y = data[:, :2], data[:, 2]

    # 计算代价与梯度
    m = X.shape[0]
    X = np.hstack((np.ones((m, 1)), X))

    # 初始化参数
    theta_initial = np.zeros_like(X[0])

    # 计算并打印初始参数对应的代价与梯度
    cost = cost_function(theta_initial, X, y)
    grad = cost_gradient(theta_initial, X, y)
    print("Cost at initial theta (zeros): ", cost)
    print("Gradient at initial theta (zeros): \n", grad)

    # 使用 scipy.optimize.fmin_cg 函数(梯度下降法) 优化模型参数
    args = (X, y)
    maxiter = 200
    ret = op.fmin_cg(cost_function,
                 theta_initial,
                 fprime=cost_gradient,
                 args=args,
                 maxiter=maxiter,
                 full_output=True)
    # 其中 cost_function 为代价函数，theta为需要优化的参数初始值，fprime=cost_gradient 给出了代价函数的梯度，args=(X, y) 给出了需要优化的函数与对应的梯度计算所需要的其他参数，maxiter=200 给出了最大迭代次数，full_output=True 则指明该函数除了输出优化得到的参数 theta_opt 外，还会返回最小的代价函数值 cost_min 等
    theta_opt, cost_min, _, _, _ = ret
    print("Cost at theta found by fmin_cg: ", cost_min)
    print("theta_op: \n", theta_opt)

    # 绘制分类面
    plot_decision_boundary(theta_opt, X, y)
    plt.show()

    # 预测 (45,80) 数据点的分类情况
    x_test = np.array([1, 45, 85.0])
    prob = sigmoid(np.dot(theta_opt, x_test))
    print('probability of class 1 is: ', prob)

    # 计算在训练集上的分类正确率
    p = predict(theta_opt, X)
    print("Train Accuracy: ", np.mean(p == y)*100.)
    
logistic_regression()

• 数据可视化的结果如下：
  

• 现在考虑一个非线性 logistic 回归，依然是二分类（考虑正则化）

• 其中特征变换函数解释如下：创建更多的特征是充分挖掘数据中的信息的一种有效手段。在函数 map_feature 中，我们将数据映射为其 6 阶多项式的所有项
  $$map\_feature(x)=\left[\begin{array}{c}1\\ x_1\\ x_2\\ x_1^2\\ x_1{x}_2\\ x_2^2\\ x_1^3\\ ⋮\\ x_1{x}_2^5\\ x_2^6\end{array}\right]$$

import numpy as np
import scipy.optimize as op
import matplotlib.pyplot as plt

def plot_data(X, y):		# 绘制所给数据的散点图
    plt.figure()
    pos = y==1		# 将绘制点分为两类
    neg = y==0
    plt.plot(X[pos,0],X[pos,1],'r+',markersize=5, markeredgewidth=1,label='Admitted')
    plt.plot(X[neg,0],X[neg,1],'bo',markersize=3, markeredgewidth=1,label='Not admitted')
    plt.xlabel("Exam 1 Score")
    plt.ylabel("Exam 2 Score")
    plt.legend(loc='upper right')
    
def plot_decision_boundary(theta, X, y):
    """绘制分类面"""
    plot_data(X[:, 1:], y)
    d = X.shape[1]
    if d <= 3:
        plot_x = np.array([np.min(X[:, 1])-2, np.max(X[:, 1])+2])
        plot_y = -1.0 / theta[2]*(theta[1]*plot_x + theta[0])
        plt.plot(plot_x, plot_y, 'm-', label="Decision Boundary")
        plt.xlim([30, 100])
        plt.ylim([30, 100])
    else:
        n_grid = 50
        u = np.linspace(-1, 1.5, n_grid)
        v = np.linspace(-1, 1.5, n_grid)
        z = np.zeros((n_grid, n_grid))
        for i in range(n_grid):
            for j in range(n_grid):
                uu, vv = np.array([u[i]]), np.array([v[j]])
                z[i, j] = np.dot(map_feature(uu, vv), theta)
        z = z.T
        CS = plt.contour(u, v, z, linewidths=2, levels=[0.0], colors=['m'])
        CS.collections[0].set_label('Decision boundary')
    plt.legend()
    
def sigmoid(z):
    """计算 sigmoid 函数"""

    z = np.asarray(z)
    g = np.zeros_like(z)
    g = 1/(1 + np.exp(-z))

    return g

def map_feature(X1, X2, degree=6):
    """特征变换函数，将低维非线性变为高维线性"""
    
    m = len(X1)
    assert len(X1) == len(X2)
    n = int((degree+2)*(degree+1)/2)
    out = np.zeros((m, n))
    idx = 0
    for i in range(degree+1):
        for j in range(i+1):
            out[:, idx] = np.power(X1, i-j)*np.power(X2, j)
            idx += 1

    return out

def cost_function_reg(theta, X, y, lmb):
    """带有正则项的代价函数"""
    
    m = 1.0*len(y)
    J = 0
    one1 = np.ones_like(y)
    J = np.sum(-y.transpose() @ np.log(sigmoid(X @ theta)) - (one1 - y).transpose() @ np.log(one1 - sigmoid(X @ theta)))/m + lmb/(2*m)*np.sum(theta*theta)

    return J

def cost_gradient_reg(theta, X, y, lmb):
    """带有正则项的代价函数的梯度"""

    m = 1.0*len(y)
    grad = np.zeros_like(theta)
    for i in range(len(grad)):
        condi = 0 if i ==0 else (lmb/m)*theta[i]
        grad[i]=(1/m)*np.sum((sigmoid(X @ theta)-y)*X[:,i])+condi

    return grad

def logistic_regression_reg(lmb=1.0):
    """建立逻辑回归模型"""

    # 加载数据
    data = np.loadtxt("LR_data2.txt", delimiter=",")
    X, y = data[:, :2], data[:, 2]

    # 计算具有正则项的代价与梯度
    # 注意map_feature会自动加入一列 1
    X = map_feature(X[:, 0], X[:, 1])
    
    # 初始化参数
    theta_initial = np.zeros_like(X[0, :])

    # 计算并打印初始参数对应的代价与梯度
    cost = cost_function_reg(theta_initial, X, y, lmb=lmb)
    grad = cost_gradient_reg(theta_initial, X, y, lmb=lmb)
    print("Cost at initial theta (zeros): ", cost)
    print("Gradient at initial theta (zeros): \n", grad)

    # 使用 scipy.optimize.fmin_cg 优化模型参数
    args = (X, y, lmb)
    maxiter = 200
    ret = op.fmin_cg(cost_function_reg,
                 theta_initial,
                 fprime=cost_gradient_reg,
                 args=args,
                 maxiter=maxiter,
                 full_output=True)
    theta_opt, cost_min, _, _, _ = ret
    print("Cost at theta found by fmin_cg: ", cost_min)
    print("theta_op: \n", theta_opt)

    # 绘制分类面
    plot_decision_boundary(theta_opt, X, y)
    plt.title("lambda = " + str(lmb))
    plt.show()

    # 计算在训练集上的分类正确率
    pred = predict(theta_opt, X)
    print("Train Accuracy: ", np.mean(pred == y)*100)

logistic_regression_reg(lmb=1.0) 

• 数据可视化结果如下：
  

数据集下载：

• 链接：https://pan.baidu.com/s/1kxqdtHtseJbk29QNcDgiAA?pwd=wm1r
  os): \n", grad)