第十章 bootstrap 方法
第十章 bootstrap 方法
§1 非参数 bootstrap 方法
设总体的分布 F 未知,但已经有一个容量为 n 的来自分布 F 的数据样本,自这一样本按放回抽样的方法抽取一个容量为 n 的样本,这种样本称为 bootstrap 样本或称为自助样本。相继地、独立地自原始样本中取很多个 bootstrap 样本,利用这些样本对总体 F 进行统计推断。这种方法称为非参数 bootstrap 方法,又称自助法。这一方法可以用于当人们对总体知之甚少的情况,它是近代统计中的一种用于数据处理的重要实用方法。这种方法的实现需要在计算机上作大量的计算,随着计算机威力的增长,它已成为一种流行的方法。
bootstrap 方法是 Efron 在 20 世纪 70 年代后期建立的.
(一) 估计量的标准误差的 bootstrap 估计
在估计总体未知参数 $ \theta $ 时,人们不但要给出 $ \theta $ 的估计 $ \hat{\theta} $,还需指出这一估计 $ \hat{\theta} $ 的精度。通常我们用估计量 $ \hat{\theta} $ 的标准差 $ \sqrt{D(\hat{\theta})} $ 来度量估计的精度。估计量 $ \hat{\theta} $ 的标准差 $ \sigma_{\hat{\theta}} = \sqrt{D(\hat{\theta})} $ 也称为估计量 $ \hat{\theta} $ 的标准误差。
设 $ X_{1}, X_{2}, \cdots, X_{n} $ 是来自以 $ F(x) $ 为分布函数的总体的样本, $ \theta $ 是我们感兴趣的未知参数,用 $ \hat{\theta} = \hat{\theta}(X_{1}, X_{2}, \cdots, X_{n}) $ 作为 $ \theta $ 的估计量,在应用中 $ \hat{\theta} $ 的抽样分布常是很难处理的,这样, $ \sqrt{D(\hat{\theta})} $ 常没有一个简单的表达式,不过我们可以用计算机模拟的方法来求得 $ \sqrt{D(\hat{\theta})} $ 的估计. 为此,自 F 产生很多容量为 n 的样本(例如 B 个),对于每一个样本计算 $ \hat{\theta} $ 的值,得 $ \hat{\theta}_{1}, \hat{\theta}_{2}, \cdots, \hat{\theta}_{B} $,则 $ \sqrt{D(\hat{\theta})} $ 可以用
$$ \hat{\sigma}_{\theta}=\sqrt{\frac{1}{B-1}\sum_{i=1}^{B}(\hat{\theta}_{i}-\overline{\theta})^{2}} $$
来估计,其中 $ \overline{\theta} = \frac{1}{B} \sum_{i=1}^{B} \hat{\theta}_{i} $ 。然而 F 常常是未知的,这样就无法产生模拟样本,不能得到(1.1)式的结果,需要另外的方法。
现在设分布 F 未知, $ x_{1}, x_{2}, \cdots, x_{n} $ 是来自 F 的样本值, $ F_{n} $ 是相应的经验分
布函数. 当 $n$ 很大时, $F_n$ 接近 $F$ (见第六章 § 3 格里汶科定理). 我们用 $F_n$ 代替上一段中的 $F$, 在 $F_n$ 中抽样. 在 $F_n$ 中抽样, 就是在原始样本 $x_1, x_2, \cdots, x_n$ 中每次随机地取一个个体作放回抽样. 如此得到一个容量为 $n$ 的样本 $x_1^*, x_2^*, \cdots, x_n^*$. 这就是第一段中所说的 bootstrap 样本. 用 bootstrap 样本按上一段中计算估计 $\hat{\theta}(x_1, x_2, \cdots, x_n)$ 那样求出 $\theta$ 的估计 $\hat{\theta}^* = \hat{\theta}(x_1^*, x_2^*, \cdots, x_n^*)$, 估计 $\hat{\theta}^*$ 称为 $\theta$ 的 bootstrap 估计. 相继地、独立地抽得 $B$ 个 bootstrap 样本, 以这些样本分别求出 $\theta$ 的相应的 bootstrap 估计如下:
boostrap 样本 1 $ x_{1}^{*1}, x_{2}^{*1}, \cdots, x_{n}^{*1} $,bootstrap 估计 $ \hat{\theta}_{1}^{*} $
boostrap 样本 2 $ x_{1}^{*2}, x_{2}^{*2}, \cdots, x_{n}^{*2} $,bootstrap 估计 $ \hat{\theta}_{2}^{*} $
$ \vdots \quad \vdots \quad \vdots $
boostrap 样本 B $ x_{1}^{*B}, x_{2}^{*B}, \cdots, x_{n}^{*B} $,bootstrap 估计 $ \hat{\theta}_{B}^{*} $ 则 $ \hat{\theta} $ 的标准误差 $ \sqrt{D(\hat{\theta})} $,就以
$$ \hat{\sigma}_{\theta}=\sqrt{\frac{1}{B-1}\sum_{i=1}^{B}(\hat{\theta}_{i}^{*}-\overline{\theta}^{*})^{2}} $$
来估计,其中 $ \overline{\theta}^{*}=\frac{1}{B}\sum_{i=1}^{B}\hat{\theta}_{i}^{*} $.(1.2)式就是 $ \sqrt{D(\hat{\theta})} $的bootstrap估计.
综上所述得到求 $ \sqrt{D(\hat{\theta})} $的bootstrap估计的步骤是
$ 1^{\circ} $ 自原始数据样本 $ \boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n}) $ 按放回抽样的方法,抽得容量为 n 的样本 $ \boldsymbol{x}^{*}=(x_{1}^{*},x_{2}^{*},\cdots,x_{n}^{*}) $ (称为 bootstrap 样本);
$ 2^{\circ} $ 相继地、独立地求出 B 个 $ (B \geqslant 1000) $ 容量为 n 的 bootstrap 样本, $ x^{*i} = (x_{1}^{*i}, x_{2}^{*i}, \cdots, x_{n}^{*i}), i = 1, 2, \cdots, B $. 对于第 i 个 bootstrap 样本,计算 $ \hat{\theta}_{i}^{*} = \hat{\theta}(x_{1}^{*i}, x_{2}^{*i}, \cdots, x_{n}^{*i}), i = 1, 2, \cdots, B $ ( $ \hat{\theta}_{i}^{*} $ 称为 $ \theta $ 的第 i 个 bootstrap 估计.)
$ 3^{\circ} $ 计算
$$ \hat{\sigma}_{\theta}=\sqrt{\frac{1}{B-1}\sum_{i=1}^{B}(\hat{\theta}_{i}^{*}-\overline{\theta}^{*})^{2}}\ , 其中 \overline{\theta}^{*}=\frac{1}{B}\sum_{i=1}^{B}\hat{\theta}_{i}^{*} $$
例 1 某种基金的年回报率是具有分布函数 F 的连续型随机变量,F 未知,F 的中位数 $ \theta $ 是未知参数。现有以下的数据(%率):
18.2 9.5 12.0 21.1 10.2
以样本中位数作为总体中位数 $ \theta $ 的估计. 试求中位数估计的标准误差的 bootstrap 估计.
解 将原始样本自小到大排序,中间一个数为12.0,得样本中位数为12.0。相继地、独立地在上述5个数据中,按放回抽样的方法取样,取B=10得到
下述10个bootstrap样本:①
| 样本1 | 9.5 | 18.2 | 12.0 | 10.2 | 18.2 |
| 样本2 | 21.1 | 18.2 | 12.0 | 9.5 | 10.2 |
| 样本3 | 21.1 | 10.2 | 10.2 | 12.0 | 10.2 |
| 样本4 | 18.2 | 12.0 | 9.5 | 18.2 | 10.2 |
| 样本5 | 21.1 | 12.0 | 18.2 | 12.0 | 18.2 |
| 样本6 | 10.2 | 10.2 | 9.5 | 21.1 | 10.2 |
| 样本7 | 9.5 | 21.1 | 12.0 | 10.2 | 12.0 |
| 样本8 | 10.2 | 18.2 | 10.2 | 21.1 | 21.1 |
| 样本9 | 10.2 | 10.2 | 18.2 | 18.2 | 18.2 |
| 样本10 | 18.2 | 10.2 | 18.2 | 10.2 | 10.2 |
对以上每个 bootstrap 样本,求得样本中位数分别为
$$ \hat{\theta}_{1}^{*}=12,0,\hat{\theta}_{2}^{*}=12,0,\hat{\theta}_{3}^{*}=10,2,\hat{\theta}_{4}^{*}=12,0,\hat{\theta}_{5}^{*}=18,2, $$
$$ \hat{\theta}_{6}^{*}=10,2,\hat{\theta}_{7}^{*}=12,0,\hat{\theta}_{8}^{*}=18,2,\hat{\theta}_{9}^{*}=18,2,\hat{\theta}_{10}^{*}=10,2, $$
于是以原始样本确定的样本中位数 $ \hat{\theta}=12.0 $作为总体中位数 $ \theta $的估计,由(1.2)式知其标准误差的bootstrap估计为
$$ \hat{\sigma}_{\hat{\theta}}=\sqrt{\frac{1}{9}\sum_{i=1}^{10}(\hat{\theta}_{i}^{*}-\overline{\theta}^{*})^{2}}=3.4579. $$
本题中取 B=10,这只是为了说明计算方法,是不能实际运用的,在实际中应取 $ B \geqslant 1000 $.
(二) 估计量的均方误差及偏差的 bootstrap 估计
设 $ \boldsymbol{X}=(X_{1},X_{2},\cdots,X_{n}) $ 是来自总体 F 的样本,F 未知, $ R=R(\boldsymbol{X}) $ 是感兴趣
先将原始样本的个体自左至右编号为1,2,3,4,5.需产生分布律为 $ P\{X=i\}=\frac{1}{5},i=1,2,3,4,5 $的5个随机数.为此先在随机数表中得到5个伪随机数
$$ \begin{array}{l} 0,21500\quad0,01011\quad0,47435\quad0,91312\quad0,12775 \end{array} $$
于是得分布 $ P\{X=i\}=\frac{1}{5}, i=1,2,3,4,5 $ 的 5 个随机数
$$ x_{1}=[5\times0.215\ 00]+1=2,\ x_{2}=[5\times0.010\ 11]+1=1,\ x_{3}=[5\times0.474\ 35]+1=3, $$
$$ x_{4}=[5\times0.913\ 12]+1=5,x_{5}=[5\times0.127\ 75]+1=1 $$
即2,1,3,5,1,于是得对应的本题的样本
$$ \begin{array}{l} 9.5 \quad 18.2 \quad 12.0 \quad 10.2 \quad 18.2 \end{array} $$
这就是这里的样本1.(参见第377页参读材料例2)
的随机变量,它依赖于样本X. 假设我们希望去估计R的分布的某些特征. 例如R的数学期望 $ E_{F}(R)^{\textcircled{1}} $,就可以按照上面所说的三个步骤 $ 1^{\circ}, 2^{\circ}, 3^{\circ} $进行,只是在 $ 2^{\circ} $中对于第i个bootstrap样本 $ x_{i}^{*}=(x_{1}^{*}, x_{2}^{*}, \cdots, x_{n}^{*}) $,计算 $ R_{i}^{*}=R(x_{i}^{*}) $代替计算 $ \theta_{i}^{*} $,且在 $ 3^{\circ} $中计算感兴趣的R的特征. 例如如果希望估计 $ E_{F}(R) $就计算
$$ E_{*}\left(R^{*}\right)\textcircled{2}=\frac{1}{B}\sum_{i=1}^{B}R_{i}^{*}\;. $$
例2(均方误差)设金属元素铂的升华热是具有分布函数 F 的连续型随机变量,F 的中位数 $ \theta $ 是未知参数,现测得以下的数据(以 kcal/mol 计 $ ^{③} $):
| 136.3 | 136.6 | 135.8 | 135.4 | 134.7 | 135.0 | 134.1 | 143.3 | 147.8 |
| 148.8 | 134.8 | 135.2 | 134.9 | 149.5 | 141.2 | 135.4 | 134.8 | 135.8 |
| 135.0 | 133.7 | 134.4 | 134.9 | 134.8 | 134.5 | 134.3 | 135.2 |
以样本中位数 $ M = M(\mathbf{X}) $ 作为总体中位数 $ \theta $ 的估计,试求均方误差 MSE = $ E\left[(M - \theta)^2\right] $ 的 bootstrap 估计.
解 将原始样本自小到大排序,左起第13个数为135.0,左起第14个数为135.2,于是样本中位数为 $ \frac{1}{2}(135.0+135.2)=135.1 $。以135.1作为总体中位数 $ \theta $的估计,即 $ \hat{\theta}=135.1 $。取 $ R=R(X)=(M-\hat{\theta})^{2} $,需估计 $ R(X) $的均值 $ E\left[(M-\hat{\theta})^{2}\right] $。
相继地、独立地抽取10 000个bootstrap样本如下:
样本1
| 133.2 | 134.1 | 134.1 | 134.1 | 134.8 | 134.8 | 134.8 | 134.9 | 134.9 |
| 134.9 | 135.0 | 135.2 | 135.2 | 135.4 | 135.4 | 135.8 | 135.8 | 136.3 |
| 136.3 | 136.6 | 136.6 | 141.2 | 143.3 | 143.3 | 147.8 | 148.8 |
得样本中位数为135.3
$$ \bullet \quad \bullet $$
样本1000
| 134.3 | 134.5 | 134.5 | 134.5 | 134.7 | 134.8 | 134.8 | 134.8 | 134.8 |
| 134.8 | 134.9 | 134.9 | 134.9 | 134.9 | 135.0 | 135.4 | 135.4 | 135.4 |
| 135.4 | 135.4 | 135.8 | 136.6 | 146.5 | 146.5 | 147.8 | 148.8 |
得样本中位数为134.9
对于用第i个样本计算
$$ R_{i}^{*}=R(x^{*_{i}})=(M_{i}^{*}-\hat{\theta})^{2}=(M_{i}^{*}-135.1)^{2},i=1,2,\cdots,10000. $$
即有对于样本1 $ (M_{1}^{*}-135.1)^{2}=(135.3-135.1)^{2}=0.04 $
$ \vdots $
对于样本 10 000 $ (M_{10 000}^{*}-135.1)^{2}=(134.9-135.1)^{2}=0.04 $.
用这10000个数的平均值
$$ \frac{1}{10\ 000}\sum_{i=1}^{10\ 000}(M_{i}^{*}\ -135.1)^{2}=0.07 $$
近似 $ E\left[(M-\theta)^{2}\right] $,即得 MSE $ \left[(M-\theta)^{2}\right] $ 的 bootstrap 估计为 0.07.
例3(偏差)设 $ \boldsymbol{X}=(X_{1},X_{2},\cdots,X_{n}) $ 是来自总体 F 的样本, $ \hat{\theta}=\hat{\theta}(X_{1},X_{2},\cdots,X_{n}) $ 是参数 $ \theta $ 的估计量。 $ \theta $ 的估计 $ \hat{\theta} $ 关于 $ \theta $ 的偏差定义为
$$ b=E(\hat{\theta}-\theta)=E(\hat{\theta})-\theta. $$
当 $ \theta $是 $ \theta $的无偏估计时b=0.
试在例2中,以样本中位数M=M(X)作为总体F的中位数 $ \theta $的估计,求偏差 $ b=E(M-\theta) $的bootstrap估计.
由例2知原始样本的中位数为135.1.以135.1作为总体中位数R= $ \theta $的估计,即 $ \hat{\theta}=135.1 $,取 $ R=R(X)=M-\hat{\theta} $,需估计 $ R(X) $的均值 $ E(M-\hat{\theta}) $。对于例2中第i个样本计算
$$ R_{i}^{*}=R(x^{*i})=(M_{i}^{*}-\hat{\theta})=(M_{i}^{*}-135.1),i=1,2,\cdots,10\ 000. $$
即有对于样本1 $ M_{1}^{*}-135.1=0.02 $.
$$ \vdots\quad\vdots $$
对于样本 10 000 $ M_{10\ 000}^{*}-135.1=-0.02 $.
将上述 10 000 个数取平均值得到偏差 b 的 bootstrap 估计为
$$ \begin{aligned}b^{*}=&\frac{1}{10\ 000}\sum_{i=1}^{10\ 000}(M_{i}^{*}-135.1)=\frac{1}{10\ 000}\sum_{i=1}^{10\ 000}M_{i}^{*}-135.1\\ =&135.14-135.1=0.04.\end{aligned} $$
(三) bootstrap 置信区间
下面介绍一种求未知参数 $ \theta $ 的 bootstrap 置信区间的方法.
设 $ \boldsymbol{X}=(X_{1},X_{2},\cdots,X_{n}) $ 是来自总体 F 容量为 n 的样本, $ \boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n}) $ 是一个已知的样本值. F 中含有未知参数 $ \theta,\hat{\theta}=\hat{\theta}(X_{1},X_{2},\cdots,X_{n}) $ 是 $ \theta $ 的估计量. 现在来求 $ \theta $ 的置信水平为 $ 1-\alpha $ 的置信区间.
相继地、独立地从样本 $ \boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n}) $ 中抽出 B 个容量为 n 的 bootstrap 样本,对于每个 bootstrap 样本求出 $ \theta $ 的 bootstrap 估计: $ \hat{\theta}_{1}^{*},\hat{\theta}_{2}^{*},\cdots,\hat{\theta}_{B}^{*} $. 将它们自小到大排序,得
$$ \hat{\theta}_{(1)}^{*}\leqslant\hat{\theta}_{(2)}^{*}\leqslant\cdots\leqslant\hat{\theta}_{(B)}^{*}. $$
取 $ R(\mathbf{X}) = \hat{\theta} $,用对应的 $ R(\mathbf{X}^*) = \hat{\theta}^* $ 的分布作为 $ R(\mathbf{X}) $ 的分布的近似,求出 $ R(\mathbf{X}^*) $ 的分布的近似分位数 $ \hat{\theta}_{a/2}^* $ 和 $ \hat{\theta}_{1-\alpha/2}^* $ 使
$$ P\{\hat{\theta}_{a/2}^{*}<\hat{\theta}^{*}<\hat{\theta}_{1-a/2}^{*}\}=1-\alpha, $$
于是近似地有
$$ \mathrm{P}\{\hat{\theta}_{\alpha/2}^{*}<\theta<\hat{\theta}_{1-\alpha/2}^{*}\}=1-\alpha. $$
记 $ k_{1}=\left[B\times\frac{\alpha}{2}\right],k_{2}=\left[B\times\left(1-\frac{\alpha}{2}\right)\right] $,在(1.5)式中以 $ \hat{\theta}_{(k_{1})}^{*} $ 和 $ \hat{\theta}_{(k_{2})}^{*} $ 分别作为分位数 $ \hat{\theta}_{a/2}^{*},\hat{\theta}_{1-a/2}^{*} $ 的估计,得到近似等式
$$ P\{\hat{\theta}_{(k_{1})}^{*}<\theta<\hat{\theta}_{(k_{2})}^{*}\}=1-\alpha. $$
于是由上式就得到 $ \theta $ 的置信水平为 $ 1-\alpha $ 的近似置信区间
$$ (\hat{\theta}_{(k_{1})}^{*},\hat{\theta}_{(k_{2})}^{*}) $$
这一区间称为 $ \theta $ 的置信水平为 $ 1-\alpha $ 的 bootstrap 置信区间。这种求置信区间的方法称为分位数法。
例4 在例2中
(1)以样本中位数作为总体中位数 $ \theta $ 的估计求 $ \theta $ 的置信水平为 0.95 的 bootstrap 置信区间;
(2)以样本20%截尾均值作为总体20%截尾均值 $ \mu_{t} $的估计,求 $ \mu_{t} $的置信水平为0.95的bootstrap置信区间.
解 n=26, B=10000,原始样本以及10000个模拟bootstrap样本见例2.
(1)对于每一个 bootstrap 样本算出中位数 $ M_{1}^{*}, M_{2}^{*}, \cdots, M_{10}^{*} \quad 000 $. 将它们自小到大排序得到
$$ \begin{aligned}M_{(1)}^{*}\leqslant&M_{(2)}^{*}\leqslant\cdots\leqslant M_{(250)}^{*}\leqslant M_{(251)}^{*}\leqslant\cdots\leqslant M_{(9\ 750)}^{*}\leqslant M_{(9\ 751)}^{*}\leqslant\cdots\leqslant M_{(10\ 000)}^{*}.\\B=&10\ 000,1-\alpha=0.95,\alpha=0.05,\end{aligned} $$
由
$$ k_{1}=\left[10\ 000\times\frac{0.05}{2}\right]=250,k_{2}=\left[10\ 000\times(1-\frac{0.05}{2})\right]=9\ 750. $$
置信区间为
$$ (M_{(250)}^{*},M_{(9750)}^{*})=(134.8,135.8). $$
(2)对于例2中的10 000个bootstrap样本中的每一个,算出样本20%截尾均值: $ \overline{x}_{t1}^{*} $, $ \overline{x}_{t2}^{*} $,…, $ \overline{x}_{t10 000}^{*} $,将它们自小到大排序得到
$$ \begin{aligned}\overline{x}_{t(1)}&\leqslant\overline{x}_{t(2)}^{\ast}\leqslant\cdots\leqslant\overline{x}_{t(250)}^{\ast}\leqslant\overline{x}_{t(251)}^{\ast}\leqslant\cdots\leqslant\overline{x}_{t(9750)}^{\ast}\\&\leqslant\overline{x}_{t(9751)}^{\ast}\leqslant\cdots\leqslant\overline{x}_{t(1000)}^{\ast}.\end{aligned} $$
按分位数法由(1.7)式得到20%截尾均值的一个置信水平为0.95的bootstrap
置信区间为
$$ (\overline{x}_{t(250)}^{*},\overline{x}_{t(9750)}^{*})=(134,85,136,92). $$
例 5 有 30 窝仔猪出生时各窝猪的存活只数为
| 9 | 8 | 10 | 12 | 11 | 12 | 7 | 9 | 11 | 8 | 9 | 7 | 7 | 8 | 9 | 7 |
| 9 | 9 | 10 | 9 | 9 | 9 | 12 | 10 | 10 | 9 | 13 | 11 | 13 | 9 |
以样本均值 $ \overline{x} $ 作为总体均值 $ \mu $ 的估计,以样本标准差 s 作为总体标准差 $ \sigma $ 的估计,按分位数法求 $ \mu $ 以及 $ \sigma $ 的置信水平为 0.90 的 bootstrap 置信区间.
解 相继地、独立地自原始样本数据用放回抽样的方法,得到10 000个容量均为30的bootstrap样本:
样本1
| 8 | 8 | 10 | 12 | 7 | 11 | 11 | 8 | 10 | 12 | 7 | 9 | 10 | 8 | 9 | 11 | 10 |
| 13 | 9 | 9 | 9 | 10 | 8 | 13 | 8 | 9 | 9 | 7 | 10 | 8 | ||||
| ⋮ | ⋮ | |||||||||||||||
| 样本 10 000 | ||||||||||||||||
| 9 | 10 | 7 | 10 | 9 | 7 | 9 | 7 | 10 | 7 | 9 | 9 | 13 | 11 | 12 | 10 | |
| 12 | 12 | 10 | 9 | 8 | 11 | 9 | 9 | 9 | 11 | 12 | 11 | 12 | 9 | |||
对上述每个 bootstrap 样本算出样本均值 $ \overline{x}_{i}^{*}(i=1,2,\cdots,10000) $,将 10 000 个 $ \overline{x}_{i}^{*} $ 按自小到大排序,左起第 500 位为 $ \overline{x}_{(500)}=9.03 $,左起第 9 500 位为 $ \overline{x}_{(9500)}=10.038 $。于是按 (1.7) 式得 $ \mu $ 的一个置信水平为 0.90 的 bootstrap 置信区间为
$$ (\overline{x}_{(500)}^{*},\quad\overline{x}_{(9500)}^{*})=(9,03,10,038). $$
对上述 10 000 个 bootstrap 样本的每一个算出标准差 $ s_{i}^{*} $ ( $ i=1,2,\cdots,10\,000 $),将 10 000 个 $ s_{i}^{*} $ 按自小到大排序。左起第 500 位为 $ s_{(500)}^{*}=1.35 $,左起第 9 500 位为 $ s_{(9\,500)}^{*}=1.98 $,于是按 (1.7) 式得 $ \sigma $ 的一个置信水平为 0.90 的 bootstrap 置信区间为
$$ (s_{(500)}^{*},~s_{(9~500)}^{*})=(1.35,~1.98). $$
(四)用 bootstrap-t 法求均值 $ \mu $ 的 bootstrap 的置信区间
设 $ \boldsymbol{X}=(X_{1},X_{2},\cdots,X_{n}) $ 是来自总体 F 的容量为 n 的样本, $ \boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n}) $ 是一个已知的样本值。均值 $ \mu $ 和方差 $ \sigma^{2} $ 均为未知参数,我们要利用样本值 x 来估计 $ \mu $。
考虑函数
$$ t=\frac{\bar{X}-\mu}{S/\sqrt{n}} $$
在第七章§5中假设总体F具有正态分布,此时t的分布与参数 $ \mu $无关,它是一个枢轴量而且有 $ t\sim t(n-1) $,利用枢轴量t,就能求得 $ \mu $的置信区间。现在,总体F不具有正态分布,但可证t仍是一个枢轴量。然而t的分布就不是 $ t(n-1) $分布,这样就不能按第七章的方法求得 $ \mu $的置信区间了,下面我们用bootstrap方法来求 $ \mu $的近似置信区间。
以原始样本 $ \boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n}) $ 的样本均值 $ \overline{x}=\frac{1}{n}\sum_{i=1}^{n}x_{i} $ 作为 $ \mu $ 的估计,考虑与 t 相应的枢轴量
$$ W^{*}=\frac{\overline{X}^{*}-\overline{x}}{S^{*}/\sqrt{n}}, $$
此处 $ \overline{X}^{*} $, $ S^{*} $分别为与 $ \overline{X} $,S相应的bootstrap样本均值与样本标准差。用 $ W^{*} $的分布近似t的分布,求出 $ W^{*} $的近似分位数 $ w_{a/2}^{*} $, $ w_{1-a/2}^{*} $,使
$$ P\{w_{a/2}^{*}<\frac{\overline{X}^{*}-\overline{x}}{S^{*}/\sqrt{n}} 于是近似地有 $$ \begin{align*}P\{w_{a/2}^{*}<\frac{X-\mu}{S/\sqrt{n}} 或即 将 $ W^{\star} $的B个bootstrap值自小到大排序 $$ \mathcal{W}_{(1)}^{*}\leqslant\mathcal{W}_{(2)}^{*}\leqslant\cdots\leqslant\mathcal{W}_{(B)}^{*}, $$ 记 $ k_{1}=\left[B\times\frac{\alpha}{2}\right],k_{2}=\left[B\times\left(1-\frac{\alpha}{2}\right)\right] $,在(1.11)式中以 $ w_{(k_{1})}^{*} $ 和 $ w_{(k_{2})}^{*} $ 分别作为分位数 $ w_{\alpha/2}^{*},w_{1-\alpha/2}^{*} $ 的估计,由(1.11)式得到近似等式 $$ P\{\overline{X}-w_{(k_{2})}^{*}\frac{S}{\sqrt{n}}<\mu<\overline{X}-w_{(k_{1})}^{*}\frac{S}{\sqrt{n}}\}=1-\alpha. $$ 由(1.12)式得到 $ \mu $ 的置信水平为 $ 1-\alpha $ 的 bootstrap 置信区间 $$ \bigl(\overline{{X}}-\mathcal{W}_{(k_{2})}^{*}\frac{S}{\sqrt{n}},\quad\overline{{X}}-\mathcal{W}_{(k_{1})}^{*}\frac{S}{\sqrt{n}}\bigr). $$ 这一方法称为 brotstrap-t 法. 例 6 在例 5 中用 bootstrap-t 法求 $ \mu $ 的置信水平为 0.90 的置信区间. 解 原始样本以及 10 000 个模拟 bootstrap 样本见例 5. 在原始样本中 n=30, $ \overline{x}=9.53 $, s=1.72, $ s^{2}=2.95 $. 对于第 i 个 $ (i=1,2,\cdots,10\ 000=B) $ bootstrap 样本,求出它的均值 $ \overline{x}_{i}^{*} $ 和 样本标准差 $ s_{i}^{*} $,从而得到 $ w^{*} $ 的第 i 个值 $$ w_{i}^{*}=\frac{\overline{x}_{i}^{*}-\overline{x}}{s_{i}^{*}/\sqrt{n}},\quad i=1,2,\cdots,1000. $$ 其中 $ \overline{x} $ 是由原始样本确定的样本均值. 将 $ w_{i}^{*} $ 自小到大排序得到 $$ \begin{array}{r}{\varpi_{(1)}^{*}\leqslant\varpi_{(2)}^{*}\leqslant\cdots\leqslant\varpi_{(1000)}^{*}.}\end{array} $$ 取置信水平 1 - $ \alpha = 0.90 $,此时 $ \alpha = 0.10 $, $ \alpha/2 = 0.05 $, $ 1 - \alpha/2 = 0.95 $,取 $ k_{1} = \left[ B \times \frac{\alpha}{2} \right] = 500 $, $ k_{2} = \left[ B \times \left( 1 - \frac{\alpha}{2} \right) \right] = 9500 $,得 $ w_{(500)}^{*} = -1.7813 $, $ w_{(9500)}^{*} = 1.6299 $,于是按 (1.13) 式得到 $ \mu $ 的置信水平为 0.90 的 bootstrap-t 置信区间为 $$ \left(9.53-1.629\ 9\times\frac{1.72}{\sqrt{30}},\quad9.53+1.781\ 3\times\frac{1.72}{\sqrt{30}}\right)=\quad(9.018\ 2,\quad10.089\ 4). $$ 用非参数 bootstrap 法来求参数的近似置信区间的优点是,不需要对总体分布的类型作任何的假设,而且可以适用于小样本,且能用于各种统计量(不限于样本均值)。 以上介绍的 bootstrap 方法,没有假设所研究的总体的分布函数 F 的形式,bootstrap 样本是来自已知的数据(原始样本),所以称之为非参数 bootstrap 方法. 假设所研究的总体的分布函数 $ F(x;\beta) $ 的形式已知,但其中包含未知参数 $ \beta $ ( $ \beta $ 可以是向量). 现在已知有一个来自 $ F(x;\beta) $ 的样本 $$ X_{1},X_{2},\cdots,X_{n}. $$ 利用这一样本求出 $ \beta $ 的最大似然估计 $ \hat{\beta} $. 在 $ F(x;\beta) $ 中以 $ \hat{\beta} $ 代替 $ \beta $ 得到 $ F(x;\hat{\beta}) $, 接着在 $ F(x;\hat{\beta}) $ 中产生容量为 n 的样本 $$ X_{1}^{*},X_{2}^{*},\cdots,X_{n}^{*}\sim F(x;\hat{\beta})^{\textcircled{1}}. $$ 这种样本可以产生很多个,例如产生 B 个 $ (B \geqslant 1000) $,就可以利用这些样本对总体进行统计推断,其做法与非参数 bootstrap 方法一样。这种方法称为参数 bootstrap 法。 例 1 已知某种电子元件的寿命(以 h 计)服从韦布尔分布,其分布函数为 $$ F(x)=\{\begin{aligned}&1-e^{-(x/\eta)^{\beta}},&x>0,\\ &0,& 其他 ,\end{aligned}.\quad\beta>0,\eta>0. $$ 概率密度为 $$ f(x)=\{\begin{aligned}&\frac{\beta}{\eta^{\beta}}x^{\beta-1}e^{-(x/\eta)^{\beta}},&x>0,\\ &0,& 其他 ,\end{aligned}. $$ 已知参数 $ \beta=2 $ 。今有样本 142.84 97.04 32.46 69.14 85.67 114.43 41.76 163.07 108.22 63.28 (1)确定参数 $ \eta $的最大似然估计. (2)对于时刻 $ t_{0}=50 $ ,求可靠性 $ R(50)=1-F(50)=\mathrm{e}^{-(50/\eta)^{2}} $ 的置信水平分别为 0.95,0.90 的 bootstrap 单侧置信下限. 解 (1)设有样本 $ x_{1}, x_{2}, \cdots, x_{n} $,似然函数为(已将 $ \beta = 2 $ 代入) $$ L=\prod_{i=1}^{n}\frac{2}{\eta^{2}}x_{i}\mathrm{e}^{-(x_{i}/\eta)^{2}}=\frac{2^{n}}{\eta^{2n}}\left(\prod_{i=1}^{n}x_{i}\right)\mathrm{e}^{-(\sum_{i=1}^{n}x_{i}^{2}/\eta^{2}}, $$ $$ \ln L=C+\left[-2n\ln\eta-\frac{1}{\eta^{2}}\sum_{i=1}^{n}x_{i}^{2}\right],\quad C 为常数 . $$ 令 $ \frac{d}{d\eta}\ln L=0 $得 $$ \frac{-2n}{\eta}+\frac{2}{\eta^{3}}\sum_{i=1}^{n}x_{i}^{2}=0, $$ $$ \hat{\eta}=\sqrt{\frac{\sum_{i=1}^{n}x_{i}^{2}}{n}}. $$ 以数据代入得 $ \eta $ 的最大似然估计为 $ \hat{\eta}=100.0696 $. (2)先说一下如何产生韦布尔分布的随机数. 设 $ U \sim U(0,1) $,令 U = 1 - $ \mathrm{e}^{-(X/\eta)^{2}} $,解得 $$ X=\eta[-\ln(1-U)]^{1/2}. $$ 因 $ 1-U \sim U(0,1) $,故 $$ X=\eta[-\ln U]^{1/2} $$ 也具有参数 $ \beta=2 $, $ \eta $ 的韦布尔分布。以 $ \hat{\eta}=100.069 $ 6 作为 $ \eta $,按 $ X=100.069\ 6[-\ln U]^{1/2} $ 就能产生韦布尔分布的随机数。以 $ F(x,\hat{\eta})=F(x,100.069\ 6) $ 为分布函数产生 5 000 个容量为 10 的 bootstrap 样本: 样本1 $ x_{1}^{*1}, x_{2}^{*1}, \cdots, x_{10}^{*1} $,得 $ \eta $的bootstrap估计 $$ \eta_{1}^{*}=\sqrt{\frac{\sum_{i=1}^{10}(x_{i}^{*1})^{2}}{10}}. $$ $$ \bullet \quad \bullet $$ 样本 5 000 $ x_{1}^{*5000}, x_{2}^{*5000}, \cdots, x_{10}^{*5000} $,得 $ \eta $ 的 bootstrap 估计 $$ \eta_{5000}^{*}=\sqrt{\frac{\sum_{i=1}^{10}(x_{i}^{*5000})^{2}}{10}}. $$ 将以上 5 000 个 $ \eta_{i}^{*} $ 自小到大排序,取左起第 250 位,得 $$ \hat{\eta}_{(250)}^{*}=73.257\ 36, $$ 取左起第500位得 $$ \hat{\eta}_{(500)}^{*}=79.036\ 52. $$ 于是在 t=50 时,可靠性 R(50) 的置信水平为 0.95 的 bootstrap 单侧置信下限为 $$ \mathrm{e}^{-(50/\bar{\eta}_{(250)}^{*})^{2}}=0.6276. $$ 在 t=50 时,可靠性 R(50) 的置信水平为 0.90 的 bootstrap 单侧置信下限为 $$ \mathrm{e}^{-(50/\hat{\eta}_{(500)}^{*})^{2}}=0.6702. $$ 例2 据 Hardy−Weinberg 定律,若基因频率处于平衡状态,则在一总体中个体具有血型 M、MN、N 的概率分别是 $ (1-\theta)^{2}, 2\theta(1-\theta), \theta^{2} $,其中 $ 0<\theta<1 $。据1937年对香港地区的调查有以下的数据: (1) 求 $ \theta $ 的最大似然估计 $ \hat{\theta} $; (2)求 $ \theta $的置信水平为0.90的bootstrap置信区间. 解 分别记 $ x_{1}, x_{2}, x_{3} $ 为具有血型为 M, MN, N 的人数,记 $ x_{1} + x_{2} + x_{3} = n $. 似然函数为 $$ \begin{aligned}L=&\left[(1-\theta)^{2}\right]^{x_{1}}\left[2\theta(1-\theta)\right]^{x_{2}}\left[\theta^{2}\right]^{x_{3}}\\=&2^{x_{2}}\theta^{x_{2}+2x_{3}}(1-\theta)^{2x_{1}+x_{2}},\end{aligned} $$ $$ \ln L=x_{2}\ln2+(x_{2}+2x_{3})\ln\theta+(2x_{1}+x_{2})\ln(1-\theta). $$ 令 $ \frac{\mathrm{d}}{\mathrm{d}\theta}\ln L=\frac{x_{2}+2x_{3}}{\theta}+\frac{-(2x_{1}+x_{2})}{1-\theta}=0 $ 解得 $$ \hat{\theta}=\frac{x_{2}+2x_{3}}{2x_{1}+2x_{2}+2x_{3}}=\frac{x_{2}+2x_{3}}{2n} $$ 以数据 $ x_{1}=342, x_{2}=500, x_{3}=187, n=1029 $ 代入得到 $ \hat{\theta}=0.4247 $。以 $ \hat{\theta} $ 代替 $ \theta $,得到 $ (1-\theta)^{2}=0.331, 2\theta(1-\theta)=0.489, \theta^{2}=0.180 $。于是血型的近似分布律为 以(2.1)为分布律产生1 000个bootstrap样本,从而得到 $ \theta $的1 000个bootstrap估计 $ \hat{\theta}_{1}^{*} $, $ \hat{\theta}_{2}^{*} $,…, $ \hat{\theta}_{1}^{*} $。将这1 000个数按自小到大的次序排序得到 $$ \hat{\theta}_{(1)}^{*}\leqslant\hat{\theta}_{(2)}^{*}\leqslant\cdots\hat{\theta}_{(50)}^{*}\leqslant\cdots\leqslant\hat{\theta}_{(950)}^{*}\leqslant\cdots\leqslant\hat{\theta}_{(1000)}^{*}. $$ 取 $ \left(\hat{\theta}_{(50)}^{*},\hat{\theta}_{(950)}^{*}\right)=(0.378\ 5,0.409\ 6) $为 $ \theta $的置信水平为0.90的bootstrap置信区间. 设 $ x=(x_1,x_2,\cdots,x_n) $ 是来自分布函数为 F 的总体的样本,F 未知。 $ R(x) $ 是 x 的函数, $ F_n $ 是相应的经验分布函数。假如我们感兴趣的是 $ R(x) $ 的某些特征,例如 R 的均值或中位数。非参数 bootstrap 方法的第一步是用已知的经验分布函数 $ F_n $ 代替 F,在 $ F_n $ 中抽样,得到数据样本 $ x^*=(x_1^*,x_2^*,...,x_n^*) $,然后计算 $ R(x^*) $ 的均值或中位数,作为所需求的均值或中位数的估计(bootstrap 估计)。通常的情况, $ R(x^*) $ 的分布过于复杂,不能用解析的方法计算得到 $ R(x^*) $ 的特征,而需要采用模拟的方法。在参数 bootstrap 方法中 $ F=F(x;\beta) $ 的形式已知,但包含未知参数 $ \beta $。先利用样本 x 求出 $ \beta $ 的最大似然估计 $ \hat{\beta} $,以 $ F(x;\hat{\beta}) $ 代替 F,在 $ F(x;\hat{\beta}) $ 中抽样得到数据样本 $ x^* $,然后计算 $ R(x^*) $ 的均值或中位数,作为所需求的均值或中位数的 bootstrap 估计。非参数和参数 bootstrap 方法可用于当人们对总体知之甚少的情况,它们是近代统计中的一种用于数据处理的重要实用方法。§2 参数 bootstrap 方法
血型 M MN N 人数 342 500 187 共 1 029 血型 M MN N 概率 0.331 0.489 0.180 小结