第十一章 在数理统计中应用 Excel 软件
第十一章 在数理统计中应用 Excel 软件
§1 概述
(一) 计算机技术在数理统计中的应用
随着现代科学技术的迅猛发展,人类社会已开始进入一个利用和开发信息资源的信息社会。信息数据数量大、范围广、变化快,传统的人工处理手段无法适应社会。经济高速发展对统计提出的要求,也难以提高数据处理的速度和精度。计算机技术在数理统计中的应用,主要是在统计信息的存贮和检索、统计资料的分析和检验等方面的应用,解决了统计工作中的难题。
不仅是在实际的技术和经济工作中要将计算机技术应用于数理统计,在学习概率论与数理统计课程的阶段,同样也需要应用计算机技术。掌握了计算机技术在数理统计中的应用以后,读者的分析和研究问题的能力将极大地提高,研究问题的规模、分析计算的效率将极大地提高。
(二) 在数理统计研究中应用 Excel 软件
功能强大的统计分析软件有 SAS(Statistical Analysis Software)、SPSS(原名为 Statistical Package for the Social Science,2000 年改为 Statistical Product and Service Solutions)等等,但是所有这些专业软件往往系统庞大、结构复杂,大多数非统计专业人员难以运用自如,而且价格昂贵,是一般人难以承受的。
微软(Microsoft)公司推出的办公软件包 Office,得到了广泛的应用,Excel 是 Office 的重要成员之一。Excel 是一个功能多、技术先进、使用方便的表格式数据综合管理和分析系统,它采用电子表格方式进行数据处理,工作直观方便;它提供了丰富的函数,可以进行数据处理、统计分析和决策辅助;还具有较好的制图功能。
只要读者使用的计算机上安装了 Office,随之就有了 Excel,不需要另加投资,而 Excel 的使用是不难学会的。
Office 有不同的版本,例如 Office 2000 和 Office XP. Excel 也有不同的版本,例如 Excel 2000,Excel 2002 和 Excel 2003. 这些版本大同小异,相互之间兼容性好.
本书中,应用 Excel 处理数理统计问题。
计算机开机后,单击显示屏左下角的“开始”按钮,然后单击“所有程序”按钮,若弹出的菜单有“Microsoft Excel”,即可调用。
启动 Excel 后就会打开 Excel 的用户界面窗口,如图 11-1 所示。该窗口自上而下有标题栏、菜单栏、常用工具栏、格式工具栏、编辑栏、工作表区、工作表标签、水平滚动条和状态栏。

工作表区由单元格组成,每个单元格由列标和行号标识。工作表区的最上面一行为列标,用 A,…,Z,AA,…,AZ,BA,…BZ,…,IA,…,IV 表示,最多可使用 256 列。工作表区左边一列为行号,用 1,2,…,65 536 表示,最多可使用 65 536 行。单元格“A1”表示单元格位于 A 列第 1 行。单元格区域则规定为矩形,例如,“A1:F5”表示一矩形区域,A1 和 F5 为其主对角线两端的单元格。每张工作表有一个标签与之对应,例如,“sheet 1”。工作表隶属于工作簿,一个工作簿最多可由 255 个不同的工作表组成。
(三) Excel 的分析工具库
检查 Excel 的“工具”菜单,看是否已安装了分析工具。如果在“工具”菜单中没有“数据分析”项,则需调用“加载宏”来安装“分析工具库”。
“工具”菜单中有了“数据分析”命令项,单击它,就出现“数据分析”对话框,其中有19个模块,它们分别属于五类:
- 基础分析:(1)随机数发生器;(2)抽样;(3)描述统计;(4)直方图;(5)排位与百分比排位.
- 检验分析: (6) t 检验, 平均值的成对两样本分析; (7) t 检验, 双样本等方差假设; (8) t 检验, 双样本异方差假设; (9) Z 检验, 双样本平均差检验; (10) F 检验, 双样本方差.
- 相关、回归:(11) 相关系数;(12) 协方差;(13) 回归.
4.方差分析:(14)方差分析,单因素方差分析;(15)方差分析,可重复双因素分析;(16)方差分析,无重复双因素分析.
- 其他分析工具: (17) 移动平均; (18) 指数平滑; (19) 傅里叶分析.
在本书中,只讲述 Excel 在几个问题上的应用。
§2 箱线图
利用 Excel 函数
$$ Q U A R T I L E(a r r a y,q u a r t i l e) $$
就能直接得到箱线图的5个点:MIN, $ Q_{1} $, $ Q_{2} $(中位数), $ Q_{3} $和MAX.函数QUARTILE(array,quartile)的功能是返回数据集的四分位数.这一函数有两个参数array和quartile,参数array是需要求四分位数的数组或数字型单元格区域;参数quartile取0或1或2或3或4,依次表示数据集的MIN, $ Q_{1} $, $ Q_{2} $(中位数), $ Q_{3} $,MAX $ ^{①} $.例如array取A1:A10,quartile取3,则函数
$$ Q U A R T I L E(A1:A10,3) $$
表示返回数组 A1~A10 的 $ Q_{3} $. 函数
$$ Q U A R T I L E(A1:A10,4) $$
表示返回数组 A1~A10 的 MAX.
例 1 分别画出第六章§2 例3 中女子、男子组肺活量的箱线图.
解 打开 Excel 工作表,将“女子组”以及“男子组”分别键入单元格 A1 和 E1 以及 B1 和 G1,将数据分别键入单元格区域 A2:B26。将“MIN, $ Q_{1} $, $ Q_{2} $, $ Q_{3} $, MAX”分别键入单元格 D2:D6。在单元格 E2, E3, …, E6, G2, G3, …, G6,分别键入有关 Excel 函数式如下。
$$ \mathrm{E2}=\mathrm{QUARTILE}\left(\mathrm{A2}:\mathrm{A26,0}\right),\quad\mathrm{G2}=\mathrm{QUARTILE}\left(\mathrm{B2}:\mathrm{B26,0}\right), $$
$$ \mathrm{E3}=\mathrm{QUARTILE}\left(\mathrm{A2}:\mathrm{A26,1}\right),\quad\mathrm{G3}=\mathrm{QUARTILE}\left(\mathrm{B2}:\mathrm{B26,1}\right), $$
$$ \mathrm{E4}=\mathrm{QUARTILE}\left(\mathrm{A2}:\mathrm{A26,2}\right),\quad\mathrm{G4}=\mathrm{QUARTILE}\left(\mathrm{B2}:\mathrm{B26,2}\right), $$
$$ \mathrm{E5}=\mathrm{QUARTILE}\left(\mathrm{A2}:\mathrm{A26,3}\right),\quad\mathrm{G5}=\mathrm{QUARTILE}\left(\mathrm{B2}:\mathrm{B26,3}\right), $$
$$ \mathrm{E6}=\mathrm{QUARTILE}\left(\mathrm{A2}:\mathrm{A26,4}\right),\quad\mathrm{G6}=\mathrm{QUARTILE}\left(\mathrm{B2}:\mathrm{B26,4}\right), $$
即得所需结果如图 11-2 中的表格所示.
由所得数据即可由手工画出箱线图如图6-4所示.
我们也可以借助 Excel 的绘图工具栏,使用鼠标作出箱线图。
以女子的肺活量数据作箱线图为例. 先在草稿纸上画一草图, 考虑将箱线图的“线”放在

E列和F列之间的分隔线上.在这根线上确定箱线图的5个点的位置:将第8行下面的线与EF列分隔线相交的点定为2.7,这就是MIN的位置.每一行的高度定为0.2,于是, $ Q_{1},Q_{2},Q_{3} $和MAX的位置就相应地确定.
在绘图工具栏中单击“矩形”图标,将鼠标箭头放到矩形拟放的位置单击左键,即产生一个矩形。在矩形的四个角点和四条边的中点共有8个控点,单击矩形的边线中间的控点(此时出现两箭头的光标)就可以在纵或横方向缩放,单击矩形角点上的控点(此时出现两箭头的斜光标)可以在纵横向同步缩放;单击矩形内部的点(此时出现四箭头的光标)可以移动矩形的位置。经调整,矩形的下边和上边分别落在拟放的位置上 $ ^{①} $。矩形的宽度不拘,应使矩形关于EF列分隔线对称。
在绘图工具栏中单击直线图标,将光标指向要画的直线的起点,单击左键,按住鼠标拖到直线终点,松开,就画成了一条直线。画上 $ Q_{2} $ 线, $ Q_{0}-Q_{1} $ 线和 $ Q_{3}-Q_{4} $ 线,箱线图就完成了。
§3 假设检验
(一) 假设检验问题 p 值的求法
例1 求第八章§2例1的检验问题的p值.这一问题使用的是t检验法.样本容量为n=16,是单边检验.已由样本得到检验统计量 $ t=\frac{\overline{X}-\mu_{0}}{s/\sqrt{n}} $的观察值为
$$ t_{0}=\frac{241.5-225}{98.725\ 9/\sqrt{16}}=0.668\ 5. $$
下面用 Excel 来求 p 值。打开一个 Excel 工作表,单击“插入”,对于弹出的菜单单击“函数”,对接着弹出的菜单单击“统计”,再对弹出的菜单单击“TDIST”,然后单击“确定”,即弹出对话框。
对于显示的对话框,键入 x=0.6685, df(自由度)=15, Tails(尾数)=1 $ ^{①} $,即得 TDIST(0.6685,15,1)=0.256985562,这就是所需求的 p 值. ☐
例2 求第八章§3例1的检验问题的p值.这一问题使用的是 $ \chi^{2} $检验法.样本容量为n=26,是双边检验.已由样本得到检验统计量 $ \chi^{2}=\frac{(n-1)S^{2}}{\sigma_{0}^{2}} $的观察值为
$$ \chi_{0}^{2}=\frac{25\times9\ 200}{5\ 000}=46. $$
现用 Excel 来求 p 值。打开一个 Excel 工作表,单击“插入”,对于弹出的菜单单击“函数”,对弹出的菜单单击“统计”,再对弹出的菜单单击“CHIDIST”,然后单击“确定”。即弹出对话框。
对显示的对话框键入 x=46, df(自由度)=25, 即得 CHIDIST(46,25)=0.006417833, 这是单边检验的 p 值. 本题是双边检验, 故所求的 p 值 =2\times0.006417833=0.012835666.
例 3 求第八章习题第 19 题的检验问题的 p 值. 这一问题使用的是 F 检验法,是单边检验. 已由样本得到统计量 $ F = s_{1}^{2} / s_{2}^{2} $ 的观察值为
$$ F_{0}=s_{1}^{2}/s_{2}^{2}=1.60 $$
可得 p 值为 FDIST(1.60,59,39)=0.060 518 8.
(二)两个等方差正态总体 $ N(\mu_{1}, \sigma^{2}) $, $ N(\mu_{2}, \sigma^{2}) $ 均值差的检验(t 检验)
设 $ X_{1}, X_{2}, \cdots, X_{n_{1}} $ 是来自正态总体 $ N(\mu_{1}, \sigma^{2}) $ 的样本, $ Y_{1}, Y_{2}, \cdots, Y_{n_{2}} $ 是来自正态总体 $ N(\mu_{2}, \sigma^{2}) $ 的样本,两样本独立。 $ \mu_{1}, \mu_{2}, \sigma^{2} $ 均未知,现用 Excel 来求解假设检验问题
$$ \begin{aligned}&H_{0}:\mu_{1}{\leqslant}\mu_{2},\quad H_{1}:\mu_{1}>\mu_{2};\\&H_{0}:\mu_{1}{\geqslant}\mu_{2},\quad H_{1}:\mu_{1}<_{\mu_{2}};\\&H_{0}:\mu_{1}=\mu_{2},\quad H_{1}:\mu_{1}{\neq}\mu_{2}.\\ \end{aligned} $$
举例来说明.
例4 在两批电阻器中分别随机地取6只,测得以下的电阻值(以Ω计)
| A 批(x) | 0.140 | 0.138 | 0.143 | 0.142 | 0.144 | 0.137 |
| B 批(y) | 0.135 | 0.140 | 0.142 | 0.136 | 0.138 | 0.140 |
设两批电阻器电阻分别来自总体 $ N(\mu_{1}, \sigma^{2}) $, $ N(\mu_{2}, \sigma^{2}) $, $ \mu_{1}, \mu_{2} $, $ \sigma^{2} $ 均未知。两样本独立。试取 $ \alpha = 0.05 $,检验假设
$$ H_{0}:\mu_{1}=\mu_{2}\;,\quad H_{1}:\mu_{1}\neq\mu_{2}. $$
解 用 Excel 求解的操作步骤如下.
1° 打开 Excel 工作表,将数据输入单元格 A1:A7 和 B1:B7.
2°依次单击“工具”,“数据分析”,“t一检验:双样本等方差假设”和“确定”,跳出对话框.
$ 3^{\circ} $ 在对话框中键入变量1的范围 A1:A7,键入变量2的范围B1:B7;在假定均值差空格中键入“0”;单击“标志”,确认 $ \alpha=0.05^{①} $,单击“确定”,跳出一个新工作表如下:
| A(x) | B(y) | |
| 平均 | 0.140 667 | 0.138 5 |
| 方差 | 7.87E-06 | 7.1E-06 |
| 观测值 | 6 | 6 |
| 合并方差 | 7.48E-06 | |
| 假设平均差 | 0 | |
| df | 10 | |
| t Stat | 1.371 845 | |
| P(T<=t)单 | 0.100 051 | |
| t 单尾临界 | 1.812 461 | |
| P(T<=t)双 | 0.200 102 | |
| t 双尾临界 | 2.228 139 |
$ 4^{\circ} $ 结果分析 可以用两种方法来判别检验的结果.
(1)临界值法 这是双边检验. 因此需将“t Stat”(t 统计量)与“t 双边临界值”的大小进行比较. 现在 t 统计量的值为 1.371845, 它小于 t 双边临界值 2.228139. 故以 $ \alpha=0.05 $ 的显著性水平接受 $ H_{0} $.
(2)p 值法 由于双边检验的 p 值为 0.200 102,大于 0.05,故接受 $ H_{0} $. ☐
§4 方差分析
(一)单因素方差分析
例 1 在 7 个不同实验室中测量某种扑尔敏药片的扑尔敏有效含量(以 mg
计). 得到以下的结果(Lab 表示实验室):
| Lab 1 | Lab 2 | Lab 3 | Lab 4 | Lab 5 | Lab 6 | Lab 7 |
| 4.13 | 3.86 | 4.00 | 3.88 | 4.02 | 4.02 | 4.00 |
| 4.07 | 3.85 | 4.02 | 3.88 | 3.95 | 3.86 | 4.02 |
| 4.04 | 4.08 | 4.01 | 3.91 | 4.02 | 3.96 | 4.03 |
| 4.07 | 4.11 | 4.01 | 3.95 | 3.89 | 3.97 | 4.04 |
| 4.05 | 4.08 | 4.04 | 3.92 | 3.91 | 4.00 | 4.10 |
| 4.04 | 4.01 | 3.99 | 3.97 | 4.01 | 3.82 | 3.81 |
| 4.02 | 4.02 | 4.03 | 3.92 | 3.89 | 3.98 | 3.91 |
| 4.06 | 4.04 | 3.97 | 3.90 | 3.89 | 3.99 | 3.96 |
| 4.10 | 3.97 | 3.98 | 3.97 | 3.99 | 4.02 | 4.05 |
| 4.04 | 3.95 | 3.98 | 3.90 | 4.00 | 3.93 | 4.06 |
设各样本分别来自正态总体 $ N(\mu_{i}, \sigma^{2}) $,i=1,2,\cdots,7,各样本相互独立。试取显著性水平 $ \alpha=0.05 $ 检验各实验室测量的扑尔敏的有效含量的均值是否有显著差异。
解 我们画出各实验室测量结果的箱线图如图11-3所示. 从图上可以看出各实验室的测量结果的大致情况.

下面用 Excel 来求解. 具体步骤如下:
$ 1^{\circ} $ 建立给定问题的原假设和备择假设
$$ H_{0}:\mu_{1}=\mu_{2}=\cdots=\mu_{7},\quad H_{1}:\mu_{1},\mu_{2},\cdots,\mu_{7} 不全相等 . $$
2° 打开 Excel 工作表,将数据输入 A1:G11.
3°依次单击“工具”,“数据分析”,“方差分析:单因素方差分析”和“确定”,跳出对话框.
$ 4^{\circ} $ 在对话框中键入变量的输入范围“A1:G11”,单击“标志位于第一行”。规定 $ \alpha=0.05 $,单击“确定”,显示结果。有两张表,前一张是各实验室的均值、方差等的汇总,后一张是本题的方差分析表:
| 差异源 | SS | df | MS | F | P-value | F crit |
| 组间 | 0.124 737 | 6 | 0.020 79 | 5.660 069 | 9.45E-05 | 2.246 408 |
| 组内 | 0.231 4 | 63 | 0.003 673 | |||
| 总计 | 0.356 137 | 69 |
5° 结果分析 临界值法. F=5.660 069 大于 F crit(即 F 临界值)=2.246 408,故拒绝 $ H_{0} $,认为各实验室测量的结果有显著差异.
p 值法. p 值 = 9.45E-0.5 远小于 $ \alpha = 0.05 $. 故拒绝 $ H_{0} $,且知差异是非常显著的.
(二) 双因素无重复试验的方差分析
我们用 Excel 来求解第九章 §2 例 3 的双因素无重复试验的方差分析问题.
做法如下:
$ 1^{\circ} $ 建立给定问题的原假设和备择假设.
按题意需在显著性水平 $ \alpha=0.05 $下检验:在不同时间下颗粒状物含量的均值有无显著差异,在不同地点下颗粒状物含量的均值有无显著差异,即需检验假设(见第九章§2(2.25),(2.26)式).
$$ \begin{aligned}&H_{01}:\alpha_{1}=\alpha_{2}=\cdots=\alpha_{r}=0,\quad\\ &H_{11}:\alpha_{1},\alpha_{2},\cdots,\alpha_{r} 不全为零 .\\&H_{02}:\beta_{1}=\beta_{2}=\cdots=\beta_{s}=0,\quad\\ &H_{12}:\beta_{1},\beta_{2},\cdots,\beta_{s} 不全为零 .\\ \end{aligned} $$
2° 打开 Excel 工作表. 将数据输入 A1:F5.
3°依次单击“工具”,“数据分析”,“方差分析:无重复双因素分析”和“确定”,跳出对话框.
$ 4^{\circ} $ 在对话框中键入变量的输入范围“A1:F5”,单击“标志”,规定 $ \alpha=0.05 $,
单击“确定”. 显示结果有两张表,后面的一张是本题的方差分析表,如下所示:
| 差异源 | SS | df | MS | F | P-value | F crit |
| 行 | 1 182.95 | 3 | 394.316 7 | 10.722 41 | 0.001 033 | 3.490 295 |
| 列 | 1 947.5 | 4 | 486.875 | 13.239 29 | 0.000 234 | 3.259 167 |
| 误差 | 441.3 | 12 | 36.775 | |||
| 总计 | 3 571.75 | 19 |
$ 5^{\circ} $ 结果分析 临界值法. F=10.722 41 大于 F crit=3.490 295; F=13.239 29 大于 F crit=3.259 167, 故拒绝 $ H_{01} $ 及 $ H_{02} $. 即认为不同时间下颗粒状物含量的均值有显著差异;认为不同地点下颗粒状物含量的均值有显著差异. 也可以用 p 值来判别. □
(三) 双因素等重复试验的方差分析
现在用 Excel 来求解第九章 §2 例 2 的双因素等重复试验的方差分析问题.
做法如下:
$ 1^{\circ} $ 建立给定问题的原假设和备择假设.
按题意需在显著性水平 $ \alpha=0.05 $ 下检验. 热处理温度、时间以及这两者的交互作用对产品强度是否有显著的影响,即需检验假设(见第九章(2.6),(2.7),(2.8)式)
$$ \begin{aligned}&\{\begin{aligned}H_{01}&:\alpha_{1}=\alpha_{2}=\cdots=\alpha_{r}=0,\\H_{11}&:\alpha_{1},\alpha_{2},\cdots,\alpha_{r} 不全为零 .\end{aligned}.\\&\{\begin{aligned}H_{02}&:\beta_{1}=\beta_{2}=\cdots=\beta_{s}=0,\\H_{12}&:\beta_{1},\beta_{2},\cdots,\beta_{s} 不全为零 .\end{aligned}.\\&\{\begin{aligned}H_{03}&:\gamma_{11}=\gamma_{12}=\cdots=\gamma_{rs}=0,\\H_{13}&:\gamma_{11},\gamma_{12},\cdots,\gamma_{rs} 不全为零 .\end{aligned}.\end{aligned} $$
2° 打开 Excel 工作表,将数据输入 A1:C5.
3°依次单击“工具”,“数据分析”,“方差分析:可重复双因素方差分析”和“确定”,跳出对话框.
$ 4^{\circ} $ 在对话框中键入变量的输入范围“A1:C5”,再键入“每一样本的行数”为“2”。规定 $ \alpha=0.05 $,单击“确定”,即显示本题的方差分析表如下:
| 差异源 | SS | df | MS | F | P-value | F crit |
| 样本 | 1.62 | 1 | 1.62 | 1.408696 | 0.300945 | 7.708647 |
| 列 | 11.52 | 1 | 11.52 | 10.01739 | 0.03402 | 7.708647 |
| 交互 | 54.08 | 1 | 54.08 | 47.02609 | 0.002367 | 7.708647 |
| 内部 | 4.6 | 4 | 1.15 | |||
| 总计 | 71.82 | 7 |
5 $ ^{~} $ 结果分析 由于 p 值=0.300 945 大于 $ \alpha=0.05 $,故接受 $ H_{01} $,认为时间对强度的影响不显著;由于 p 值=0.034 02 小于 $ \alpha=0.05 $,故拒绝 $ H_{02} $,认为温度对强度影响显著;由于 p 值=0.002 367 小于 $ \alpha=0.05 $,故拒绝 $ H_{03} $,认为交互作用的影响显著. □
§5 一元线性回归
我们以例题来说明用 Excel 求解一元线性回归问题的做法.
例1 将冰晶放入一容器内,容器内维持固定的温度 $ (-5\,^{\circ}\mathrm{C}) $和固定的湿度. 观察自冰晶放入的时刻开始计算的时间 T(以s计)和晶体生长的轴向长度 A(以 $ \mu $m计),得到43对观察数据如下.
| T | 50 | 60 | 60 | 70 | 70 | 80 | 80 | 90 | 90 | 90 | 95 | 100 | 100 | 100 | 105 | 105 |
| A | 19 | 20 | 21 | 17 | 22 | 25 | 28 | 21 | 25 | 31 | 25 | 30 | 29 | 33 | 35 | 32 |
| T | 110 | 110 | 110 | 115 | 115 | 115 | 120 | 120 | 120 | 120 | 125 | 130 | 130 | 135 | 135 | 135 |
| A | 30 | 28 | 30 | 31 | 36 | 30 | 36 | 25 | 28 | 28 | 31 | 32 | 34 | 35 | 36 | 35 |
| T | 140 | 140 | 145 | 150 | 150 | 155 | 155 | 160 | 160 | 160 | 160 | 165 | 170 | 170 | 180 | 180 |
| A | 26 | 33 | 31 | 36 | 33 | 41 | 33 | 40 | 30 | 37 | 32 | 35 | 38 | 39 | 40 | 38 |
设题目符合回归模型所要求的条件.
(1)画出散点图;
(2) 求线性回归方程 $ \hat{A} = \hat{a} + \hat{b}T $;
(3)检验假设 $ H_{0}: b = 0, H_{1}: b \neq 0 $.
解 $ 1^{\circ} $ 打开 Excel 工作表,将数据输入单元格 A1:A44 和 B1:B44.
2°依次单击“插入”,“图表”,“XY散点图”,“下一步”,弹出对话框.在框中“数据区域”键入“A1:B44”,在“系列产生在”认定“列”。单击“下一步”,弹出“图表选项”对话框,在这一对话框的“图表标题”键入“冰晶长度一时间”,在“X
轴”键入“A2:A44”,在“Y轴”键入“B2:B44”。单击“完成”,即显示散点图,如图11-4所示。

3° 重新打开 Excel 工作表,将数据输入单元格 A1:B44. 依次单击“工具”,“数据分析”,“回归”和“确定”。弹出对话框。在“Y值输入区域”框键入“B1:B44”,在“X值输入区域”框键入“A1:A44”,单击“标志”,认定置信水平为 95%,“输出选项”选定“新工作表组”,单击“确定”,即得计算结果表格。输出的表格共三张,本题的结果载于最后一张表格上,如下所示:
| Coefficient | 标准误差 | t Stat | P-value | Lower 95% | Upper 95% | 下限 95% | 上限 95% | |
| Intercept | 14.41065 | 2.149358 | 6.70463 | 4.31E-08 | 10.06993 | 18.75136 | 10.06993 | 18.75136 |
| T | 0.130768 | 0.0176 | 7.430176 | 4.1E-09 | 0.095225 | 0.166312 | 0.095225 | 0.166312 |
我们得到以下的结果:
(1)表中 Coefficient 一栏中载有 Intercept:14.410 65, T:0.130 768. 它们分别是 a, b 的估计,即 $ \hat{a}=14.410\ 65, \hat{b}=0.130\ 768 $. 于是得 A 关于 T 的回归方程为
$$ \hat{A}=14.41065+0.130768T. $$
(2) 表中 p-value 一栏中载有 T:4.1E-09, 这是关于 b 的双边检验 $ H_{0}: b=0, H_{1}: b \neq 0 $ (见第九章 §3(3.19)式)
的p值.由于4.1E-09< $ \alpha $=0.05,故拒绝 $ H_{0} $,认为回归效果是显著的.
(3)表中95%下限一栏中载有T:0.095 225;95%上限一栏中载有T:0.166 312.这表示b的置信水平为0.95的置信区间为
$$ (0.095\;225,0.166\;312). $$
§ 6 bootstrap 方法、宏、VBA
我们在 Excel 环境中求解 bootstrap 问题.
Excel 虽然功能强大,但不可能直接处理各种各样的具体问题。很多情况下需要读者自编称为“宏(Macro)”的程序以解决问题。
“宏”是包括了一连串指令的一个小程序。编写宏要使用VBA(Visual Basic for Application)语言,VBA是Office软件包中的标准语言。
本书没有篇幅介绍宏和 VBA 本身,请读者参阅有关书籍,本书将围绕 bootstrap 方法应用的例子作有关的说明。
现在给出一个宏以解决如下的问题.
问题 设金属元素铂的升华热是具有分布函数 F 的连续型随机变量,F 的中位数 $ \theta $ 是未知参数,测得以下的数据(以 kcal/mol 计):
| 133.2 | 134.1 | 134.3 | 134.4 | 134.5 | 134.7 | 134.8 | 134.8 | 134.8 |
| 134.9 | 134.9 | 135.0 | 135.0 | 135.2 | 135.2 | 135.4 | 135.4 | 135.8 |
| 135.8 | 136.3 | 136.6 | 141.2 | 143.3 | 146.5 | 147.8 | 148.8 |
(数据已经过排序)
(1)以样本中位数 M 作为总体中位数 $ \theta $ 的估计 $ \hat{\theta} $
(i) 求估计量 $ \hat{\theta} $ 的标准误差 $ \sigma_{\hat{\theta}} = \sqrt{D(\hat{\theta})} $ 的 bootstrap 估计.
(ii) 求均方误差 $ MSE = E\left[(M - \theta)^2\right] $ 的 bootstrap 估计.
(iii)求偏差 $ b = E(M - \theta) $ 的 bootstrap 估计.
(2)以样本中位数作为总体中位数 $ \theta $ 的估计,求总体中位数 $ \theta $ 的 bootstrap 置信区间;以样本 20% 截尾均值作为总体 20% 截尾均值 $ \mu_{t} $ 的估计,求截尾均值 $ \mu_{t} $ 的 bootstrap 置信区间.
以下是求解上述问题的宏:
Dim i As Integer, j As Integer, m As Integer, n As Integer
Dim Shu As Single, Temp As Single, Sum As Double, Squ As Double
For i = 1 To 10000
For j = 1 To 26
Shu = Rnd * 26
If Shu <= 1 # Then
Cells(i, j + 5) = 133.2: GoTo 5
ElseIf Shu <= 2 # Then
Cells(i, j + 5) = 134.1: GoTo 5
ElseIf Shu <= 3# Then
Cells(i, j + 5) = 134.3: GoTo 5
ElseIf Shu <= 4# Then
Cells(i, j + 5) = 134.4: GoTo 5
ElseIf Shu <= 5# Then
Cells(i, j + 5) = 134.5: GoTo 5
ElseIf Shu <= 6# Then
Cells(i, j + 5) = 134.7: GoTo 5
ElseIf Shu <= 9# Then
Cells(i, j + 5) = 134.8: GoTo 5
ElseIf Shu <= 11# Then
Cells(i, j + 5) = 134.9: GoTo 5
ElseIf Shu <= 13# Then
Cells(i, j + 5) = 135#: GoTo 5
ElseIf Shu <= 15# Then
Cells(i, j + 5) = 135.2: GoTo 5
ElseIf Shu <= 17# Then
Cells(i, j + 5) = 135.4: GoTo 5
ElseIf Shu <= 19# Then
Cells(i, j + 5) = 135.8: GoTo 5
ElseIf Shu <= 20# Then
Cells(i, j + 5) = 136.3: GoTo 5
ElseIf Shu <= 21# Then
Cells(i, j + 5) = 136.6: GoTo 5
ElseIf Shu <= 22# Then
Cells(i, j + 5) = 141.2: GoTo 5
ElseIf Shu <= 23# Then
Cells(i, j + 5) = 143.3: GoTo 5
ElseIf Shu <= 24# Then
Cells(i, j + 5) = 146.5: GoTo 5
ElseIf Shu <= 25# Then
Cells(i, j + 5) = 147.8: GoTo 5
Else
Cells(i, j + 5) = 148.8
5 End If
Next j
8 For j = 7 To 31
For n = 6 To j - 1
If Cells(i, j) < Cells(i, n) Then
Temp = Cells(i, j)
For m = j - 1 To n Step - 1
Cells(i, m + 1) = Cells(i, m)
Next m
Cells(i, n) = Temp: GoTo 10
End If
Next n
10 Next j
Cells(i, 1) = (Cells(i, 18) + Cells(i, 19)) / 2
Cells(i, 3) = WorksheetFunction.Average(Cells(i, 11), Cells(i, 12), Cells(i, 13), Cells(i, 14), Cells(i, 15), Cells(i, 16), Cells(i, 17), Cells(i, 18), Cells(i, 19), Cells(i, 20), Cells(i, 21), Cells(i, 22), Cells(i, 23).Cells(i, 24), Cells(i, 25), Cells(i, 26))
Cells(i, 2) = Cells(i, 1)
Cells(i, 4) = Cells(i, 3)
Next i
20 Sum = 0 #
For i = 1 To 10000
Sum = Sum + Cells(i, 1)
Next i
Cells(2, 5) = Sum / 10000
Squ = 0 #
For i = 1 To 10000
Squ = Squ + (Cells(i, 1) - Cells(2, 5)) ~ 2
Next i
Cells(4, 5) = Squ / 9999
Cells(6, 5) = Sqr(Cells(4, 5))
Cells(8, 5) = Cells(2, 5) - 135.1
Squ = 0 #
For i = 1 To 10000
Squ = Squ + (Cells(i, 1) - 135.1) ~ 2
Next i
30 Cells(10, 5) = Squ / 10000
End Sub
对这个宏,说明如下:
- 对于一个 Excel 的工作表,设计为:使用工作表的第 1 至第 10 000 行、第 1 至第 31 列(A,B,…,Z,AA,…,AE 列)。用单元格 F1:AE10 000 存放 bootstrap 样本,单元格 A1:A10 000 存放由 bootstrap 样本求出的中位数。单元格 C1:C10 000 存放由 bootstrap 样本求出的截尾均值。单元格 E2、E4、E6、E8、E10 分别存放标准误差、均方误差、偏差的 bootstrap 估计等结果。
- 在宏中出现的变量名和数组要用 Dim 语句来声明 (Dim 是 dimension 的缩写): Integer 表示整数值,从 -32768 至 32767;Long 表示大整数值,从 -2147483648 到 2147483647;Single 表示单精度浮点数,负数从 -3.402823E38 到 -1.401298E-45,正数从 1.401298E-45 到 3.402823E38;Double 表示双精度浮点数,负数从 -1.79769313486232E308 到 -4.94065645841247E-324,正数从 4.94065645841247E-324 到 1.79769313486232E308。
- 使用循环语句,以 i 计工作表中的行(1~10 000),以 j 计工作表中的列(6~31),以 3 标注的程序行至以 5 标注的程序行,藉随机数 Rnd 所处的范围来产生 bootstrap 样本,存放在单元格 Cells(i, j+5) 中.
本例的原样本,n=26,134.8 出现3次,134.9,135.0,135.2,135.4,135.8 各出现2次,其余都只出现1次. 原样本的中位数为135.1,原样本20%截尾均值为135.2875.
在 Excel 环境中,(0,1) 上的均匀分布随机数为 RAND,但在 VBA 环境中,(0,1) 上的均匀分布随机数为 Rnd,在宏中需使用 Rnd.
将 Rnd 乘以 26,得 Shu. 对于某一个确定的 $ j(j=1,2,\cdots,26) $,若 $ 0 < \text{Shu} \leq 1 $,则将 133.2 赋予 Cells(i,j+5);若 $ 1 < \text{Shu} \leq 2 $,则将 134.1 赋予 Cells(i,j+5),…,若 $ 6 < \text{Shu} \leq 9 $,则将 134.8 赋予 Cells(i,j+5),…,若 $ 25 < \text{Shu} < 26 $,则将 148.8 赋予 Cells(i,j+5)。数字之后有 #号,表示该数为浮点数。
接着,对下一个j,用同样的做法对 Cells(i,j+5)赋值.
这样继续做下去,就对第i行的26个单元格 Cells(i,j+5) 赋值,这就得到一个容量为26的 bootstrap 样本.
- 以 8 标注的程序行到以 10 标注的程序行,对自 Cells(i,6) 至 Cells(i,31) 之间的 bootstrap 样本自左至右,从小到大排序。例如,得排序后的样本:
| 133.2 | 134.1 | 134.3 | 134.3 | 134.3 | 134.7 | 134.8 | 134.8 | 134.8 |
| 134.8 | 134.8 | 134.9 | 135.2 | 135.2 | 135.4 | 135.4 | 135.4 | 135.8 |
135.8 136.6 141.2 141.2 146.5 148.8 148.8 148.8 又如,得样本
134.3 134.4 134.7 134.7 134.8 134.9 134.9 135.0 135.0
135.2 135.2 135.4 135.4 135.4 135.8 135.8 136.3 141.2
141.2 141.2 143.3 143.3 143.3 146.5 148.8 148.8
- 排序后,第13个数 Cells(i,18) 和第14个数 Cells(i,19) 的均值就是 bootstrap 样本的中位数,赋予 Cells(i,1),即是放入 A 列.
- 排序后,求20%截尾均值,它就是自 Cells(i,11) 至 Cells(i,26) 的均值。这里,要使用 Average 函数。Average 函数在 VBA 中没有,在 Excel 中才有,在宏中使用(即是在 VBA 环境中使用)要在其前面加上“Worksheet Function.”(凡是在 VBA 中没有而在 Excel 中才有的函数,在 VBA 环境中使用,都要如此处理。)
求得的20%截尾均值,赋予Cells(i,3),即是放入C列.
- 将 A 列数据作一备份,放入 B 列;将 C 列数据作一备份,放入 D 列。这是为了以后 A 列和 C 列数据按列自小到大排序,数据与所处的行的 bootstrap 样本不再对应。
- 以20标注的程序行至以30标注的程序行,对A列数据作计算. 算出10000个bootstrap样本的中位数 $ M_{i}^{*}(i=1,2,\cdots,10000) $的均值 $ \overline{M}^{*}=\frac{1}{10000}\sum_{i=1}^{10000}M_{i}^{*}=135.1438 $. 将它放入单元格E2.计算各中位数与135.1438之差的平方和,除以9999,得0.068133,放入单元格E4.计算0.068133的平方根得0.261023,放入单元格E6,于是E6中的数是
$$ \sqrt{\frac{1}{9\ 999}\sum_{i=1}^{10000}(M_{i}^{*}-\overline{M}^{*})^{2}}=0.261023, $$
这就是我们要求的以样本中位数作为总体中位数的估计的标准误差的 bootstrap 估计.
计算 135.1438 与 135.1 之差得 0.043783,将它放入单元格 E8,计算各中位数与 135.1 之差的平方的均值,得 0.070043,放入单元格 E10。于是 E10 中的数是
$$ \frac{1}{10\ 000}\sum_{i=1}^{10\ 000}(M_{i}^{*}-135.1)^{2}=0.070\ 043. $$
它是以样本中位数作为总体中位数的估计,其均方误差的 bootstrap 估计.
E8 中的数是
$$ \frac{1}{10\ 000}\sum_{i=1}^{10\ 000}M_{i}^{*}\ -135.1=0.043\ 783. $$
它是以样本中位数作为总体中位数的估计,其偏差的 bootstrap 估计.(以上内容请参见第10章§1例1、例2、例3).
- 程序运行结束,运算结果呈现以后,依次将光标指向 A 列列标和工具“ $ \underline{\text{A↓}} $”,单击鼠标左键,根据屏幕提示,在“扩展选定区域”和“以当前选定区域排序”二者之间选定后者,再单击“排序”,即可将列数据自上而下按从小到大排序.然后,依次将光标指向 C 列列标和工具“ $ \underline{\text{A↓}} $”,作同上的操作,即可将 C 列数据自上而下按从小到大排序.
排序以后,可以读出下列数据:
| 行号 | 1 | 250 | 500 | 1 000 | 5 000 | 9 000 | 9 500 | 9 750 | 10 000 |
| 中位数 | 134.6 | 134.8 | 134.85 | 134.9 | 135.1 | 135.4 | 135.6 | 135.8 | 141.8 |
| 截尾均值 | 134.55 | 134.85 | 134.9 | 134.9571 | 135.2571 | 136.3286 | 136.6286 | 136.9214 | 139.7786 |
从而得
| 置信水平 | 0.95 | 0.90 | 0.80 | |
| 置信区间 | 中位数 | (134.8, 135.8) | (134.85, 135.6) | (134.9, 135.4) |
| 截尾均值 | (134.85, 136.9214) | (134.9, 136.6286) | (134.9571, 136.3286) | |
下面是求解本书第十章 §2 例 1(2) 的宏.
Sub Macro7()
Dim i As Integer, j As Integer, k(1 To 10) As Double, kk As Double
For i = 1 To 5000
kk = 0 #
For j = 1 To 10
k(j) = 100.0696 * Sqr(-WorksheetFunction.Ln(Rnd))
Cells(i, 3 + j) = k(j)
kk = kk + k(j) * k(j)
Next j
Cells(i, 1) = Sqr(kk / 10)
Cells(i, 2) = Cells(i, 1)
Next i
End Sub
读了前面的说明以后,这里不需要作更多的说明。
bootstrap 样本容量为 10,共产生 5000 组.
样本1:59.098879.329373.91327111.4051109.506550.55522206.721652.3322145.3297758.67922
$ \eta^{*}=96.5288 $
样本 5 000: 45.357 72 66.256 3 108.903 47.129 5 163.166
90.971 14 83.443 9 52.360 29 43.647 44 81.822 15
$ \eta_{5}^{*}000 = 85.867 12 $
程序运行结束,结果显示以后,点击 A 列列标,再点击工具“ $ \frac{A}{2}\downarrow $”将 A 列数据自上到下,从小到大排序,得:
| 行号 | 1 | 250 | 500 | 2 500 | 5 000 |
|---|---|---|---|---|---|
| $ \eta^{*} $ | 47.274 63 | 73.257 36 | 79.036 52 | 98.675 6 | 177.349 5 |
即到 $ \hat{\eta}_{(250)}^{*}=73.25736\hat{\eta}_{(500)}^{*}=79.03652 $。这就是本书第十章§2例1(2)的中间结果。
本章参考文献
[1] John Walkenbach 等著. Excel 2002 宝典. 牛力等译. 北京:电子工业出版社,2001.
[2] John Walkenbach 著. Excel 2002 公式与函数应用宝典. 路晓村等译. 北京:电子工业出版社,2002.
[3] Gini Counter 等著. Excel 2002 从入门到精通(中文版). 魏江力等译. 北京:电子工业出版社,2002.
[4] M. C. Martin 等著. Excel 2000 从入门到精通(中文版). 惠林等译. 北京:电子工业出版社,2000.
[5] Paul McFedries 著,Office 2000 VBA 编程技术. 韩松等译. 北京:电子工业出版社,2000.