8.2 自助法和最大似然方法

8.2.1 示例:函数平滑

自助法通过在训练数据上的抽样为评估不确定性提供了一个直接的计算方式。本节演示自助法在一个简单的一维平滑问题中的应用,并说明其与最大似然的联系。

训练数据记为 $\mathbf{Z}=\{z_1,z_2,\dots,z_N\}$,其中 $z_i=(x_i,y_i)$,$i=1,2,\dots,N$。这里的 $x_i$ 是一维的输入变量,$y_i$ 是连续或分类的输出变量。作为一个例子,考虑图 8.1 的左图中的 $N=50$ 个数据点。

**图 8.1**:左图为函数平滑例子中的数据。右图为七个 B-样条基函数。垂直的虚线标记了三个结点的位置。
图 8.1:左图为函数平滑例子中的数据。右图为七个 B-样条基函数。垂直的虚线标记了三个结点的位置。

假设我们决定用三次样条来拟合数据,以 $X$ 取值的四分位数作为结点位置。这是一个七维的函数线性空间,可以表述为诸如 B-样条基函数的线性展开(见第 5.9.2 节):

$$\mu(x) = \sum_{j=1}^7 \beta_j h_j(x) \tag{8.1}$$

式中的 $h_j(x)$,$j=1,2,\dots,7$ 为图 8.1 的右图中展示的七个函数。可以将 $\mu(x)$ 理解为代表了条件期望 $\operatorname{E}(Y|X=x)$。

令 $\mathbf{H}$ 为 $ N\times7$ 的矩阵,它的第 ij 个元素为 $h_j(x_i)$。通常会通过在训练集上最小化误差平方和来获得 $\beta$ 的估计,其解为:

$$\hat{\beta} = (\mathbf{H}^T\mathbf{H})^{-1} \mathbf{H}^T\mathbf{y} \tag{8.2}$$
**图 8.2**:左上图为数据的 B-样条平滑曲线。右上图为 B-样条平滑加减 1.96 倍标准误差的区间带。左下图为十次自助法重复计算的 B-样条平滑曲线。右下图为基于自助法分布产生的 B-样条平滑曲线和 95% 的标准误差区间带。
图 8.2:左上图为数据的 B-样条平滑曲线。右上图为 B-样条平滑加减 1.96 倍标准误差的区间带。左下图为十次自助法重复计算的 B-样条平滑曲线。右下图为基于自助法分布产生的 B-样条平滑曲线和 95% 的标准误差区间带。

图 8.2 的左上图给出了对应的拟合结果 $\hat{\mu}(x)=\sum_{j=1}^7\hat{\beta}h_j(x)$。

$\hat{\beta}$ 的协方差矩阵估计为:

$$\widehat{\operatorname{Var}}(\hat{\beta}) = (\mathbf{H}^T\mathbf{H})^{-1} \hat{\sigma}^2 \tag{8.3}$$

其中噪声方差的估计为 $\hat{\sigma}^2=\sum_{i=1}^N(y_i-\hat{\mu}(x_i))^2/N$。令 $h(x)^T=(h_1(x),h_2(x),\dots,h_7(x))$,则预测值 $\hat{\mu}(x)=h(x)^T\hat{\beta}$ 的标准误差为:

$$\widehat{\operatorname{se}}[\hat{\mu}(x)] = [h(x)^T (\mathbf{H}^T\mathbf{H})^{-1} h(x)]^{\frac{1}{2}} \hat{\sigma} \tag{8.4}$$

图 8.2 的右上图中画出了 $\hat{\mu}(x)\pm1.96\cdot\widehat{\operatorname{se}}[\hat{\mu}(x)]$。由于标准正态分布函数在 97.5% 对应的值是 1.96,所以这对应了大约 $100-2\times2.5%=95%$ 的 $\mu(x)$ 逐点置信区间带。

以下将在这个例子中应用自助法。先对训练数据进行有放回的抽样,得到 $B$ 组大小都为 $N=50$ 的训练集,其中的样本元素为 $z_i=(x_i,y_i)$。对每个自助集 $\mathbf{Z}^*$ 通过三次样条拟合得到 $\hat{\mu}^*(x)$;图 8.2 的左下图展示了十个这样生成的拟合曲线。在每个 $x$ 点处,用 $B=200$ 个自助样本可从其百分位数生成 95% 的逐点置信区间带,即在每个 $x$ 点处找出第 $2.5%\times 200=5$ 大和小的拟合值。图 8.2 的右下图画出了这个置信区间带。其与右上图中看起来相似,但在两端点处稍微宽一些。

实际上最小二乘估计的表达式 8.2 和 8.3、自助法和最大似然法之间存在紧密的联系。进一步假设模型的误差服从高斯分布:

$$\begin{align} Y &= \mu(X) + \varepsilon; \varepsilon \sim \mathcal{N}(0, \sigma^2) \\ \mu(x) &= \sum_{j=1}^7 \beta_j h_j(x) \tag{8.5}\end{align}$$

上述的自助法从训练数据中有放回的抽样,这也称为 非参数(nonparametric) 自助法。由于它直接使用原始数据而不是特定的参数模型来生成新数据集,也可理解为这个方法与模型无关。考虑另一种方式,称为 参数(parametric) 自助法,将预测值与高斯噪声相加来模拟新的输出变量:

$$y_i^* = \hat{\mu}(x_i) + \varepsilon_i^*; \varepsilon_i^* \sim \mathcal{N}(0, \hat{\sigma}^2); i = 1,2,\dots,N \tag{8.6}$$

重复 $B$ 次这个过程,例如 $B=200$。结果产生的自助集为 $(x_1,y_1^*)$,……,$(x_N,y_N^*)$,然后在每个自助集上重新计算 B-样条平滑。随着自助样本数量趋于无穷 $B\rightarrow\infty$,这种方法产生的置信区间带会与右上图中的最小二乘完全相等。通过自助样本 $\mathbf{y}^*$ 估计出的函数为 $\hat{\mu}^*(x)=h(x)^T(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{y}^*$,其分布为:

$$\hat{\mu}^*(x) \sim \mathcal{N}(\hat{\mu}(x), h(x)^T (\mathbf{H}^T\mathbf{H})^{-1} h(x) \hat{\sigma}^2) \tag{8.7}$$

注意这个分布的均值就等于最小二乘估计,而标准差也与近似表达式 8.4 一样。

8.2.2 最大似然推断

在上面的例子中,参数自助法之所以与最小二乘的结果一致,是因为模型(式 8.5)中的高斯误差是加性的。一般来说,参数自助法不一定与最小二乘的结果一致,但会与最大似然的结果一致,论述如下。

先对观测样本指定一个概率密度或质量函数:

$$z_i \sim g_\theta(z) \tag{8.8}$$

这个表达式中的 $\theta$ 代表了控制着 $Z$ 分布的一个或多个未知参数,这种形式被称为 $Z$ 的 参数模型(parametric model)。例如,若 $Z$ 服从均值为 $\mu$ 方差为 $\sigma^2$ 的正态分布,则:

$$\theta = (\mu, \sigma^2) \tag{8.9}$$

并且

$$g_\theta(z) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{1}{2} (z-\mu)^2 / \sigma^2} \tag{8.10}$$

最大似然方法基于 似然函数(likelihood function)

$$L(\theta; \mathbf{Z}) = \prod_{i=1}^N g_\theta(z_i) \tag{8.11}$$

即在模型 $g_\theta$ 下得到观测数据样本的概率。严格的似然函数定义中需要有一个正的乘子,在这里令其等于一,不会影响结果。将数据 $\mathbf{Z}$ 视为常数,则可将 $L(\theta;\mathbf{Z})$ 视为 $\theta$ 的函数。

$L(\theta;\mathbf{Z})$ 的对数可写为:

$$\begin{align} l(\theta; \mathbf{Z}) &= \sum_{i=1}^N l(\theta; z_i) \\ &= \sum_{i=1}^N \log g_\theta(z_i) \tag{8.12}\end{align}$$

有时会将其简写为 $\ell(\theta)$。这个表达式被称为对数似然度,每个值 $\ell(\theta;z_i)=\log g_\theta(z_i)$ 被称为对数似然成分(component)。最大似然方法寻找可以最大化 $\ell(\theta;\mathbf{Z})$ 的参数值 $\theta=\hat{\theta}$。

似然函数也可用来评估 $\hat{\theta}$ 的准确性。下面介绍几个定义。定义 得分函数(score function) 为:

$$\dot{\ell}(\theta; \mathbf{Z}) = \sum_{i=1}^N \dot{\ell}(\theta; z_i) \tag{8.13}$$

其中 $\dot{\ell}(\theta;z_i)=\partial\ell(\theta;z_i)/\partial\theta$。假设似然函数在其参数空间的内部取到最大值,则 $\dot{\ell}(\hat{\theta};\mathbf{Z})=0$。信息矩阵(information matrix) 为:

$$\mathbf{I}(\theta) = -\sum_{i=1}^N \frac{\partial^2 \ell(\theta; z_i)}{\partial\theta\partial\theta^T} \tag{8.14}$$

当取值在 $\theta=\hat{\theta}$ 点时,这也被称为 观测信息量(observed information)费雪信息量(Fisher information),或期望信息量,定义为:

$$\mathbf{i}(\theta) = \operatorname{E}_\theta[\mathbf{I}(\theta)] \tag{8.15}$$

最后,令 $\theta_0$ 代表参数 $\theta$ 的真实值。

一个常规的结论是,随着 $N\rightarrow\infty$,最大似然估计的样本分布会趋向于正态分布:

$$\hat{\theta} \rightarrow \mathcal{N} (\theta_0, \mathbf{i}(\hat{\theta}_0)^{-1}) \tag{8.16}$$

这里的样本是从 $g_{\theta_0}(z)$ 的独立抽样。这意味着 $\hat{\theta}$ 的样本分布可近似为:

$$\mathcal{N}(\hat{\theta}, \mathbf{i}(\hat{\theta})^{-1}) \text{ 或 } \mathcal{N}(\hat{\theta}, \mathbf{I}(\hat{\theta})^{-1}) \tag{8.17}$$

其中 $\hat{\theta}$ 代表从观测数据中得出的最大似然估计。

对应的 $\hat{\theta}_j$ 标准误差的估计为:

$$\sqrt{\mathbf{i}(\hat{\theta})_{jj}^{-1}} \text{ 和 } \sqrt{\mathbf{I}(\hat{\theta})_{jj}^{-1}} \tag{8.18}$$

$\theta_j$ 的置信区间点可用表达式 8.17 中的任一个近似生成,分别可写为:

$$\theta_j - z^{(1-\alpha)} \cdot \sqrt{\mathbf{i}(\hat{\theta})_{jj}^{-1}} \text{ 或 } \theta_j - z^{(1-\alpha)} \cdot \sqrt{\mathbf{I}(\hat{\theta})_{jj}^{-1}}$$

其中 $z^{(1-\alpha)}$ 是标准正态分布的 $1-\alpha$ 百分位数。使用卡方分布近似,可从似然函数推导出更准确的置信区间:

$$2[\ell(\hat{\theta}) − \ell(\theta_0)] \sim \chi_p^2 \tag{8.19}$$

其中 $p$ 为参数 $\theta$ 中元素的个数。得到的 $1-2\alpha$ 置信区间,是所有满足 $2[\ell(\hat{\theta}) − \ell(\theta_0)] \sim \chi_p^2$ 的 $\theta_0$ 的集合,其中 $\chi_p^2$ 是 $p$ 自由度卡方分布的 $1-2\alpha$ 百分位数。

回到函数平滑的例子中,考察应用最大似然的结果。参数为 $\theta=(\beta,\sigma^2)$。对数似然度为:

$$\ell(\theta) = -\frac{N}{2} \log \sigma^2 2\pi -\frac{1}{2\sigma^2} \sum_{i=1}^N (y_i - h(x_i)^T\beta)^2 \tag{8.20}$$

通过条件 $\partial\ell/\partial\beta=0$ 和 $\partial\ell/\partial\sigma^2=0$ 解出最大似然估计:

$$\begin{align} \hat{\beta} &= (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{y} \\ \hat{\sigma}^2 &= \frac{1}{N} \sum (y_i - \hat{\mu}(x_i))^2 \end{align}\tag{8.21}$$

这与式 8.2 和 8.3 中通常的最小二乘估计是一样的。

参数 $\theta=(\beta,\sigma^2)$ 的信息矩阵是分块对角矩阵,对应着 $\beta$ 的子矩阵是:

$$\mathbf{I}(\beta) = (\mathbf{H}^T\mathbf{H}) / \sigma^2 \tag{8.22}$$

因此方差的估计 $(\mathbf{H}^T\mathbf{H})/\hat{\sigma}^2$ 与最小二乘的估计(式 8.3)是一样的。

8.2.3 自助法与最大似然

自助法本质上是非参数或参数最大似然的计算机实现。自助法对最大似然表达式的优势在于其允许在无法得到表达式的场景中计算标准误差和其他量值的最大似然估计。

在之前的例子中,假设不预先确定 B-样条中的结点位置,而是通过交叉验证自适应地选择结点的数量和位置。将结点的集合和他们的位置记为 $\lambda$。那么标准误差和置信区间带应体现对 $\lambda$ 的自适应选择,然而从解析上无法得出表达式。通过自助法,可对每个自助样本计算含有自适应结点的 B-样条平滑。得到函数曲线的百分位数中同时包含了来自目标变量中的噪声以及结点选择 $\hat{\lambda}$ 的变动。在这个特定的例子中,是否预先确定结点的置信区间带差异不大(未在图中展示)。但在自适应选择更多的其他问题中,对置信区间的影响会比较大。

下一页