Introduction

这是一本概率统计进阶版本的知识点理解书籍，相关内容基于作者个人学习工作中的理解整理。各位读者如有相关问题和建议可以在相应的github库中以issue形式提交，作者欢迎各种批评指正，但鉴于个人时间有限，目前更新周期不定，请谅解。作者个人详情可登录作者的个人网站了解

Contributing

本电子书可以在作者的个人github上找到相关源码，社区维护依赖于issue，具体修正由作者本人管理操作。

License

本书开源代码和文档使用GPLv3协议，商业用途必须联系作者，否则自负相关法律责任。Copyright by 施华。

NO.1.tip: 条件期望

背景

  许多故事的开始似乎总有那样的身影让我们似曾相识......

  1800年代Francis Galton第一次用回归一词命名了一种朝向均值的现象模型，计量经济学中为了分析条件概率也常常使用其对应的一阶矩进行建模分析，这些数学模型在概率统计中拥有着一个共同的名字——条件期望。从经典的回归到计量经济学，从现在的深度学习到强化学习，它一直出现在每一个故事的许多角落，或开始、或中章、或结尾。下面让我们一起缓缓揭开条件期望的面纱。

源与流

  首先，我们使用测度论中的Radon Nikodym定理给出条件期望的严格定义。设$(X,\mathcal{F},P)$是一个概率空间，$f$是它上面积分存在的随机变量。又设$\mathcal{G}$是$\mathcal{F}$的子$\sigma$-域，即$\mathcal{G}$是一个$X$上的$\sigma$域，而且$\mathcal{G}\subset \mathcal{F}$。对每个$A\in \mathcal{G}$，令$$\phi (A)={\int }_{A}fdp$$ 易见$\phi$是$\mathcal{G}$上的符号测度，它和限制在$\mathcal{G}$上的测度$P$之间有关系$\phi \ll P$。因此，由Radon Nikodym定理知，存在$(X,\mathcal{F},P)$上$a.s.$意义下唯一的可测函数$E(f|\mathcal{G})$，使对每个$A\in \mathcal{F}$有$$\phi (A)={\int }_{A}E(f|\mathcal{G})dp$$(这个可测函数之所以记为$E(f|\mathcal{G})$，是因为它既与$f$有关，又与$\mathcal{G}$有关)。利用$E(f|\mathcal{G})$的性质，可以把条件期望公理化地定义如下：

  设$f$为概率空间$(X,\mathcal{F},P)$上积分存在的随机变量。称$E(f|\mathcal{G})$为$f$关于$\mathcal{F}$的子$\sigma$域$\mathcal{G}$的条件期望，如果
  (1) $E(f|\mathcal{G})$是$(X,\mathcal{F},P)$上积分存在的可测函数；
  (2) 对任何$A\in \mathcal{G}$，有$${\int }_{A}E(f|\mathcal{G})dp={\int }_{A}fdp$$

  从通俗意义上理解，条件期望即是一种符号测度对一种测度的Radon Nikodym导数，而严格的Radon Nikodym导数定义如下：

  设$\phi$是测度空间$(X,\mathcal{F},\mu )$上的符号测度。如果存在$a.e.$意义下唯一的可测函数使下式成立，则称$f$为$\phi$对$\mu$的Radon Nikodym导数$$\rho (A)={\int }_{A}fd\mu ,\forall A\in \mathcal{F}$$

  从几何意义上理解，条件期望是$X$到$L^2$空间的一种正交投影算子，因此它也会具有正交投影的特殊性质：
  (1) 正交投影算子可以有效地将向量从一个低维空间映射到另一个低维空间，而不会破坏其所表示的特征信息；
  (2) 正交投影算子的作用是可逆的，双向投影结果完全一致。

  条件期望正因它具有如此良好的性质，才让我们在许多场景中使用它来建模现实。那么，接下来问题来了，我们如何明确一个条件期望呢？换一个角度思考条件期望，比如$E(X|Y)$，在观察$Y$之前，并不知道$E(X|Y=y)$的值，所以条件期望本质上是一种随机变量。循着这个思路，我们知道，为了确定一个随机变量，我们往往需要对其概率分布建模，而在统计学中，我们通常会假设样本来自一个给定的$f(X|\Theta )$。我们将这样的函数称为概率密度函数，然后采用各种方法去估计参数$\Theta$。由此引出一个关键概念——估计量，对于这个估计量，有多种方法求得。常见的方法有矩估计、极大似然估计、贝叶斯估计和EM算法。

  在此，故事已开始，精彩待纷呈。后续，我们会开始进入到估计量的四大估计方法的传奇旅程。

NO.2.tip: 矩估计

背景

  矩法也许是最早的求点估计量的方法，至少可以追溯到19世纪末的KarlPearson。此法的优点是使用很简单，从而几乎总是可以求出估计值。尽管令人遗憾的是在很多情况下，矩法导出的估计量还需要改进，但是在其他方法难以实施的时候，它仍然不失为一个很好的工作起点。另外，它也可以作为其他需要循环几次的算法的初始值。

源与流

  设$X_1,...,X_n$是来$f(x|{\theta }_1,...,{\theta }_n)$为其概率密度函数的总体的样本。矩法估计量这样得到的：令前$k$阶的样本句子矩与相应的前$k$阶总体矩相等，这样就得到一个联立方程组，求解之，就得到矩统计量。更清楚地说，我们定义 $$\{\begin{array}{l}{m}_1=\frac{1}{n}{\sum }_{i=1}^n{X}_i^{{}^{\prime }},{\mu }_1^{{}^{\prime }}=E{X}^{{}^{\prime }},\\ m_2=\frac{1}{n}{\sum }_{i=1}^n{X}_i^2,{\mu }_2^{{}^{\prime }}=E{X}^2,\\ ...\\ m_k=\frac{1}{n}{\sum }_{i=1}^n{X}_i^k,{\mu }_k^{{}^{\prime }}=E{X}^k,\end{array}$$

  在典型的情况下，总体矩${\mu }_j^{{}^{\prime }}$是参数${\theta }_1,...,{\theta }_k$的一个函数，可以记作${\mu }_j^{{}^{\prime }}$。于是$({\theta }_1,...,{\theta }_k)$的矩法估计量$({\theta }_1,...,{\theta }_k)$就可以通过求解下面的关于$({\theta }_1,...,{\theta }_k)$的方程组 $$\{\begin{array}{l}{m}_1={\mu }_1^{{}^{\prime }}({\theta }_1,...,{\theta }_k),\\ m_2={\mu }_2^{{}^{\prime }}({\theta }_1,...,{\theta }_k),\\ ...\\ m_k={\mu }_k^{{}^{\prime }}({\theta }_1,...,{\theta }_k),\end{array}$$ 得到。

  值得一提的是，在计量经济学的世界中，矩法已是大名鼎鼎的明星。它是工具变量法（Instrumental Variable,IV）和广义矩估计法（Generalized Moment Method,GMM）的基础。在矩估计法中关键是利用了由随机干扰项的条件均值零假设所推出的非条件零均值特性，以及随机干扰项各解释变量间同期不相关特性$$E(X_i^{{}^{\prime }}{\mu }_i)=0$$作为总体矩条件。如果某个解释变量与随机干扰项相关，只要能找到1个工具变量，仍然可以构成一组矩条件，这就是工具变量法。如果存在大于$k+1$个变量与随机干扰项不相关，可以构成一组包含大于$k+1$个方程的矩条件，这就是广义矩估计法。

NO.3.tip: 最大似然估计

背景

  最大似然估计是统计中的一个重要理论框架，作为一门优化技术，他用最朴素的表达--凸优化--描绘了机器学习的万里长卷

源与流

  首先，我们了解一下什么是似然函数。考虑$R^m$上的一族概率分布，每个概率分布对应一个向量$x\in R^n$，概率密度为$p_x(.)$。固定$y\in R^m$，$p_x(y)$可以看成$x$的函数，成为似然函数。事实上，考虑似然函数的对数更为方便，我们称为对数-似然函数，用$l$表示，即$$l(x)=log{p}_x(y)$$ 通常对参数$x$都有一些曰肃，可以作为x的先验知识，或者是似然函数的定义域。

  现在考虑如下问题，根据观测到的服从分布的一个样本$y$，估计参数$x$的值。最大似然估计就是一种广为采用的方法，其中参数$x$的估计为 $${\hat{x}}_{ml}=argma{x}_x{p}_x(y)=argma{x}_{x}l(x)$$ 即选择使得似然（或者对数似然）函数在$y$的观测值处最大的那个参数值作为$x$的估计只。如果我们有关于$x$的先验信息，比如说$x\in C\subseteq R^n$，我们可以显示添加约束$xinC$，或者隐式添加约束，当$x\notin C$时，定义$p_x(y)$为$D$。

  求取参数向量$x$的最大似然估计问题可以表述如下： $$\begin{array}{r}maximizel(x)=log{p}_x(y)\\ subjecttox\in C\end{array}$$ 其中$x\in C$给出了参数向量$x$的先验信息或者其他约束。在这个优化问题中，向量$x\in R^n$（在概率密度中是参数）是优化变量，向量$y\in R^m$（观测到样本）是问题参数。

  对于$y$的每个观测值，如果对数-似然函数$l$都是凹的，且集合$C$可以表述为线性等式约束以及凸不等式约束的组合，那么最大似然估计问题是凸优化问题。事实上，很多这一类问题都满足这个条件。对这些问题，可以通过凸优化确定最大似然估计。

  现在，我们已经清楚了最大似然估计的优化表达，如何应用到具体算法场景中呢？不要着急，我们需要加入一个叫$IID$（独立同分布）的假设，才能让足底啊似然估计真正落地。

  考虑线性测量模型 $$y_i{\alpha }_i^T+v_i,i=1,...,m$$ 其中，$x\in R^m$是待估计参数向量，$y_i\in R$是测量量或观测量，$v_i$是测量误差或者噪声。假设$v_i$独立同分布$(IID)$，在$R$上具有概率密度$p$，此时似然函数为 $$p_x(y)=\prod _{i=-1}^{m}p(y_i-{\alpha }_{i^T}x)$$ 因此对数-似然函数为 $$l(x)=log{p}_x(y)=\sum _i^{m}logp(y_i-{\alpha }_{i^T}x)$$ 最大似然估计是下面优化问题的任意一个最优点： $$maximize\sum _i^{m}logp(y_i-{\alpha }_{i^T}x)$$ 其中优化变量为$x$。如果概率密度$p$是对数-凹的，此优化问题是凸优化问题。

  在有了$IID$假设后，我们还不能去优化似然函数，因为其中的概率分布$p$还没有明确的解析表达。接下来，我们通过为噪声指定具体分布的方式来确定似然函数的明确解析表达，进而凸优化问题可解。常见的分布有以下两种：

  Gauss噪声：当$v_i$是Gauss噪声，均值为$0$，方差为${\sigma }^2$时，概率密度函数为$p(z)=(2\pi {\sigma }^2{)}^{-\frac{1}{2}}{e}^{-\frac{z^2}{2{\sigma }^2}}$，则对数-似然函数为： $$l(x)=-(\frac{m}{2})log(2\pi {\sigma }^2)-\frac{1}{2{\sigma }^2}\Vert Ax-y{\Vert }_2^2$$ 其中矩阵$A$的行向量为${\alpha }_1^T,...,{\alpha }_m^T$。因此$x$的最大似然估计为$x_{ml}=argma{x}_x\Vert Ax-y{\Vert }_2^2$，也即$L_2-$范数逼近问题的解。

  Laplace噪声：当$v_i$服从Laplace分布时，即具有概率密度$p(z)=(\frac{1}{2a})e^{-\frac{|z|}{a}}$，其中$a>0$。$x$的最大似然估计为$\hat{x}=argma{x}_x\Vert Ax-y\Vert$，也即$L_1-$范数逼近问题的解。

  在这里，我们也可以从最大似然估计的角度理解$L_1-$范数逼近在大误差情况下的鲁棒性。可以将$L_1-$范数逼近问题看成是概率密度函数是Laplace函数的最大似然估计；$L_2-$范数逼近问题可以看成是噪声服从Gauss分布的最大似然估计。Laplace密度函数比Gauss密度函数具有更大的截尾，即对于Laplace密度函数，$v_i$取值很大的概率远远大于Gauss密度函数的情况。因此，相应的最大似然估计方法将会接收更多大残差。

NO.4.tip: 贝叶斯估计

背景

  很多情况下，我们都会遇到预测方差过高问题，它的一个常见表现就是过拟合。想要缓解甚至处理方差过高问题，我们要了解的是参数不确定性。以下是参数不确定性的定义以及它对预测方差的影响：

  假设模型参数化线性模型 $$\hat{f}(x;\theta )=\varphi (x{)}^T\theta$$ $\varphi (x)$是基函数（高次多项式、神经网络特征等）$\theta$待估参数。对固定$x_0$，预测 $$\hat{f}(x_0;D)=\varphi (x_0{)}^T\hat{\theta }(D)$$ $\hat{\theta }(D)$是由随机数据集$D$得到的参数估计，本身是随机变量，参数不确定性即为 $$Var(\hat{\theta })$$

  有了参数不确定性的定义，我们就可以推导其对预测方差的影响： $$Va{r}_D[\hat{f}(x_0)]=Var[\varphi (x_0{)}^T\hat{\theta }]={\varphi }_{x_0}^{T}cov(\hat{\theta })\varphi (x_0)$$ $cov(\hat{\theta })$是参数协方差矩阵，对角线为各参数方差，刻画了参数整体不确定性。我们发现，预测方差是参数协方差矩阵的二次型：

（1）只要参数存在不确定性（$cov(\hat{\theta })$非零，特征值大），预测方差一定上升

（2）参数越不确定（$cov(\hat{\theta })$越大），同等特征下预测波动越强，模型方差越高

  注意，$cov(\hat{\theta })$是协方差矩阵，对矩阵的不确定性度量一般从以下四个角度入手：

（1）迹，所有参数自身方差想家，仅累加单参数波动，忽略参数间相关程度

（2）行列式，表征参数联合搞死分布置信椭球的体积，数值越大代表整体联合不确定性越高

（3）最大特征值（范数），捕捉参数空间波动最剧烈的单一方向，用来识别不稳定主成分

（4）Frobenious范数，整合全部方差与两两参数协方差，完整度量协方差矩阵整体波动能量

源与流

  现在，我们知道了饿预测方差过高很大程度上来源于参数不确定性，那么我们就可以从这里入手。实际上，我们有两种思路：

（1）一种是缩小参数估计的范围，从而间接降低参数不确定性；

（2）一种是对参数不确定性完整建模，从而直接控制参数不确定性，这种思想源自于统计学的贝叶斯理论，是一种标准贝叶斯方法

  首先，我们来看下第一种，缩小参数估计范围的思路。这种思想的两个典型例子就是Lasso回归和Ridge回归。

Ridge回归： $$\begin{array}{r}\widehat{{\beta }_{ridge}}=argmi{n}_{\beta }\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^P{x}_{ij}{\beta }_j{)}^2\\ subjectto:\sum _{j=1}^P{\beta }_j^2\le t\end{array}$$ Lasso回归： $$\begin{array}{r}\widehat{{\beta }_{lasso}}=argmi{n}_{\beta }\sum _{i=1}^N(y_i-{\beta }_0-\sum _{j=1}^P{x}_{ij}{\beta }_j{)}^2\\ subjectto:\sum _{j=1}^P|{\beta }_j|\le t\end{array}$$

  接下来，我们再开看第二种思路，对参数不确定性完整建模--标准贝叶斯方法。我们先给出贝叶斯公式： $$p(\theta ｜y)=\frac{p(y|\theta )p(\theta )}{p(y)}$$ 其中$p(\theta )$是关联的先验概率密度函数。我们对于以上狮子不再仅仅是要使分子最大，还要将$P(\theta ｜y)$作为一个整体来使用。这里，秘密大部分在于分母$p(y)$，它基本上是归一化常数： $$p(y)=\int p(y|\theta )p(\theta )d\theta$$

  正如我们看到的，在$p(y)$中隐藏的信息远远超过了知识为了计算$p(\theta ｜y)$的需要。上式的计算难点在于，一般来说，积分的估计没有解析计算方法。在这种情况下，人们需要借助近似技术来获得所需的信息。为了这个目的，可以使用很多种方法，更具体地说，我们将考虑以下这些方法：

（1）拉普拉斯近似法

（2）变分近似法

（3）蒙特卡洛积分求值技术

NO.5.tip： 变分推理

背景

  数学上，变分法就是求使积分泛函取极值的未知函数，相当于函数版的求导找最值。统计上，考虑一个具有未知（潜在）变量 $z$、已知变量 $x$ 和固定参数 $\theta$ 的模型。（如果参数未知，可以将这些参数添加到 $z$ 中，稍后我们将展开讨论。）我们假设先验为 $p_{\theta }(z)$，似然为 $p_{\theta }(x|z)$，因此未归一化联合概率为 $p_{\theta }(x,z)=p_{\theta }(x|z)p_{\theta }(z)$，后验为：

$$p_{\theta }(z|x)=\frac{p_{\theta }(x,z)}{p_{\theta }(x)}$$

  由于我们在函数（即分布 $q$）上最小化，所以这被称为变分方法。在实践中，我们选择一个参数族 $\mathcal{Q}$，在此处我们使用 $\psi$，称为变分参数，来指定我们使用的是族中的哪一个成员。对于给定的 $x$，我们可以使用以下公式来计算最佳变分参数：

$$\begin{array}{rl}{\psi }^{\ast }& =\arg \underset{\psi }{\min }{D}_{KL}(q_{\psi }(z)\Vert p_{\theta }(z|x))\\ & =\arg \underset{\psi }{\min }{\mathbb{E}}_{q_{\psi }(z)}\left[\log q_{\psi }(z)-\log \left(\frac{p_{\theta }(x|z)p_{\theta }(z)}{p_{\theta }(x)}\right)\right]\\ & =\arg \underset{\psi }{\min }{\mathbb{E}}_{q_{\psi }(z)}\left[\log q_{\psi }(z)-\log p_{\theta }(x|z)-\log p_{\theta }(z)\right]+\log p_{\theta }(x)\\ & =\arg \underset{\psi }{\min }{\mathbb{E}}_{q_{\psi }(z)}\left[\log q_{\psi }(z)-\log p_{\theta }(x|z)-\log p_{\theta }(z)\right]\end{array}$$

  一般而言，最后一项 $\log p_{\theta }(x)=\int p_{\theta }(x,z)dz$ 难以计算。幸运的是，该项独立于 $\psi$，因此我们可以忽略该项。结果，我们留下了第一个项，可以记作：

$$\mathcal{L}(\theta ,\psi |x)={\mathbb{E}}_{q_{\psi }(z)}\left[-\log p_{\theta }(x,z)+\log q_{\psi }(z)\right]$$

  最小化这个目标将最小化 KL 散度，使我们的近似接近真正的后验。

源与流

  因为KL散度非负，则

$$\begin{array}{rl}{D}_{KL}(q_{\psi }(z)\Vert p_{\theta }(z|x))& =\mathcal{L}(\theta ,\psi |x)+\log p_{\theta }(x)\ge 0\\ \mathcal{L}(\theta ,\psi |x)& \ge -\log p_{\theta }(x)\end{array}$$

第一项的负值被称为证据下界函数：

$$\mathcal{L}(\theta ,\psi |x)\triangleq {\mathbb{E}}_{q_{\psi }(z)}[\log p_{\theta }(x,z)]-\log q_{\psi }(z)]=\text{证据下界}$$

  "证据下界"这个名字的产生是因为：

$$\mathcal{L}(\theta ,\psi |x)\le \log p_{\theta }(x)$$

  其中，$\log p_{\theta }(x)$ 被称为"证据"。因此，最大化相对于 $\psi$ 的证据下界将最小化原始 KL 散度，因为 $\log p_{\theta }(x)$ 相对于 $\psi$ 是一个常数。（请注意：我们使用符号 $\mathcal{L}$ 来表示证据下界，而不是 $\mathcal{L}$，因为我们希望最大化 $\mathcal{L}$，但最小化 $\mathcal{L}$。）

  我们可以将证据下界重写如下：

$$\mathcal{L}(\theta ,\psi |x)={\mathbb{E}}_{q_{\psi }(z)}[\log p_{\theta }(x,z)]+\mathbb{H}(q_{\psi }(z))$$

  我们可以将其解释为：

证据下界 = 期望对数联合概率 + 熵

  第二项鼓励后验具有最大熵，而第一项鼓励后验成为联合最大后验估计配置。

  我们也可以将证据下界重写为：

$$\begin{array}{rl}\mathcal{L}(\theta ,\psi |x)& ={\mathbb{E}}_{q_{\psi }(z)}[\log p_{\theta }(x|z)+\log p_{\theta }(z)-\log q_{\psi }(z)]\\ & ={\mathbb{E}}_{q_{\psi }(z)}[\log p_{\theta }(x|z)]-D_{KL}(q_{\psi }(z)\Vert p_{\theta }(z))\end{array}$$

  我们可以将其解释为：

证据下界 = 期望对数似然 - 从后验到先验的 KL 散度

  KL 散度项的作用就像一个正则化因子，用于防止后验与先验偏离过多。

  存在有两种主要的方法来选择变分后验 $q_{\psi }(z|x)$ 的形式。在第一种方法中，我们选择一种方便的函数形式，例如多元高斯分布，然后使用基于梯度的方法优化证据下界。这被称为固定形式变分推理。另一种方法是进行均值场假设，即后验因子分解：

$$q_{\psi }(z)=\prod _{j=1}^J{q}_{{\psi }_j}(z_j)$$

  其中，$q_j(z_j)=q_{{\psi }_j}(z_j)$ 是第 $j$ 组变量的后验。我们不需要为每个 $q_j$ 指定函数形式。相反，通过以坐标上升的方式，一次一个地最大化相对于每组变分参数的证据下界，来推导出最优分布形式。因此，这被称为自由形式变分推理。

  下界最大化（即标准变分推理）存在的一个问题是，我们正在最小化 $D_{KL}(q\Vert p)$，这会导致迫零性行为，这意味着 $q(z|x)$ 往往过于紧凑（过于自信），以避免 $q(z|x)>0$ 但 $p(z|x)=0$ 的情况，这将导致无限的 KL 惩罚。

  尽管迫零性对于一些多模态后验（例如，混合模型）来说可能是理想的行为，但对于许多单模态后验来说（例如，贝叶斯逻辑回归或者具有对数凹似然的高斯过程）就不那么合理了。避免这个问题的一种方法是最小化 $D_{KL}(p\Vert q)$，这是避零性，这往往会导致广泛的后验分布，从而避免过度自信。在本节中，我们讨论期望传播。期望传播可以被视为 $D_{KL}(p\Vert q)$ 的局部近似。

NO.6.tip： 马尔可夫链蒙特卡罗方法

背景

  我们经常需要计算随机变量的某个函数的期望值 $\mathbb{E}[f(\mathbf{X})]$。这需要计算以下积分：

$$\mathbb{E}[f(\mathbf{x})]=\int f(\mathbf{x})p(\mathbf{x})\mathrm{d}\mathbf{x}$$

  其中，$\mathbf{x}\in {\mathbb{R}}^n$，$f:{\mathbb{R}}^n\to {\mathbb{R}}^m$，$p(\mathbf{x})$ 是 $\mathbf{X}$ 的目标分布。在低维度（例如3维）中，我们可以使用数值积分有效地计算上述积分。数值积分（自适应地）计算网格，然后评估网格上每个点的函数。但这无法扩展到更高的维度。

  另一种方法是抽取多个随机样本 ${\mathbf{x}}_n\sim p(\mathbf{x})$，然后计算：

$$\mathbb{E}[f(\mathbf{x})]\approx \frac{1}{N_s}\sum _{n=1}^{N_s}f({\mathbf{x}}_n)$$

  这被称为蒙特卡罗积分。与数值积分相比，蒙特卡罗积分的优势在于，函数只在概率不可忽略的地方进行评估，因此不需要均匀覆盖整个空间。特别是，可以证明，精度原则上与 $\mathbf{x}$ 的维度无关，仅取决于样本数量 $N_s$。关键在于，我们首先需要一种生成样本 ${\mathbf{x}}_n\sim p(\mathbf{x})$ 的方法。

  有几种常见的确定性的采样方式：

逆概率变换：用随机数"反查"概率分布的累积函数，把均匀随机变量变成服从目标分布的样本，从而用样本均值近似任意复杂函数的积分。

拒绝采样：在目标分布难以直接采样时，用一个容易采样的"建议分布"生成候选点，然后按目标分布与建议分布的比值概率接受或拒绝这些点，最终保留下来的样本就服从目标分布，进而用于蒙特卡洛积分。

重要性采样：不直接从目标分布采样，而是从一个更容易采样的"建议分布"中抽取样本，然后给每个样本赋予一个权重（目标密度与建议密度之比）来修正偏差，让积分估计更集中在函数值大的区域，从而降低方差、提高效率。

  确定性方法用规则网格铺满空间，误差 $O(N^{-k/d})$ 随维度指数恶化；蒙特卡洛用随机样本估计期望，误差 $O(N^{-1/2})$ 与维度无关，只取决于样本数N——前者靠"覆盖"换精度，后者靠"概率"换自由。

源与流

  因此，马尔可夫链蒙特卡罗方法就是这样一种流行的从高维分布中采样的方法。马尔可夫链蒙特卡罗方法背后的基本思想是在状态空间 $\mathcal{X}$ 上构建一个马尔可夫链，其平稳分布是感兴趣的目标密度 $p^{\ast }(x)$。（在贝叶斯上下文中，这通常是一个后验，$p^{\ast }(x)\propto p(x\mid \mathcal{D})$，但马尔可夫链蒙特卡罗方法可以应用于从任何类型的分布生成样本。）也就是说，通过从链中抽取（相关的）样本 $x_0,x_1,x_2,\cdots$，我们可以相对于 $p^{\ast }$ 进行蒙特卡罗积分。

  Metropolis-Hastings 算法的基本思想是，在每一步中，我们都建议从当前状态 $x$ 移动到概率为 $q(x^{\prime }\mid x)$ 的新状态 $x^{\prime }$，其中 $q$ 被称为提议分布（也称为核）。用户可以自由使用他们想要的任何类型的提议分布，但须遵守我们在下面解释的一些条件。这使得 Metropolis-Hastings 算法成为一种非常灵活的方法。

  在提议移动到 $x^{\prime }$ 之后，我们根据某种公式决定是接受还是拒绝这一提议，该公式确保在每个状态花费的长期时间与 $p^{\ast }(x)$ 成比例。如果提议被接受，则新状态为 $x^{\prime }$，否则新状态与当前状态 $x$ 相同（即重复样本）。

  如果提议是对称的，那么 $q(x^{\prime }\mid x)=q(x\mid x^{\prime })$，接受概率由以下公式给出：

$$A=\min \left(1,\frac{p^{\ast }(x^{\prime })}{p^{\ast }(x)}\right)$$

  我们看到，如果 $x^{\prime }$ 比 $x$ 更有可能，那么我们肯定会移动到 $x^{\prime }$（因为 $\dfrac{p^{}(x')}{p^{}(x)} > 1$），但如果 $x^{\prime }$ 不太可能，我们仍然可能移动到 $x^{\prime }$，这取决于相对概率。因此，我们不是贪婪地只进入更可能的状态，而是偶尔允许"下坡"进入不太可能的状态。可以证明，该过程可以确保在每个状态 $x$ 下花费的时间等于 $p^{\ast }(x)$。

  如果提议分布是不对称的，那么 $q(x^{\prime }\mid x)\ne q(x\mid x^{\prime })$，我们需要 Hastings 修正（Hastings correction），由以下公式给出：

$$\begin{array}{rl}A& =\min (1,\alpha )\\ \alpha & =\frac{p^{\ast }(x^{\prime })q(x\mid x^{\prime })}{p^{\ast }(x)q(x^{\prime }\mid x)}=\frac{p^{\ast }(x^{\prime })/q(x^{\prime }\mid x)}{p^{\ast }(x)/q(x\mid x^{\prime })}\end{array}$$

  这种修正是为了补偿提议分布本身（而不仅仅是目标分布）可能有利于某些状态的事实。

  Metropolis-Hastings 算法之所以是一种有用算法的一个重要原因是，在评估 $\alpha$ 时，我们只需要知道达到归一化常数的目标密度。特别地，假设 $p^{\ast }(x)=\frac{1}{Z}\tilde{p}(x)$，其中 $\tilde{p}(x)$ 是未归一化分布，$Z$ 是归一化常数。于是：

$$\alpha =\frac{(\tilde{p}(x^{\prime })/Z)q(x\mid x^{\prime })}{(\tilde{p}(x)/Z)q(x^{\prime }\mid x)}$$

  所以 $Z$ 可以消除。因此，即使 $Z$ 未知，我们也可以从 $p^{\ast }$ 中采样。

  Metropolis-Hastings 方法的主要问题是需要选择提议分布，以及接受率可能较低。还有一种经典方法吉布斯采样，该方法利用图模型的条件独立性来自动创建一个好的提议分布，接受概率为1。

NO.7.tip: 指数族分布

背景

  指数分布在统计结构中占用重要地位，他得成员具有许多共同的重要性质，对这些性质的一般性讨论可以带来不少启发。

  给定参数 $\mathbf{\eta }$，关于 $\mathbf{x}$ 的指数族分布可以定义为以下形式的分布的集合：

$$p(\mathbf{x}\mid \mathbf{\eta })=h(\mathbf{x})g(\mathbf{\eta })\exp \left\{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\text{\hspace{1em}(1)}$$

其中 $\mathbf{x}$ 可以是标量或矢量，它可以是离散的或连续的。$\mathbf{\eta }$ 则称为分布的自然参数，$\mathbf{u}(\mathbf{x})$ 是 $\mathbf{x}$ 的某个函数。函数 $g(\mathbf{\eta })$ 可以解释为用来确保分布被归一化的系数，因此满足

$$g(\mathbf{\eta })\int h(\mathbf{x})\exp \left\{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{u}(\mathbf{x})\right\}\mathrm{d}\mathbf{x}=1\text{\hspace{1em}(2)}$$

其中 $\mathbf{x}$ 若为离散变量，则积分改为求和。

源与流

  首先考虑伯努利分布：

$$p(x\mid \mu )=\mathrm{Bern}(x\mid \mu )={\mu }^x(1-\mu {)}^{1-x}\text{\hspace{1em}(3)}$$

  将式 (3) 的右侧表示为对数的指数，可得

$$\begin{array}{rl}p(x\mid \mu )& =\exp \{x\ln \mu +(1-x)\ln (1-\mu )\}\\ & =(1-\mu )\exp \left\{\ln \left(\frac{\mu }{1-\mu }\right)x\right\}\end{array}\text{\hspace{1em}(4)}$$

  与式 (1) 比较后，我们可以确定

$$\eta =\ln \left(\frac{\mu }{1-\mu }\right)\text{\hspace{1em}(5)}$$

  对 $\mu$ 进行求解可得 $\mu =\sigma (\eta )$，其中

$$\sigma (\eta )=\frac{1}{1+\exp (-\eta )}\text{\hspace{1em}(6)}$$

  它又称为sigmoid 函数。因此，我们可以使用以下标准表达式来描述伯努利分布：

$$p(x\mid \eta )=\sigma (-\eta )\exp (\eta x)\text{\hspace{1em}(7)}$$

  这里使用了等式 $1-\sigma (\eta )=\sigma (-\eta )$，其易由式 (6) 证得。与式 (1) 比较可得

$$u(x)=x\text{\hspace{1em}(8)}$$

$$h(x)=1\text{\hspace{1em}(9)}$$

$$g(\eta )=\sigma (-\eta )\text{\hspace{1em}(10)}$$

  接下来考虑多项分布，对于单个观测 $\mathbf{x}$，其形式为

$$p(\mathbf{x}\mid \mathbf{\mu })=\prod _{k=1}^M{\mu }_k^{x_k}=\exp \left\{\sum _{k=1}^M{x}_k\ln {\mu }_k\right\}\text{\hspace{1em}(11)}$$

  其中 $\mathbf{x}=(x_1,\cdots ,x_M{)}^{\mathrm{T}}$。式 (11) 同样可以写成标准表达式：

$$p(\mathbf{x}\mid \mathbf{\eta })=\exp ({\mathbf{\eta }}^{\mathrm{T}}\mathbf{x})\text{\hspace{1em}(12)}$$

  其中 ${\eta }_k=\ln {\mu }_k$，并且我们定义 $\mathbf{\eta }=({\eta }_1,\cdots ,{\eta }_M{)}^{\mathrm{T}}$。与式 (1) 比较可得

$$\mathbf{u}(\mathbf{x})=\mathbf{x}\text{\hspace{1em}(13)}$$

$$h(\mathbf{x})=1\text{\hspace{1em}(14)}$$

$$g(\mathbf{\eta })=1\text{\hspace{1em}(15)}$$

  注意，参数 ${\eta }_k$ 不是独立的，因为参数 ${\mu }_k$ 受到下式的约束：

$$\sum _{k=1}^M{\mu }_k=1\text{\hspace{1em}(16)}$$

  因此，对于给定的任意 $M-1$ 个参数 ${\mu }_k$，剩余参数的值被固定。在某些情况下，通过仅用 $M-1$ 个参数表达分布，可以方便地去除这一约束。这可以通过使用关系式 (16) 将 ${\mu }_M$ 用剩余的 $\{{\mu }_k\}$（其中 $k=1,\cdots ,M-1$）表示，从而留下 $M-1$ 个参数来实现。请注意，这些剩余参数仍然受到下式的约束：

$$0⩽{\mu }_k⩽1,\sum _{k=1}^{M-1}{\mu }_k⩽1\text{\hspace{1em}(17)}$$

  利用约束式 (16)，多项分布变为

$$\begin{array}{rl}& \exp \left\{\sum _{k=1}^M{x}_k\ln {\mu }_k\right\}\\ & =\exp \left\{\sum _{k=1}^{M-1}{x}_k\ln {\mu }_k+\left(1-\sum _{k=1}^{M-1}{x}_k\right)\ln \left(1-\sum _{k=1}^{M-1}{\mu }_k\right)\right\}\\ & =\exp \left\{\sum _{k=1}^{M-1}{x}_k\ln \left(\frac{{\mu }_k}{1-{\sum }_{j=1}^{M-1}{\mu }_j}\right)+\ln \left(1-\sum _{k=1}^{M-1}{\mu }_k\right)\right\}\end{array}\text{\hspace{1em}(18)}$$

  现在我们确定了

$$\ln \left(\frac{{\mu }_k}{1-{\sum }_j{\mu }_j}\right)={\eta }_k\text{\hspace{1em}(19)}$$

  可以通过首先对式 (19) 两边的 $k$ 求和，然后重新排列和反向代入来求解 $k$：

$${\mu }_k=\frac{\exp ({\eta }_k)}{1+{\sum }_j\exp ({\eta }_j)}\text{\hspace{1em}(20)}$$

  它又称为softmax 函数或归一化指数。因此在这种表示中，多项分布将采用以下形式：

$$p(\mathbf{x}\mid \mathbf{\eta })={\left(1+\sum _{k=1}^{M-1}\exp ({\eta }_k)\right)}^{-1}\exp ({\mathbf{\eta }}^{\mathrm{T}}\mathbf{x})\text{\hspace{1em}(21)}$$

  这是带有参数向量 $\mathbf{\eta }=({\eta }_1,\cdots ,{\eta }_{M-1}{)}^{\mathrm{T}}$ 的指数族分布的标准形式。与式 (1) 比较可得

$$\mathbf{u}(\mathbf{x})=\mathbf{x}\text{\hspace{1em}(22)}$$

$$h(\mathbf{x})=1\text{\hspace{1em}(23)}$$

$$g(\mathbf{\eta })={\left(1+\sum _{k=1}^{M-1}\exp ({\eta }_k)\right)}^{-1}\text{\hspace{1em}(24)}$$

  最后考虑高斯分布。对于一元高斯分布，我们有

$$p(x\mid \mu ,{\sigma }^2)=\frac{1}{(2\pi {\sigma }^2{)}^{1/2}}\exp \left\{-\frac{1}{2{\sigma }^2}(x-\mu {)}^2\right\}\text{\hspace{1em}(25)}$$

$$=\frac{1}{(2\pi {\sigma }^2{)}^{1/2}}\exp \left\{-\frac{1}{2{\sigma }^2}{x}^2+\frac{\mu }{{\sigma }^2}x-\frac{1}{2{\sigma }^2}{\mu }^2\right\}\text{\hspace{1em}(26)}$$

  经过一些简单的重排，就可以得到标准指数族分布的形式：

$$\mathbf{\eta }=\left(\begin{array}{c}\mu /{\sigma }^2\\ -1/2{\sigma }^2\end{array}\right)\text{\hspace{1em}(27)}$$

$$\mathbf{u}(x)=\left(\begin{array}{c}x\\ x^2\end{array}\right)\text{\hspace{1em}(28)}$$

$$h(x)=(2\pi {)}^{-1/2}\text{\hspace{1em}(29)}$$

$$g(\mathbf{\eta })=(-2{\eta }_2{)}^{1/2}\exp \left(\frac{{\eta }_1^2}{4{\eta }_2}\right)\text{\hspace{1em}(30)}$$

  指数族分布在参数估计时具有很实际的性质，即它天然给出最小充分统计量（这来源于费曼分解定理），做参数估计时，只需要保存充分统计量，不需要存全部原始数据，大数据场景极致节省内存。

NO.8.tip： 核密度估计

背景

  参数话的密度估计方法的一个关键局限是，所选择的密度模型可能是不适合数据分布的，这可能导致预测性能不佳。而非参数化的密度估计方法对分布的形式几乎不作任何假设。

  直方图是非参数的一个良好灵感来源。标准直方图简单将 $x$ 划分为互不相交的宽度为 ${\Delta }_i$ 的分箱，然后计算落在不同分箱中的观测值 $x$ 的数量 $n_i$。为了将此计数值转换为归一化的概率密度，我们可以将其简单地除以观测值总数 $N$ 和分箱宽度 ${\Delta }_i$ 的乘积，以获得落入每个分箱的概率值：

$$p_i=\frac{n_i}{N{\Delta }_i}(1)$$

可以看出，$\int p(x)\mathrm{d}x=1$。这种方法就给出了密度 $p(x)$ 的一种模型，在每个分箱的宽度上，密度是恒定的。通常情况下，我们选择具有相同宽度的分箱，${\Delta }_i=\Delta$。

  直方图方法也给密度估计带来了两条重要的经验。首先，要在特定位置估计概率密度，就应该考虑落在该位置某个局部邻域内的数据点。其次，为了获得良好的结果，平滑参数的值既不应该太大也不应该太小。

  假设观测值是从某个 $D$ 维空间（这里将其视为欧氏空间）中的某个未知概率密度 $p(x)$ 中产生的，我们希望估计 $p(x)$ 的值。根据先前关于局部性的讨论，考虑包含 $x$ 的一个小区域 $\mathcal{R}$。与该区域相关联的概率质量为

$$P={\int }_{\mathcal{R}}p(x)\mathrm{d}x(2)$$

  假设我们已经收集了一个包含 $N$ 个从 $p(x)$ 中产生的观测值的数据集。因为每个数据点在 $\mathcal{R}$ 内的概率为 $P$，所以 $\mathcal{R}$ 内的点的总数 $K$ 将服从二项分布：

$$\mathrm{Bin}(K\mid N,P)=\frac{N!}{K!(N-K)!}{P}^K(1-P{)}^{N-K}(3)$$

  可以看出，落入该区域的点的平均占比为 $\mathbb{E}[K/N]=P$。类似地，可以看出，此均值对应的方差为 $\mathrm{var}[K/N]=P(1-P)/N$。对于较大的 $N$，该分布在均值处将变得十分尖锐，因此

$$K\approx NP(4)$$

  然而，如果我们同时假设区域 $\mathcal{R}$ 足够小，以至于概率密度 $p(x)$ 在该区域内大致保持恒定，则有

$$P\approx p(x)V(5)$$

  其中 $V$ 是区域 $\mathcal{R}$ 的体积。结合式 (4) 和式 (5)，可以得到密度估计的形式为

$$p(x)=\frac{K}{NV}(6)$$

  请注意，式 (6) 的有效性依靠两个相互矛盾的假设：区域 $\mathcal{R}$ 足够小，以至于密度在该区域内近似恒定；但区域 $\mathcal{R}$ 又足够大（与该区域内密度的值相关），使得落入该区域的点的数量 $K$ 足以使二项分布变得十分尖锐。

  我们可以通过两种不同的方式利用式 (6)。一种是固定 $K$ 并根据数据确定 $V$ 的值，这将引出后面讨论的 $K$ 近邻技术；另一种是固定 $V$ 并根据数据确定 $K$ 的值，从而引出核密度估计技术。可以证明，当 $N\to \infty$ 时，只要 $V$ 随着 $N$ 以合适的速率缩小，而 $K$ 随着 $N$ 以合适的速率增长，$K$ 近邻密度估计和核密度估计都会收敛到真实概率密度。

源与流

  下面详细讨论核密度估计技术。首先，对于希望确定概率密度的点 $x$，将区域 $\mathcal{R}$ 取为以点 $x$ 为中心的超小立方体。为了计算落入该区域的点的数量 $K$，定义如下函数：

$$k(u)=\{\begin{array}{ll}1,& |u_i|\le 1/2,i=1,\dots ,D\\ 0,& \text{其他}\end{array}(7)$$

  它表示以原点为中心的单位立方体。函数 $k(u)$ 是一个核函数，在这里，它又称为 Parzen 窗。根据式 (7)，如果数据点 $x_n$ 位于以点 $x$ 为中心、边长为 $h$ 的立方体内，那么量 $k\left(\frac{x-x_n}{h}\right)$ 将为 1，否则为 0。因此，位于该立方体内的数据点的总数为

$$K=\sum _{n=1}^{N}k\left(\frac{x-x_n}{h}\right)(8)$$

  将式 (8) 代入式 (6)，可得点 $x$ 处估计的密度为

$$p(x)=\frac{1}{N}\sum _{n=1}^N\frac{1}{h^D}k\left(\frac{x-x_n}{h}\right)(9)$$

  这里使用了 $D$ 维空间中超小立方体的体积公式 $V=h^D$。利用函数 $k(u)$ 的对称性，我们可以不再将其视为以点 $x$ 为中心的单个立方体，而是视为以 $N$ 个数据点 $x_n$ 为中心的 $N$ 个立方体的总和，进而重新解释这个方程。

  就目前情况来看，核密度估计 [式 (9)] 将遇到与直方图方法相同的一个问题，即存在人为的不连续性，此处体现在立方体的边界上。如果我们选择更平滑的核函数，则可以获得更平滑的密度模型。我们通常选择使用高斯函数，导出的核密度模型为

$$p(x)=\frac{1}{N}\sum _{n=1}^N\frac{1}{(2\pi h^2{)}^{D/2}}\exp \left\{-\frac{\Vert x-x_n{\Vert }^2}{2h^2}\right\}(10)$$

  其中 $h$ 代表所使用的高斯分量的标准差。因此，我们得到密度模型的方式是，首先在每个数据点上放置一个高斯分布，并将整个数据集的贡献相加，然后除以 $N$ 来正确地归一化密度。正如预期的那样，我们可以看到参数 $h$ 充当了平滑参数的角色，并且 $h$ 较小时对噪声敏感，与 $h$ 较大时过度平滑之间存在权衡。同样，类似于直方图密度估计中的分箱选择或曲线拟合中所使用多项式的阶数选择，$h$ 的优化是模型复杂性的问题。

  我们可以选择任何其他核函数 $k(u)$，只要满足以下两个条件即可。

$$k(u)\ge 0(11)$$

$$\int k(u)\mathrm{d}u=1(12)$$

  这两个条件确保了得到的概率分布在任何地方都是非负的，并且积分为 1。由式 (9) 给出的这类密度模型称为核密度估计器或 Parzen 估计器。它们有一个很大的优点，就是在"训练"阶段不涉及计算，因为只需要将训练集存储起来即可。然而，这同时也是这类密度模型的一个巨大缺点，因为评估密度的计算成本会随着数据集的增大而线性增长。

NO.9.tip： K近邻密度估计

背景

  使用核方法进行密度估计的一个难点在于，对所有核来说控制核宽度的参数 $h$ 是固定的。在数据密度较高的区域，较大的 $h$ 值可能导致过度平滑，并且可能淡化从数据中本应提取出的结构。然而，减小 $h$ 值可能会导致在其他数据密度较低的区域产生噪声估计。因此，$h$ 的最佳选择可能取决于数据空间中的位置。这个问题可以通过使用最近邻方法进行密度估计来解决。

源与流

  因此，下面我们重新回到局部密度估计的一般结果： $$p(x)=\frac{K}{NV}(P=p(x)V,K=NP)\text{\hspace{1em}(1)}$$

  我们不使用固定 $V$ 值并从数据中确定 $V$ 值的方式，而是考虑固定 $K$ 值并使用数据找到适当的 $V$ 值。为此，对于希望估计密度 $p(\mathbf{x})$ 的点 $\mathbf{x}$，考虑以点 $\mathbf{x}$ 为中心的小球，并允许小球的半径增长，直至恰好包含 $K$ 个数据点。然后，密度 $p(\mathbf{x})$ 的估计由式（1）给出，其中 $V$ 设置为小球的体积。这种技术称为 $K$ 近邻。

  相同数据集下，我们可以选择不同参数 $K$ 时的 $K$ 近邻密度估计。我们可以看到 $K$ 值决定了模型的平滑程度，而且同样存在一个既不太大也不太小的最佳 $K$ 值。注意 $K$ 近邻产生的模型并不是一个真正的密度模型，因为其在整个空间上的积分是发散的。

NO.10.tip： AIC和BIC准则

背景

  一般来说，训练误差

$$\overline{\mathrm{err}}=\frac{1}{N}\sum _{i=1}^{N}L(y_i,\hat{f}(x_i))\text{\hspace{1em}(1)}$$

比真实的预测误差 ${\mathrm{Err}}_{\mathcal{T}}$ 小，因为同样的数据既被用于拟合，又被用于评估其误差（见习题7.9）。拟合算法往往会对训练数据本身（而不是真实的分布）产生适应性，因此训练误差 $\overline{\mathrm{err}}$ 将会是对泛化误差 ${\mathrm{Err}}_{\mathcal{T}}$ 的过于乐观的估计。

  两者产生差异的部分原因是测试观测的位置或评估点发生位置的不一样。${\mathrm{Err}}_{\mathcal{T}}$ 可以认为是样本集外的样本（extra-sample）误差，因为测试样本不必和训练样本一样。$\overline{\mathrm{err}}$ 的乐观本质在我们仅关注样本集内的样本（in-sample）误差时最容易理解

$${\mathrm{Err}}_{\mathrm{in}}=\frac{1}{N}\sum _{i=1}^N{\mathbb{E}}_{Y^0}[L(Y_i^0,\hat{f}(x_i))|\mathcal{T}]\text{\hspace{1em}(2)}$$

这里 $Y^0$ 表示我们在每个训练样本 $x_i$ 获得的 $N$ 个新响应。我们将（误差估计的）乐观程度（optimism）定义为训练误差 ${\mathrm{Err}}_{\mathrm{in}}$ 和 $\overline{\mathrm{err}}$ 的差

$$\mathrm{op}\equiv {\mathrm{Err}}_{\mathrm{in}}-\overline{\mathrm{err}}\text{\hspace{1em}(3)}$$

这个值一般都是正的，因为 $\overline{\mathrm{err}}$ 通常是对预测误差偏小的估计。最后，平均乐观程度是乐观程度在所有训练集上的期望误差

$$\omega \equiv {\mathbb{E}}_y(\mathrm{op})\text{\hspace{1em}(4)}$$

这里，训练集的预测子是固定的，期望是对训练集上响应或结果值上（即 $Y_i^0$）取的，因此我们的记号是 ${\mathbb{E}}_y$ 而不是 ${\mathbb{E}}_{\mathcal{T}}$。通常，我们可以计算出期望乐观程度 $\omega$ 而不是 $\mathrm{op}$，正如我们能够估计期望预测误差 $\mathrm{Err}$，而不能估计出条件误差 ${\mathrm{Err}}_{\mathcal{T}}$。

  对于平方误差、0-1以及其他的损失函数，可以证明如下一般性的结论

$$\omega =\frac{2}{N}\sum _{i=1}^N\mathrm{Cov}({\hat{y}}_i,y_i)\text{\hspace{1em}(5)}$$

其中 $\mathrm{Cov}$ 是协方差。因此，$\overline{\mathrm{err}}$ 低估真实误差的程度依赖于 $y_i$ 能以多强的程度影响它自己的预测值。我们对数据的拟合程度越高，$\mathrm{Cov}({\hat{y}}_i,y_i)$ 会越大，也将因此而增大乐观程度。

  由此，我们有如下重要的关系

$${\mathbb{E}}_y({\mathrm{Err}}_{\mathrm{in}})={\mathbb{E}}_y(\overline{\mathrm{err}})+\frac{2}{N}\sum _{i=1}^N\mathrm{Cov}({\hat{y}}_i,y_i)\text{\hspace{1em}(6)}$$

  如果 ${\hat{y}}_i$ 是 $d$ 维输入或基函数的线性拟合结果，此表达式还可简化。例如，对于加性误差模型 $Y=f(X)+\epsilon$

$$\sum _{i=1}^N\mathrm{Cov}({\hat{y}}_i,y_i)=d{\sigma }_{\epsilon }^2\text{\hspace{1em}(7)}$$

因此

$${\mathbb{E}}_y({\mathrm{Err}}_{\mathrm{in}})={\mathbb{E}}_y(\overline{\mathrm{err}})+2\cdot \frac{d}{N}{\sigma }_{\epsilon }^2\text{\hspace{1em}(8)}$$

式(7)是有效参数个数的基础。乐观程度随着输入维数或者使用的基函数个数 $d$ 线性增长，但是随训练集样本数目的增加而减少。式(8)对于其他的误差模型，比如二值数据与熵损失函数，也能近似的成立。

  一种估计预测误差的明显方式是估计乐观程度，然后将它加到训练误差 $\overline{\mathrm{err}}$ 上。在下一节中介绍的方法——$C_p$、AIC、BIC以及其他的方法——都是这样对一类特殊的估计起作用的，这类估计是它们的参数的线性函数。

源与流

  样本内误差的一般形式为

$${\widehat{\mathrm{Err}}}_{\mathrm{in}}=\overline{\mathrm{err}}+\hat{\omega }\text{\hspace{1em}(9)}$$

其中 $\hat{\omega }$ 是对平均乐观程度的估计。

  利用式(8)，当 $d$ 个参数通过平方误差损失函数拟合时，可计算出一个所谓的 $C_p$ 统计量：

$$C_p=\overline{\mathrm{err}}+2\cdot \frac{d}{N}{\hat{\sigma }}_{\epsilon }^2\text{\hspace{1em}(10)}$$

这里 ${\hat{\sigma }}_{\epsilon }^2$ 是对噪声方差的估计，它从低偏差模型的均方误差获得。利用这个准则，我们通过一个正比于使用的基函数数量的因子，来调整训练误差。

  Akaike 信息准则是对 ${\mathrm{Err}}_{\mathrm{in}}$ 的类似但更一般的估计，如果用的是对数似然损失函数。它依赖于与式(8)类似的关系。该关系

$$-2\cdot \mathbb{E}\log {\mathrm{Pr}}_{\hat{\theta }}(Y)\approx -\frac{2}{N}\cdot \mathbb{E}[\mathrm{loglik}]+2\cdot \frac{d}{N}\text{\hspace{1em}(11)}$$

在 $N\to \infty$ 时，渐进的成立。这里 ${\mathrm{Pr}}_{\hat{\theta }}(Y)$ 表示一族 $Y$ 的密度（包含"真实"的密度函数），$\hat{\theta }$ 是对 $\theta$ 的最大似然估计，而 $\mathrm{loglik}$ 是最大化的似然函数

$$\mathrm{loglik}=\sum _{i=1}^N\log {\mathrm{Pr}}_{\hat{\theta }}(y_i)\text{\hspace{1em}(12)}$$

例如，对 logistic 回归模型，使用二项分布的对数似然函数，有

$$\mathrm{AIC}=-\frac{2}{N}\cdot \mathrm{loglik}+2\cdot \frac{d}{N}\text{\hspace{1em}(13)}$$

对于线性回归模型，其损失函数为平方误差（其方差 ${\sigma }_{\epsilon }^2={\hat{\sigma }}_{\epsilon }^2$ 假定已知），对应的 AIC 统计量等价于 $C_p$，因此我们将它们统称为 AIC。

  为了将 AIC 用于模型选择，我们就选择给定模型中 AIC 值最小的那个。对于非线性和其他复杂的模型，需要将 $d$ 替换为其他的能够衡量模型复杂性的数值。

  给定一族模型 $f_{\alpha }(x)$，其中 $\alpha$ 是可调整的参数，用 $\overline{\mathrm{err}}(\alpha )$ 和 $d(\alpha )$ 表示训练误差以及对应模型的参数个数，则对这族模型，我们定义

$$\mathrm{AIC}(\alpha )=\overline{\mathrm{err}}(\alpha )+2\cdot \frac{d(\alpha )}{N}{\hat{\sigma }}_{\epsilon }^2\text{\hspace{1em}(14)}$$

函数 $\mathrm{AIC}(\alpha )$ 提供对测试误差曲线的估计，我们去寻找能最小化该值的可调参数 $\hat{\alpha }$。

  尽管 AIC 和 BIC 是非常相似的，BIC 源自一个非常不同的起因，模型选择的贝叶斯方法。

  假定有一族候选模型 ${\mathcal{M}}_m,m=1,\dots ,M$ 且对应的模型参数为 ${\theta }_m$，我们希望从中选择一个最佳模型。假定每个模型 ${\mathcal{M}}_m$ 的参数的先验分布为 $\mathrm{Pr}({\theta }_m|{\mathcal{M}}_m)$，则给定模型的后验概率为

$$\mathrm{Pr}({\mathcal{M}}_m|Z)\propto \mathrm{Pr}({\mathcal{M}}_m)\cdot \mathrm{Pr}(Z|{\mathcal{M}}_m)\propto \mathrm{Pr}({\mathcal{M}}_m)\cdot \int \mathrm{Pr}(Z|{\theta }_m,{\mathcal{M}}_m)\mathrm{Pr}({\theta }_m|{\mathcal{M}}_m)d{\theta }_m\text{\hspace{1em}(21)}$$

其中 $Z$ 代表训练数据 $\{x_i,y_i{\}}_1^N$。为了比较两个模型 ${\mathcal{M}}_m$ 和 ${\mathcal{M}}_{\ell }$，我们计算后验概率比

$$\frac{\mathrm{Pr}({\mathcal{M}}_m|Z)}{\mathrm{Pr}({\mathcal{M}}_{\ell }|Z)}=\frac{\mathrm{Pr}({\mathcal{M}}_m)}{\mathrm{Pr}({\mathcal{M}}_{\ell })}\cdot \frac{\mathrm{Pr}(Z|{\mathcal{M}}_m)}{\mathrm{Pr}(Z|{\mathcal{M}}_{\ell })}\text{\hspace{1em}(22)}$$

如果比值大于 1，我们就选择模型 $m$，否则我们选择模型 $\ell$。最右边的一项

$$\mathrm{BF}(Z)=\frac{\mathrm{Pr}(Z|{\mathcal{M}}_m)}{\mathrm{Pr}(Z|{\mathcal{M}}_{\ell })}\text{\hspace{1em}(23)}$$

称为贝叶斯因子（Bayes factor），即数据对后验概率比的贡献。

  一般情况下，我们假设模型的先验是均匀分布，这样 $\mathrm{Pr}({\mathcal{M}}_m)$ 是常数。我们需要一些近似 $\mathrm{Pr}(Z|{\mathcal{M}}_m)$ 的方法。通过对积分的所谓拉普拉斯近似以及另外一些简化，式(21)化简为

$$\log \mathrm{Pr}(Z|{\mathcal{M}}_m)=\log \mathrm{Pr}(Z|{\hat{\theta }}_m,{\mathcal{M}}_m)-\frac{d_m}{2}\cdot \log N+O(1)\text{\hspace{1em}(24)}$$

这里 ${\hat{\theta }}_m$ 是最大似然估计，$d_m$ 是模型 ${\mathcal{M}}_m$ 的自由参数个数。如果损失函数定义为

$$-2\log \mathrm{Pr}(Z|{\hat{\theta }}_m,{\mathcal{M}}_m)$$

  因此，选择最小 BIC 的模型等价于选择能够（近似）最大后验概率的模型。但这个框架给了我们更多的东西。如果我们对 $M$ 个模型计算 BIC，有相应的 ${\mathrm{BIC}}_m,m=1,2,\dots ,M$，那么我们可以用它们来估计每个模型 ${\mathcal{M}}_m$ 的后验概率

$$\frac{e^{-\frac{1}{2}\cdot {\mathrm{BIC}}_m}}{{\sum }_{\ell =1}^M{e}^{-\frac{1}{2}\cdot {\mathrm{BIC}}_{\ell }}}\text{\hspace{1em}(25)}$$

这样我们不仅估计了最好的模型，同时还评估了各个模型的相对好处。

  就模型选择的目的而言，在 AIC 和 BIC 之间并没有明确的倾向。BIC 作为选择的准则是渐进一致的。这也就是说，给定一族模型，包含真实的模型，BIC 选择正确模型的概率会随着样本容量 $N\to \infty$ 接近 1。但是对 AIC 却并不是这样，它在 $N\to \infty$ 的时候倾向选择过于复杂的模型。另外一方面，对有限样本，BIC 经常由于对复杂模型的严重惩罚而选择过于简单的模型。

NO.11.tip： 假设检验

背景

  假设 $X$ 是随机变量，在 $\{1,\cdots ,n\}$ 中取值，其概率密度分布和参数 $\theta \in \{1,\cdots ,m\}$ 的取值有关。对 $\theta$ 的 $m$ 个可能值，$X$ 的概率分布可以由矩阵 $P\in {\mathbf{R}}^{n\times m}$ 表征，其元素为

$$p_{kj}=\mathrm{prob}(X=k\mid \theta =j).\text{\hspace{1em}(1)}$$

矩阵 $P$ 的第 $j$ 列为对应参数值 $\theta =j$ 的概率分布。

  考虑基于 $X$ 的观测样本估计 $\theta$ 的问题。换言之，样本 $X$ 来自于 $m$ 种可能的概率分布，我们需要确定是哪个。$\theta$ 的 $m$ 个取值称为假设，我们从这些假设中猜想哪个是正确的（即产生观测样本 $X$ 的概率分布），这个问题称为假设检验。在很多情况下，某种假设对应一种常规情形，而其他假设均对应反常的事件。此时，假设检验可以看成观测到 $X$ 的某个取值，然后判断不寻常事件是否发生，如果发生了，是哪一个不寻常事件。因此，假设检验又称为检测。

  在很多情况下，假设的顺序对结果没有什么影响；对 $m$ 个不同的假设，将其任意标识为 $\theta =1,\cdots ,m$。如果 $\hat{\theta }=\theta$，其中 $\hat{\theta }$ 表示 $\theta$ 的估计值，那么我们成功地猜想出了参数值 $\theta$。如果 $\hat{\theta }\ne \theta$，我们关于参数值 $\theta$ 的猜想不正确；我们误认为 $\theta$ 是 $\hat{\theta }$。在其他情况下，假设的顺序非常重要。在这种情况下，事件 $\hat{\theta }>\theta$，即我们关于 $\theta$ 的估计值偏大，对结果是有意义的。

  给定 $X$ 的观测值，$\theta$ 的估计值是随机的。$\theta$ 的一个随机检测器是随机变量 $\hat{\theta }\in \{1,\cdots ,m\}$，其分布取决于 $X$ 的观测值。随机检测器可以定义为一个矩阵 $T\in {\mathbf{R}}^{m\times n}$，其元素为

$$t_{ik}=\mathrm{prob}(\hat{\theta }=i\mid X=k).\text{\hspace{1em}(3)}$$

可以这么理解上述定义：如果我们得到观测值 $X=k$，那么检测器以概率 $t_{ik}$ 给出估计值 $\hat{\theta }=i$。（给定观测值 $X=k$，$T$ 的第 $k$ 列，我们用 $t_k$ 表示，给出了 $\hat{\theta }$ 的概率分布。如果 $T$ 的每一列都是单位向量，那么随机检测器是确定性检测器，即 $\hat{\theta }$ 是 $X$ 的观测值的一个（确定性）函数。）

  初步看来，在估计或检测过程中引入随机因素似乎只会让估计效果更差。但是在后文我们给出例子，在例子中随机检测器比所有的确定性估计器效果都好。

  我们致力于寻找定义随机检测器的矩阵 $T$。显然，$T$ 的列 $t_k$ 必须满足（线性等式和不等式）约束

$$t_k⪰0,1^T{t}_k=1.\text{\hspace{1em}(4)}$$

  对于由矩阵 $T$ 决定的随机检测器，定义检测概率矩阵为 $D=TP$。我们有

$$D_{ij}=(TP{)}_{ij}=\mathrm{prob}(\hat{\theta }=i\mid \theta =j),\text{\hspace{1em}(5)}$$

因此，$D_{ij}$ 是当 $\theta =j$ 时，猜想为 $\hat{\theta }=i$ 的概率。$m\times m$ 检测概率矩阵 $D$ 衡量了由矩阵 $T$ 定义的随机检测器的性能。（对角元素 $D_{ii}$ 是当 $\theta =i$ 时，猜想为 $\hat{\theta }=i$ 的概率，即正确地检测出 $\theta =i$ 的概率。非对角线元素 $D_{ij}(i\ne j)$ 是 $\theta =j$ 时误判为 $\theta =i$ 的概率，即事实上 $\theta =j$，而我们的猜想为 $\hat{\theta }=i$ 的概率。如果 $D=I$，检测器是理想的：无论参数 $\theta$ 的值如何，我们的猜想都是正确的，$\hat{\theta }=\theta$。）

  $D$ 的对角线元素，排列成一个向量，称为检测概率，用 $P^{\mathrm{d}}$ 表示，即

$$P_i^{\mathrm{d}}=D_{ii}=\mathrm{prob}(\hat{\theta }=i\mid \theta =i).\text{\hspace{1em}(6)}$$

错误概率是上面概率的补，用 $P^{\mathrm{e}}$ 表示，即

$$P_i^{\mathrm{e}}=1-D_{ii}=\mathrm{prob}(\hat{\theta }\ne i\mid \theta =i).\text{\hspace{1em}(7)}$$

因为检测概率矩阵 $D$ 的每列和为 $1$，错误概率可以表示为

$$P_i^{\mathrm{e}}=\sum _{j\ne i}{D}_{ji}.\text{\hspace{1em}(8)}$$

源与流

  检测器设计问题中有很大一部分目标函数是 $D$，也是 $T$（即优化变量）的线性、仿射或者凸分片线性函数。类似地，检测器设计问题中很多约束都可以表示成 $D$ 的线性不等式。因此，很多最优检测器设计问题都可以表述成线性规划。

  最优检测器设计问题可以看成一个多准则优化问题，约束为式 (4)，$m(m-1)$ 个目标由矩阵 $D$ 的非对角线元素给出，这些元素表征了不同类型的检测错误的概率：

$$\begin{array}{ll}\text{minimize}(\text{关于}{\mathbf{R}}_{+}^{m(m-1)})& D_{ij},i,j=1,\cdots ,m,i\ne j\\ \text{subject to}& t_k⪰0,1^T{t}_k=1,k=1,\cdots ,n,\end{array}\text{\hspace{1em}(14)}$$

其中优化变量 $t_1,\cdots ,t_n\in {\mathbf{R}}^m$。因为每个目标 $D_{ij}$ 都是优化变量的线性函数，这是一个多准则线性规划问题。

  为了将这个多准则优化问题标量化，引入权系数将目标函数加权求和，

$$\sum _{i,j=1}^m{W}_{ij}{D}_{ij}=\mathrm{tr}(W^{T}D),\text{\hspace{1em}(15)}$$

其中权矩阵 $W\in {\mathbf{R}}^{m\times m}$ 满足

$$W_{ii}=0,i=1,\cdots ,m,W_{ij}>0,i,j=1,\cdots ,m,i\ne j.\text{\hspace{1em}(16)}$$

此时目标函数是 $m(m-1)$ 个误差概率的加权和。权系数 $W_{ij}$ 对应着下述事件，事实上 $\theta =j$ 但猜想为 $\hat{\theta }=i$。加权矩阵有时亦称为损失矩阵。

为了找到多准则优化问题 (14) 的一个 Pareto 最优点，我们求解标量优化问题

$$\begin{array}{ll}\text{minimize}& \mathrm{tr}(W^{T}D)\\ \text{subject to}& t_k⪰0,1^T{t}_k=1,k=1,\cdots ,n,\end{array}\text{\hspace{1em}(17)}$$

它是一个线性规划。这个线性规划对变量 $t_1,\cdots ,t_n$ 是可分的。目标函数可以描述为一系列 $t_k$ 的（线性）函数的和

$$\mathrm{tr}(W^{T}D)=\mathrm{tr}(W^{T}TP)=\mathrm{tr}(P{W}^{T}T)=\sum _{k=1}^n{c}_k^T{t}_k,\text{\hspace{1em}(18)}$$

其中 $c_k$ 是 $W{P}^T$ 的第 $k$ 列。约束是可分的（即我们可以得到关于每个 $t_i$ 的可分的约束）。因此我们可以通过分别求解下面的线性规划来求解线性规划 (17)

$$\begin{array}{ll}\text{minimize}& c_k^T{t}_k\\ \text{subject to}& t_k⪰0,1^T{t}_k=1,\end{array}\text{\hspace{1em}(19)}$$

每个线性规划都有一个简单的解析解。我们可以找到一个指标 $q$ 使得 $c_{kq}={\min }_j{c}_{kj}$。令 $t_k^{\star }=e_q$。这个最优点对应一个确定性检测器：即当已知观测值 $X=k$ 时，我们的估计值为

$$\hat{\theta }={\mathrm{argmin}}_j(W{P}^T{)}_{jk}.\text{\hspace{1em}(20)}$$

因此，对每一个非对角线元素大于零的权矩阵 $W$ 我们可以找到一个确定性检测器，使得加权求和之后的目标函数最小。这似乎意味着随机检测器并没有意义，在后文中我们将看到事实并不是这样的。多目标优化问题 (14) 的 Pareto 最优权衡表面是分片线性的；式 (20) 描述的确定性检测器对应了 Pareto 最优权衡表面的顶点。

  作为一个例子，我们考虑 $m=2$ 的特殊情形，也称为二值假设检验。随机变量 $X$ 可能服从两种随机分布中的一种，为了简单起见，这两种随机分布用 $p\in {\mathbf{R}}^n$ 和 $q\in {\mathbf{R}}^n$ 表示。大部分情况下，假设 $\theta =1$ 对应一些常规的情形，$\theta =2$ 对应一些异常事件，而我们要检测出来这些异常事件。如果 $\hat{\theta }=1$，我们称检验是阴性的（即我们认为异常事件并没有发生）；如果 $\hat{\theta }=2$，我们称检验是阳性的（即我们认为异常事件发生了）。

  通常可以将检测概率矩阵 $D\in {\mathbf{R}}^{2\times 2}$ 表示为

$$D=\left[\begin{array}{cc}1-P_{\mathrm{fp}}& P_{\mathrm{fn}}\\ P_{\mathrm{fp}}& 1-P_{\mathrm{fn}}\end{array}\right].\text{\hspace{1em}(26)}$$

这里 $P_{\mathrm{fn}}$ 是假阴性的概率（即检验是阴性的但事实上异常事件发生了），$P_{\mathrm{fp}}$ 是假阳性的概率（即检验是阳性的但事实上异常事件并没有发生），亦被称为误报概率。最优检测器设计问题是一个双准则优化问题，目标函数为 $P_{\mathrm{fn}}$ 和 $P_{\mathrm{fp}}$。

  $P_{\mathrm{fn}}$ 和 $P_{\mathrm{fp}}$ 的最优权衡曲线称为接收器工作特性（ROC），由 $p$ 和 $q$ 服从的概率分布决定。可以用 §7.3.4 中的方法，标量化双准则问题并进行求解得到 ROC。给定权矩阵 $W$，最优检测器 (20) 为

$$\hat{\theta }=\{\begin{array}{ll}1& W_{21}{p}_k>W_{12}{q}_k\\ 2& W_{21}{p}_k⩽W_{12}{q}_k\end{array}\text{\hspace{1em}(27)}$$

其中观测值为 $X=k$。称这种检测器为似然比阈值检验：如果比值 $p_k/q_k$ 大于阈值 $W_{12}/W_{21}$，检验是阴性的（即 $\hat{\theta }=1$）；否则，检验是阳性的。通过选择不同的阈值，我们可以得到不同的（确定性）Pareto 最优检测器，这些检测器分别对应了不同水平的假阳性和假阴性错误概率。上述结论称为 Neyman-Pearson 引理。