您好、欢迎来到现金彩票网!
当前位置:秒速快3计划 > 枢轴语法 >

概率统计实验(初稿)

发布时间:2019-06-07 06:38 来源:未知 编辑:admin

  . 一般说来,大量重复试验中稀有事件出现的频数 X 均服从或近似服从泊松分布.泊松 分布是由法国数学家泊松于 1837 年作为二项分布的极限分布而导出的,这在当时是一项了 不起的发现.泊松定理告诉我们:若 ( )~,XBnp,其中 n 很大, p 很小,而 np λ= 不太大 (一般说来,要求 λ ≤5)时,则 X 就近似服从泊松分布 P(λ).由此得到一个实际中常用的 近似计算公式: () ()1.! k nkkk np n npCp p ek − −−≈ 本实验就是要通过实际计算、作图和比较等 方法对上述结果在不同参数组合情形下给出直观、动态的展示. 3.2 实验的目的和要求 实验目的是让实验者在计算机上学会二项分布和泊松分布相关概率列、分布函数及分位 数的计算方法,掌握在 Excel 中动态演示二项分布和泊松分布近似关系的方法.具体要求为: 掌握计算二项分布和泊松分布相关概率分布列、分布函数的命令,学会滚动条的制作并能用 滚动条对不同的 n 和 p 进行滚动控制,从而实现对二项分布和泊松分布近似关系的动态演 示并对演示结果进行总结分析. 3.3 实验过程 例 3.1 假如生三胞胎的概率是1 10000 ,试问在 100000 次生育中:(1)恰有 10 次生 三胞胎的概率是多少? (2)至多有 10 次生三胞胎的概率是多少? 解 设 X 表示在 100000 次生育中生三胞胎的次数,则 ( )~ 100000,1 10000XB . (1)所求概率为 () ()()10 9999010 10000010 1 10000 1 1 10000 0.1251PX C== − = . (2)所求概率为 () ()() 10 100000 100000 0 10 1 10000 1 1 10000 0.5830kkk k PX C − = ≤= − =∑ . 易见,当试验次数 n 较大时,二项分布的概率值计算是非常繁琐的,此时建议利用 Excel 命令来进行计算: ¾ 概率值 () (1)kk nk nPX k Cp p −== − 可用命令“ (,, ,0)BINOMDIST k n p= ”计算; ¾ 累积概率 0 () (1) k ii ni n i PX k Cp p − = ≤= −∑ 可用命令“ (,, ,1)BINOMDIST k n p= ”计算. 23 以上面的例 3.1 为例,对于问题(1),在任一空白单元格(这里为 A1)中输入函数命 令“=BINOMDIST(10,100000,1/10000,0)”,回车后返回的计算结果为 0.125116,对 于问题(2),将上述命令的最后一个选项改为 1 即可,计算结果为 0.583040.见图 3.1. 图 3.1 二项分布和泊松分布概率计算 例 3.2 在一个大城市里,调查一种稀有的病症,假如每个人得病的概率都是 1/500,现 随机选出 1000 个人检查,试问这 1000 中有 10 个人得病的概率是多少? 解 设 X 表示在这 1000 人中得病的人数,则 X ~ B(1000,1/500),于是所求概率为 ()()()10 99010 100010 1 500 1 1 500PX C== − .右端计算很麻烦,可利用上例 3.1 介绍的函数 命令,在单元格 A3 内输入“=BINOMDIST(10,1000,1/500,0)”,回车后得到结果 0.00003717 . 注意:这里 n=1000 很大,p=1/500 很小,而 λ=np=1000(1/500)=2,大小适中,所以可 用泊松分布命令来作近似计算,即 ( ) 10 210 2 10!=0.0000381899PX e−=≈ .在 Excel 中, ¾ 概率值 ()! ,0,1,2, k k PX k e kλλ −== = 可用命令“ (,, ,0)POISSON k n p= ”计算; ¾ 累积概率 0 () ! k i i PX k e iλλ − = ≤=∑ 可用命令“ (,, ,1)POISSON k n p= ”计算. 以上面的例 3.2 为例,在 A4 中单元格中输入命令“=POISSON(10, 2, 0)”,返回的计 算结果为 0.00003819,这个值与实际值 0.00003717 很接近,见图 3.1. 例 3.3 制作滚动条演示二项分布概率如何随着 n 和 p 的变化而变化. 仿照实验一的方法,先设置对试验次数 n 的滚动条(n 的值显示在 A3 单元格内):在 Excel 界面中依次点击【开发工具】/【插入】,在出现的【表单控件】对话框中点击【滚动 条(窗体控件)】 ,然后再在工作表空白区任意点击即出现的滚动条 ,并可调整其 位置、大小,再右击此滚动条,然后在出现的对话框中设置控件格式,包括最大值、最小值、 步长、页步长等,关键是在【单元格链接】框中输入要对其值进行滚动的单元格,本例选择 图 3.2 对单元格 A3(试验次数 n)和 A8(A7 中概率 p=A8/100)的滚动条设置 24 $A$3,并设定其最小值、最大值分别为 1、30,步长为 1,页步长为 1,如图 3.2 所示: 确定后完成设置.把该滚动条移动到 A4 处,再点击滚动条的左右小黑三角就会发现 A3 中 的值(即试验次数 n)在 1~30 之间发生大小变化.下面考虑对 p 值的滚动:同上先设置一 个对 A8 的滚动条,滚动范围是 0~100,步长为 1.然后在单元格 A7 内填入“=A8/100”, 那么 A7 的值(即 p 值)就随着 A8 中数值的变化以步长 0.01=1/100 在 0~1 之间相应变化(为 工作表界面简明起见,图中已用后一个滚动条遮住了 A8 单元格),如图 3.3 所示. 图 3.3 演示二项分布概率随着 n 和 p 的变化而变化 其次,在单元格 B2 和 B3 中分别填入 0 和 1,选取单元格区域 B2:B3,将鼠标移至右下 角,出现实心十字时,向下拖放至单元格 B32 放开,就可在单元格区域 B2:B32 内得到随机 变量 X 的取值 k,k=0,1,2,…,30.如图 3.3 所示. 再次,利用函数 BINOMDIST 计算二项分布相应的概率值:先在单元格 C2 中填入的命 令“=IF( B2

  ⎪⎩ (6.2) 这就是根据样本观测值得到的经验分布函数的具体形式. 6.2 实验目的及要求 理解经验分布函数的构成,经验分布函数是样本的函数,随着样本观测值的变化而变化, 通过实验学习经验分布函数图形的绘制方法和动态演示过程.具体要求为 1. 任意产生一组随机样本,对该样本从小到大排序; 2. 利用排序后的样本作经验分布函数图形; 3. 让样本动态变化,观察相应的经验分布函数图形的变化,写出实验体会. 36 6.3 实验过程 为了说明经验分布函数图形的绘制和动态演示过程,我们通过一个实例来进行说明. 例 6.1 在 Excel 中随机产生一个服从均匀分布 (1, 6)U 的,样本容量 =4n 的随机样本.如 图 6.1,在单元格 A2 中输入产生均匀分布 (1, 6)U 随机数命令“=1+5*rand()”,再将其拖放填 充至 A5,就可在单元格区域 A2:A5 中产生了 4 个样本观测值 1234,,,x xxx,每按一次 F9 键, 这些随机数都会发生变化,这为我们进行动态演示带来方便,这里 1 3.82,x = 2 2.81,x = 3 1.32,x = 4 4.44x = . 接着我们把观测值 1234,,,x xxx从小到大排序, 在单元格区域 B2:B5 中分别使用命令“=SMALL($A$2:$A$5, k)” (k=1,2,3,4) 得到顺序样本 观测值: (1) (2) (3) (4)1.32, 2.81, 3.82, 4.44xx x x====. 在此基础上,我们可以利用条件语句和散点图绘制经验分布函数的图形.由(6.2)可知,此时 经验分布函数的表达式为 0, 1.32, 0.25, 1.32 2.81, () 0.5 2.81 3.82, 0.75, 3.82 4.44, 1, 4.44. n x x Fx x x x

  . 我们采用 t -检验. 利用 Excel 的解题过程如下: 1. 原假设 01 2H: 0μμ−≤; 2. 在工作表 A、B 列中分别输入数据; 3. 依次单击主菜单【数据】/【数据分析】/【t-检验:双样本等方差假设】,在弹出的对 线 输入相应选项,单击【确定】按钮后,输出结果如图 9.5 所示. 图 9.4 52 4. 作出判断,在显著水平为α 下,检验的拒绝域为{ }ttα≥ .从图 9.5 中可知,t 的值落 在拒绝域中( 0.05t 2.4234, t (15) 1.7531==),故拒绝原假设,可判断镍合金硬度有所提高. 图 9.5 t-检验: 双样本等方差假设 2.双正态总体方差比的 F-检验 例 9.5 甲、乙两台机床加工某种零件,零件的直径服从正态分布,总体方差反映了加 工精度,为比较两台机床的加工精度有无差别,现从各自加工的零件中分别抽取 7 件产品和 8 件产品,测得其直径为, X(甲机床):16.2 16.4 15.8 15.5 16.7 15.6 15.8; Y(乙机床):15.9 16.0 16.4 16.1 16.5 15.8 15.7 15.0. 试在显著水平 0.05α = 下判断甲、乙两台机床的加工精度有无差别? 显然,这是一个检验双总体方差是否相等的问题,待检验假设为 01 2H: vs σ σ= 11 2H:σ σ≠ . 我们采用 F—检验.利用 Excel 的解题过程如下: 1. 原假设 22 01 2H:σ σ= ; 2. 在工作表 A、B 列中分别输入数据; 3. 依次单击主菜单【数据】/【数据分析】/【F-检验:双样本方差】,在弹出的对线 输入相应选项,单击【确定】按钮后,输出结果如图 9.7 所示. 4. 确定临界值,在 E11 中键入公式“=FINV(0.025, 6, 7)”,得上临界值 5.1186, 在 E12 中键入公式“=FIVN(0.975, 6, 7)”,得下临界值 0.1756. 5. 作出判断,在显著水平为α 下,检验的拒绝域为{ } { }/2 1 /2FF FFαα−≥≤∪ .从图 9.7 中 可知,F 的值落在(0.1756, 5.1186)中,故不拒绝原假设,可认为甲、乙两台机床的加工精 度无显著差别. 53 图 9.6 图 9.7 F-检验 双样本方差分析 9.4 讨论 本实验还可以进一步的推广.从以上实验可以看出,假设检验的关键是提出原假设和备 择假设以及检验统计量.用类似的方法,我们还可进行其它假设检验,如非参数检验、比例 检验和成对样本均值检验等. 54 实验十 皮尔逊拟合检验和列联表独立性检验 10.1 实验原理 设总体 X 的分布函数 )(xF 未知, ),,,( 21 nXXX … 是来自总体 X 的一个样本,提出假 设: )()(: 00 xFxFH = ; )()(: 01 xFxFH ≠ . 此时 )(0 xF 已知,可以含有未知参数,若 )(0 xF 中含有未知参数,一般先利用最大似然估计 给出参数的一个点估计值. 若总体是离散的,则假设为 0H :总体分布律 ii pxXP == )( ( 1, 2, ;i = ip 为已知). 若总体是连续型的,则假设为 0H : X 的密度函数为 )(xf ( )(xf 已知). 当然在 ip 及 )(xf 中均可以含有未知参数,处理方法同 )(0 xF 中的一样. 将随机试验的可能结果的全体分为 k 个互斥事件 kAAA ,,, 21 ,在 0H 成立的条件下 计算 ii pAP =)( ( ki ,,2,1 = ),若将试验进行 n 次,事件 iA 出现的频率为 n ni ( nn k i i =∑ =1 , in 为落入第i 个小区间 iA 的样本值的个数),由大数定律可知,若试验次数很多,在 0H 成 立的前提下, i i pn n − 的值应该比较小,故选用统计量 ∑ = −= k i i ii np npn 1 2 2 )(χ . 若知道统计量 2χ 的分布,就可以用其样本值来判别是否能拒绝 0H 了,皮尔逊(Pearson) 以定理的形式给出了 2χ 的分布. 定理 10.1 若试验次数 n 充分大(一般 50≥n ),则当 0H 成立时,无论总体 X 服从 何种分布,统计量 ∑ = −= k i i ii np npn 1 2 2 )(χ 近似地服从自由度为 1−− rk 的 2χ 分布,其中 r 是 分布中未知参数的个数. 55 统计量 2χ 的分布已经知道了,那么在 0H 成立的条件下,对于给定的显著性水平α , 查表或由 CHIINV 函数得到 )1(2 −− rkαχ 的值,再利用样本值计算出 2χ 的值,如果 )1(22 −−

  rkαχχ ,则拒绝 0H ,否则不能拒绝 0H . 注意:统计软件中常常使用检验统计量的 p 值来判断是否拒绝原假设 统计量检验是根据事先确定的显著性水平α 确定的拒绝域来作出决策,不论检验统计 量的值是大还是小,只要它落入拒绝域就拒绝原假设 0H ,否则就不能拒绝原假设 0H .这 样,无论统计量落在拒绝域的什么位置,你也只能说犯第一类错误的概率是α ,但实际上, α 是犯第一类错误的概率的上限值,统计量落在拒绝域的不同位置,决策时犯第一类错误 的概率是不一样的.如果能把犯第一类错误的概率计算出来,就可以直接用这个概率进行决 策,而不需要管事先给定的显著性水平α .这个犯第一类错误的真实概率就是 p 值.例如, 在卡方右侧单边检验(设自由度为 n)时,若根据样本观测值计算得到的统计量的值为 0T , 则此统计量对应的 p 值就为 2 0(() )pP n Tχ=

  ,这个 p 值的大小,从图形上来看等于 0T 右侧 尾部的面积.目前的统计软件都会给出检验统计量的 p 值. 用 p 值进行决策的规则很简单:如果 p

  −− 或计算其 p 值来进行检验,这里 ( )22(1)(1)pP r cχ χ=−−

  . 当 2rc==时,列联表中只有 4 个格子,称为“四格表”,这时 2χ 简化为 2 11 22 12 21 12 12 ()nnn n n nnnnχ ⋅⋅⋅⋅ −= . 其自由度为 1. 10.2 实验目的及要求 理解 2χ 检验的基本思想,熟练掌握 2χ 检验的分析方法和基本步骤;学会针对实际背景 57 提出原假设和备择假设来检验实际问题,并根据检验结果作出符合统计学原理和实际情况的 判断和结论.掌握列联表独立性检验的原理,能够运用 Excel 对两个分类变量的独立性进行 假设检验. 10.3 实验过程 现将 2χ 检验法的基本步骤总结如下: ○1 提出原假设 )()(: 00 xFxFH = 或 0H :总体服从某种分布(注:若总体分布中 含有 r 个未知参数,则需要先用极大似然估计法估计出这些参数的值); ○2 将总体 X 的取值范围分成 k 个互不重叠的小区间 kAAA ,,, 21 ,一般 165 ≤≤ k ; ○3 计算频数 in ,即所给的 n 个样本观测值 nxxx ,,, 21 中落入第i 个区间 iA 内的个数; ○4在 0H 成立的条件下,计算 X 落入第第i 个区间 iA 的概率 ip ; ○5 计算统计量 ∑ = −= k i i ii np npn 1 2 2 )(χ 的值; ○6 对给定的显著性水平α ,查表或者在 Excel 中由 CHIINV 函数得到 2 (1)krαχ −− 的值, 若 )1(22 −−

  rkαχχ ,则拒绝 0H ,否则不能拒绝 0H . 注意: 目前统计软件中流行的方法是利用检验统计量的 p 值来进行判断:若 p

  ,所 以拒绝原假设 0H ,认为这颗骰子的六个面不均匀. 或者用 p 值来进行判断:由于 0.0011662 0.05p α=

  =,故不能拒绝原假设 0H . 例 10.3 下面是对 123 人进行关于某项政策调查所得结果的一个列联表,在给定的显著 性水平 05.0=α 的条件下,试检验个人观点与收入是否相互独立? 表 10.4 不同收入人群对某项政策的观点 反对 支持 合计 低收入 7 45 52 中收入 15 25 40 高收入 19 12 31 合计 41 82 123 60 解 要检验个人观点与收入是否相互独立,其原假设和备择假设分别为: H0:个人观点和收入相互独立 v.s. H1:个人观点和收入不独立. 下面根据 2χ 检验的基本原理,利用 Excel 创建公式进行独立性检验.其过程如图 10.3 所示.这里显著性水平为 0.05,3 代表列联表的行数,2 代表列联表的列数. 图 10.3 列联表独立性检验 10.4 讨论 1. 自 1965 年 1 月 1 日至 1971 年 2 月 9 日共 2231 天中,全世界记录到里氏震级 4 级和 4 级以上地震计 162 次,统计如下: 表 10.5 相继两次地震记录表 试检验:相继两次地震间隔的天数 X 服从指数分布( =α 0.05). 2. 从一批棉纱中随机抽取300条进行拉力试验, 结果列在下表中, 请检验棉纱的拉力强 度(单位: 公斤) X 是否服从正态分布 (0.01)α = . 61 表 10.6 棉纱拉力数据表 5648.1~34.17 138.2~18.2135334.1~20.16 318.2~04.2123720.1~06.15 1604.2~90.1112506.1~92.04 1990.1~76.110992.0~78.03 2576.1~62.19278.0~64.02 5362.1~48.18164.0~5.01 ii fxifxi 3 为了研究吸烟是否与患肺癌有关,对63 位患肺癌患者及 43 名非肺癌患者(对照组) 调查了其中的吸烟人数,得到 22× 列联表,见表 10.7 所示: 表 10.7 列联表数据 患肺癌 未患肺癌 合计 吸烟 60 32 92 不吸烟 3 11 14 合计 63 43 106 问“是否吸烟”与“是否患肺癌”两者是否统计独立? 4 某研究人员收集了亚洲、欧洲和北美洲人的 A、B、AB、O 血型资料,结果如表 10.8 所示,问不同地区的人群血型分类构成是否一样。 表 10.8 三个不同地区血型样本的频数分布 地区 A B AB O 合计 亚洲 321 369 95 295 1080 欧洲 516 86 44 388 1034 北美洲 408 106 37 444 995 合计 1245 561 176 1127 3109 62 实验十一 单因素方差分析 11.1 实验原理 方差分析是利用实验数据,分析多个因素(如品种、施肥量等)对一事物某指标(如平 均亩产量)的影响是否显著的一种统计分析方法.单因素方差分析是方差分析中最简单的一 种,即影响该指标的因素只有一个. 具体而言:设一个因素 A 有 r 个水平,在这 r 个水平(对 应 r 个总体)下分别进行试验,得到来自这 r 个总体的试验数据,现在要分析这 r 个总体的 该指标是否有显著差异. 单因素方差分析的三个基本假设: (1)每个总体都应服从正态分布.也就是说,对于因素的每一个水平,其观察数据是来 自服从正态分布总体的简单随机样本. (2)各个总体的方差 2σ 虽然未知,但要假设它们相同.也就是说,各组观察数据是从 具有相同方差的总体中抽取的. (3)观察数据是独立的. 单因素方差分析的数学模型 设因素 A 有 r 个水平 1,,rAA… ,在水平 iA 下进行 (2iinn≥ )次独立试验,试验结果如 表 11.1. 表 11.1 样本观察值 水平 样本观察值 1A 2A rA 11x 12x … 11nx 21x 22x … 22nx 1rx 2rx … rrnx 其中 ijx 表示第i 水平 iA 下进行第 j 次试验结果. 设 2~(,)ij ixNμ σ ,则方差分析的模型为 22 , , 1,2, , ; 1,2, , , ~(0,), .iij ij i ij i ij x irjn N μσ ε με εσ =+ = =⎧ ⎨ ⎩ 未知,且各 相互独立 (11.1) 对于模型(11.1),待检验的假设为 01 2 112::,,,rrHvsHμ μμ μμμ===不全相等. 63 记 2 11 ), inr ij ij SS x x == −∑∑=( 2 . 11 ), inr eiji ij SS x x == −∑∑=( 2 . 11 ). inr Ai ij SS x x == −∑∑=( 则 .eASS SS SS=+ 理论上可以证明:在 0H 为真的条件下 / ~( , )./ AA Ae ee SS dfF F df dfSS df= 所以,拒绝域为 { }(,)AeWFFdfdfα≥= .其中 ,Aedf df 分别表示 ASS 和 eSS 的自由度.对给 定α ,可作如下判断: (1)如果 (,)AeF F df dfα

  ,则拒绝原假设 0H ,认为因素 A 显著; (2)如果 (,)AeF F df dfα≤ ,则不拒绝原假设 0H ,认为因素 A 不显著. 11.2 实验目的及要求 掌握 Excel 自带的数据分析工具中“方差分析:单因素方差分析”工具的具体操作步骤; 对方差分析的结果能够作出正确合理的解释,从而判断该因素的各水平是否有显著差异,并 作出相应的决策. 11.3 实验过程 例 11.1 为了研究咖啡因对人体功能的影响,特选 30 名体质大致相同的健康的男大学 生进行手指叩击训练,此外咖啡因选三个水平: 12 30( ), 100( ), 200( ).AmgAmgAmg== = 每个水平下冲泡 10 杯水,外观无差别,并加以编号,然后让 30 位大学生每人从中任选一杯 服下,2 小时后,请每人做手指叩击,统计员记录其每分钟叩击次数,试验结果如下表: 表 11.1 30 个大学生叩击次数记录 请对上述数据进行方差分析,考查咖啡因对人体功能是否有影响. (1)建立数据文件如图 11.1 图 11.1 64 (2)从主菜单中选择【数据】/【分析】/【数据分析】,在【数据分析】对话框中选 择【方差分析:单因素方差分析】选项,如图 11.2 图 11.2 (3)单击【确定】按钮,弹出【方差分析:单因素方差分析】对线)单击确定,单因素方差分析结果如图 11.4 所示. 图 11.4 (4)结果分析:因为 F 的值 6.18 大于 F 临界值 3.35,所以可认为因素 A 显著,即可 认为咖啡因对人体功能是有显著影响的.这也可以由 0.00616 0.05p =

  , (,)()0PX iY j P===∅=,所以 (,)X Y 的联合分布列为 77 表 14.1 三项分布联合分布列的计算 Y X 0 1 2 3 0 0.001 0.009 0.027 0.027 1 0.018 0.108 0.162 0 2 0.108 0.324 0 0 3 0.216 0 0 0 该联合分布列的计算过程比较繁琐,但在 Excel 中用阶乘函数命令 FACT(k)和乘幂函数 命令 POWER(n, k)并利用 Excel 的拖放填充功能可以准确快速的计算出结果. 如图 14.1,先 在单元格区域 B2:B4 分别输入抽到一等品、二等品和三等品的概率,然后在 B5 中输入抽取 的产品件数 n,接着在单元格区域 C2:C5 和 D1:G1 中输入随机变量 X 和Y 可能的取值,然 后在单元格 D2 中输入公式: “=IF(($C2+D$1)

  3,0,(FACT($B$5)/FACT($C2)/FACT(D$1)/FACT($B$5-$C2-D$1)) *POWER($B$2,$C2)*POWER($B$3,D$1)*POWER($B$4,$B$5-$C2-D$1))”, 用此公式可计算 003003! 6 3 1(0,0) ()()()0!0!(3 0 0)! 10 10 10PX Y − −=== −− . 确定后得所求概率值 0.001. 其中 IF(logical_test,value_if_true,value_if_false)用于创建条件 函数,表示“如果 3ij+

  ,则, (,)PX iY j= = =0,否则按公式(3.1.6)计算概率”, POWER(number, power)用来计算数字的乘幂,POWER($B$2,$C2)即 06()10 . 然后鼠标点击单 元格 D2 右下角,等出现小黑十字后按住不放,拖动鼠标至 G5 放开,就可以自动计算出所 对应的全部概率值 (,)PX iY j==, 得到(,)X Y 的联合分布列. 图 14.1 三项分布计算结果 (2) 多维超几何分布 多维超几何分布是超几何分布的自然推广,它可这样来描述: 设一批产品共有N件,其中有Ni件i等品, 1, 2, ,ir= , 显然有 12 rNN N N+ ++ = . 78 现从这批产品中不放回地任意取出 n 件,记 iX 为这 n 件产品中 i 等品的件数, 1, 2, ,ir= , 则 12( , , , )rXX X 取值 12( , , , )rnn n 的概率为 12 12 1122(, ,,) , r r nn n NN N rr n N CC CPX n X n X n C== == (14.3) 其中 12 rnn n n+++= . 称此联合分布为 r 维超几何分布. 当 r = 3 时称为三维超几何分 布,它是最简单的多维超几何分布. 例 14.2 将例 14.1 中的有放回抽样改为不放回抽样,其余条件不变,求二维随机变量 (,)X Y 的联合分布列. 解 本例和例 14.1 类似,计算结果如图 14.2,先在单元格 B2 中的输入计算公式: “=IF($A2+B$1

  3,0,COMBIN(60,$A2)*COMBIN(30,B$1)*COMBIN(10,3-$A2-B$1)/COMB IN(100,3))”(请读者思考该公式的含义),拖放填充至 E5 后放开可得到联合分布列. 图 14.2 三项超几何分布计算结果 79 实验十五 定积分的近似计算 15.1 实验原理 定积分近似计算是实际问题中经常遇到的数学问题之一,在原子能研究和参数估计等问 题中会出现许多复杂的单重与多重定积分,通常,人们都选用矩形公式、辛普生公式等来完 成积分的近似计算.在许多情况下,使用这些近似计算公式虽然能够得到相当满意的结果, 但计算量随着积分重数的增加而显著的增加,以致于用计算机都无法完成. 定积分的计算是 Monte Carlo 方法引入计算数学的开端,在实际问题中,我们经常遇到 计算多重积分的复杂问题,用 Monte Carlo 方法一般都能够很有效地予以解决,尽管 Monte Carlo 方法计算结果的精度不很高,但它能很快地提供出一个低精度的模拟结果也是很有价 值的.而且,在重积分中,由于 Monte Carlo 方法的计算误差与积分重数无关,因此它比常 规方法在近似计算中更具优势. 1 随机投点法 设 ()f x 在[0,1]上连续,且有 0()1fx≤ ≤ ,求 1 0 ()s fxdx= ∫ . 设 (,)X Y 服从[0,1;0,1]上的均匀分布,则可知 X Y和 服从[0,1]上的均匀分布,且 X Y和 独 立.因为 ()11 00 0 (()) () fx pPYfX dydxfxdxs=≤ = = =∫∫ ∫ , 所以,定积分的值 s 就是事件 (())AYfX=≤ 的概率 p . 这样,根据大数定律,我们可以 用频率来近似概率. 具体做法如下:先产生 2 [0,1]n 个 上均匀的随机数: 12 12,,,;,,,nnx xxyyy.然后 ),iinxy对 对数据( , 记录满足 ()iiyfx≤ 的次数 k ,最后得 s 的估计值,即 ks n≈ . 注:为了计算一般的定积分,我们可以这样来解决.设 () b a s fxdx= ∫ . 其中 () [,]f xab在 有界,不妨设 ()LfxM≤≤.令 ()x abay= +− ,则有 11 00 ()() ( ( ))( ) ( )[( ( )) ] b a bas f x dx f a b a y b a dy M L f a b a y L L dyML −==+−−=−+−−+−∫∫ ∫ 11 00 ()( )[(())] ()baM LfabayL dyLbadyML −=− +−− +−−∫∫ 1 0 (( ))( )() ()fa b ay LM Lb a dyLbaML +− −=− − +−−∫ . 1 1 0 (( ))fa b ay Ls dyML +− −= −∫令 ,我们只需求 1s 即可. 80 2.平均值法 求 () b a s fxdx= ∫ ,其中 ()f x 在[,]ab上可积. 设 X 服从[,]ab上的均匀分布,则 ()YfX= 的数学期望为 1(( )) () b a EfX fxdxsba= =− ∫ . 所以估计 s 的值就是估计 ()fX的数学期望的值.由辛钦大数定律,可以用 ()fX的观察 值的平均值去估计 ()fX的数学期望.具体做法如下:先产生 [,]nab个 上均匀的随机数: 12,,,nx xx .然后对每个 ix 计算 ()ifx ,最后 s 的估计值,即 1 1 () n i i s fxn = ≈ ∑ . 15.2 实验目的及要求 加深对大数定律的理解,增强动手能力,掌握利用 Excel 进行随机模拟的方法,提高解 决实际问题的能力. 15.3 实验过程 例 15.1(随机投点法) 估计定积分 1 1 0 1 1 xeJ dxe −= −∫ .(当01x≤ ≤ 时, 10() 11 xefx e −≤ =≤− ). 随机投点法的具体步骤为: (1) 独立地产生 2n 个服从(0,1)上均匀分布的随机数, 12 12,,,;,,,nnx xxyyy; (2) 统计 ()iiyfx≤ 的次数 k ; (3) 用 k n 来估计 1J . 利用Excel进行估计的过程如下: (1)从主菜单中选择【数据】/【数据分析】,在【数据分析】对话框中滚动列表框,选 择“随机数发生器”选项,如图 15.1. (2) 由以上操作可得 1000 对随机数,作变换 ()iizfx= ,然后进行比较 iiyz≤ ,如图 15.2, 最后对 D l 列进行计数(COUNTIF(D2:D1002,1)),得 424k = ; (3) 1 0.424kJ n≈= . 图 15.1 图 15.2 81 例 15.2(随机投点法) 估计定积分 1 2 1 xJ edx− = ∫ . 因为 1 2 1 xJ edx− = ∫ 1 12 1 11 1 0 2( ) 2 yeeee dyeee −+ − −− − −=− +−∫ ,所以我们只需估计 1 12 1 2 1 0 yees dyee −+ − − −= −∫ , 重复上面的过程可得 2 0.363s = ,所以, 11 222( ) 2J ee s e− −=− + 2.442= . 例 15.3(平均值法) 估计定积分 1 1 0 1 .1 xeJ dxe −= −∫ (当01x≤ ≤ 时, 10() 11 xefx e −≤ =≤− ). 随机投点法的具体步骤为: (1) 独立地产生 n 个服从 (0,1) 区间上的均匀分布的随机数 12,,,nx xx ; (2) 计算 ()ifx ; (3) 用 1 1 () n i i fxn = ∑ 来估计 1J . 利用Excel进行估计的过程如下: 从主菜单中选择【数据】/【数据分析】,在【数据分析】对话框中滚动列表框,选择“随 机数发生器”选项,产生 1000 个 (0,1) 区间上的均匀分布的随机数,计算 ()ifx ,然后用 AVERAGE 命令求平均值 1 1 () n i i f xn = ∑ , 1 0.4150J ≈ ,如图 15.3. 图 15.3 图 15.4 例 15.4(平均值法) 估计定积分 1 2 1 xJ edx− = ∫ . 首先产生 (1,1)− 的均匀随机数,其它类似,如图 15.4. 2 2 1.19 2.38J ×==. 从上面所做的 1000 次模拟定积分试验的结果中可以看出,随机投点法的估计定积分 1J 和 2J 分别为 0.424,2.442,平均值法估计的定积分 1J 和 2J 分别为 0.415,2.380,精确值分别 为 0.418,2.350,平均值法估计的定积分值比随机投点法的估计定积分值精度高. 15.4 讨论 上述两种算法也适用于多重积分,不会发生原则性的困难,这正是 Monte-Carlo 方法的 优点.我们只简单地叙述一下平均值法,若 ()g x 是 m 元函数,它在 mR 中有限闭区域 D 上可 积,估计积分 12 12(, , )mmD J g x x x dx dx dx= ∫∫ .设 1(,, )nX X 服从 D 均匀分布, D 的 m 维 82 体积记为 D ,则 12 12 12 1(( , , ) ( , , )mmmD EgX X X gx x x dxdx dxD= ∫∫ . 根据此式,我们产生 n 个 D 上的 m 维均匀分布随机数 12(, , ),1,,ii imx xxi n=,计算函数 值 12(, , )ii img xx x ,用平均值 12 1 ˆ (, , ) n ii im i DJ gx x xn = = ∑ . 近似积分值 J . 另外,从以上我们可以看出,估计一个参数可能有多种方法,它们的结果可能都不同, 哪一种方法好,这是值得研究的问题. 83 实验十六 高尔顿钉板试验及其在 Excel 中的实现 16.1 实验原理 高尔顿钉板试验是英国统计学家高尔顿设计的用来验证中心极限定理的著名试验.这个 试验是在一块垂直放置的木板上钉上一排排互相平行、水平间隔相等的 n 排钉子,第一排 1 颗,第二排 2 颗,……,第 n 排 n 颗,下一排中各个钉子正好对准上面一排两相邻的钉子的 正中央.从入口处(第一排钉子的正上方)放入一个直径略小于两颗钉子间隔的小球,小球 碰到第一排的这颗钉子后,等可能地从其左边或者右边落下,接着小球会碰到第二排的钉子, 再等可能地从其左边或者右边落下.在不断下落的过程中,小球会碰到每一排钉子中的某颗 钉子,并且从其左边落下与从其右边落下的机会相等,碰到下一排钉子时也是如此,小球这 样不断下落,最后落入放置在木板最下方的一排格子内.因此,任意从第一排钉子的正上方 放下一个小球,那么此小球最后落入最下方的哪一个格子内是事先无法确定的.独立地重复 本试验 N 次(N 充分大),则试验结果表明, 最后堆积在下方格子中小球的轮廓线形状总是 近似于正态分布密度曲线,这就是独立同分布中心极限定理的一个直观形象的诠释,即独立 同分布(这里为 p = 0.5 的两点分布)随机变量和的极限分布近似于正态分布. 作为一种简洁的教学辅助试验,高尔顿钉板试验常常出现在大学及高中的教材和相关文 献中,用于中心极限定理的直观而生动的演示.下面主要讨论如何在 Excel 中对高尔顿钉板 试验进行模拟并在独立情形下将高尔顿钉板试验推广到二维情形. 16.2 实验目的及要求 掌握高尔顿钉板试验的设计原理,加深对独立同分布中心极限定理的理解,学会利用 Excel 中的随机数发生器、随机数函数命令、条件语句、模拟运算表和柱形图等来模拟演示 高尔顿钉板试验. 16.3 实验过程 一、一维高尔顿钉板试验的设计原理及模拟 1 设计原理 在本实验在 Excel 中用伯努利随机数(即 0-1 随机数)和条件语句来模拟高尔顿钉板试 验.用单元格表示小球可能下落经过的位置.如果单元格内数字为 1,则表示小球落经此单 元格;如果单元格数字为 0,则表示小球没有落经此单元格.因此小球下落的轨迹可以用类 似于图 16.1 中的一串数字“1”来表示.每放下一个小球,其下落的轨迹都在随机地发生变化. 图 16.1 小球下落的轨迹 84 为叙述明确起见,我们用(i, j)来表示 Excel 中,位于第 i 行,第 j 列的单元格.再令 1, ( , ), 0, ( , ).ij ijX ij ⎧= ⎨ ⎩ 小球落经单元格 小球未落经单元格 其中 i, j=1,2,…,n. 首先对小球下落过程中,穿过前三行(或称前三排)钉子时相应单元格的取值进行分析: 第一行:假设小球首先从单元格(1,k)处下落,则有 X1k=1. 事实上,其后的每一个小球 都从此单元格放下,因此 X1k 总是取值 1. 第二行:小球可能落经第 2 行两个单元格之一: (2,k-1) 或 (2,k+1). 小球落经哪个单元格,取决于从第 1 行下落的小球是从(1,k)的左方落下还是从其向右方落 下.由于向左或向右是随机的,等可能的,所以我们可在单元格(2,k-1)处利用 Excel 随机数函 数以概率 0.5 产生一个 0-1 随机数 X2,k-1,即 X2,k-1 取 0 或 1 的概率都是 0.5,注 意 X2,k-1 和 X2,k+1 分别位于“钉子 X1,k”的左下侧和右下侧,因此有 X2,k+1=1-x2,k-1. 在单元格(2,k-1)中产生 0-1 随机数的公式为 X2,k-1“=RANDBETWEEN(0,1)”. Excel 中函数 RANDBETWEEN(m, n)等可能地返回 m 与 n 这两个整数之间的任意一个整数. 第三行:小球可能落经三个单元格之一: (3,k-2), (3,k) 或 (3,k+2). 小球落经哪个单元格,取决于小球在第二行时落经哪个单元格以及是从其左方落下还是从其 右方落下.下面对这三个单元格可能的取值依次进行分析: (1) 小球要落经(3,k-2),前提只可能为小球落经(2,k-1)处并从其左方落下; (2) 小球要落经(3,k),有两种情况:①小球落经(2,k-1) 处并从其右方落下;②小球落经 (2,k+1)处并从其左方落下; (3)小球要落经(3,k+2),前提只能为小球落经(2,k+1)处并从其右方落下. 综合以上分析,可知这三个单元格中产生随机数的命令分别应该为: X3,k-2“=IF(X2,k-1=1,RANDBETWEEN(0,1),0)”; X3,k“=IF(X2,k-1=1,IF(X3,k-2=1,0,1),IF(X2,k+1=1,RANDBETWEEN(0,1),0))”; X3,k+2“=IF(X2,k+1=1,IF(X3,k=1,0,1),0)”. 按照上面的规律依次进行下去,我们可以写出小球在落经任何一行时,该行单元格中的 公式的具体形式.当小球下落落经第 n 行时,我们可以用下面的 Excel 函数命令和条件语句 确定第 n 行中的 n 个单元格中的公式如下: (1)当 j=k-(n-1)时: Xn,j“=IF(Xn-1,j-1=1,RANDBETWEEN(0,1),0)”; (2)当 k-(n-3)≤j≤k+(n-3)时: Xn,j“=IF(Xn-1,j-1=1, IF(Xn,j-2=1,0,1), IF(Xn-1,j+1=1,RANDBETWEEN(0,1),0))” (3)当 j=k+(n-1)时: Xn,j“=IF(Xn-1,j-1=1, IF(Xn,j-2=1,0,1),0)”. 其中 j表示列数,k≥n. 85 通过上面的算法,我们可以确定一个小球从第一行钉子上方落下(相当于从上面的单元 格(1,k)处落下),连续地逐次与各行中的某颗钉子发生碰撞,直到与第 n-1 行钉子中的某颗 钉子碰撞后再下落,再碰到第 n 行钉子中的某颗,该颗钉子所在的单元格可看成是小球下落 的最终位置,即上面所说的小球最终落入的那一排格子.由于在一次实验中,小球最后只能 落入一个单元格内,所以,最后第 n 行钉子所在的 n 个单元格中,有 n-1 个数字是 0,只有 一个数字是 1,而数字 1 所在的单元格就是小球最终落入得格子. 为了记录这次试验的结果,并用模拟运算表独立地重复模拟这个试验,我们需在第 n 行 单元格(如图 16.2 中的 B7:N7)之下再另取一行单元格(如图 16.2 中的 B10:N10)来重新记 录一次上述试验结果,即有公式:B10“=B7”, D10“= D7”,…, N10“=N7”. 为了独立地重复模拟这个试验,我们需要借助 Excel 中的模拟运算表.模拟运算表是一 个单元格区域,它可显示一个或多个公式中替换不同值时的结果.有两种类型的模拟运算表: 单输入模拟运算表和双输入模拟运算表.单输入模拟运算表中,用户可以对一个变量键入不 同的值从而查看它对一个或多个公式的影响.双输入模拟运算表中,用户对两个变量输入不 同值,而查看它对同一个公式的影响. 本试验是从第一排钉子上方放下一个小球开始的,也相当于从单元格(1,k)中生成一个数 字 1 开始的,其后小球的下落过程会自动完成.所以只要在单元格(1,k)中连续地生成数字 1, 就可不断地重复这个试验.设模拟 100 次小球下落过程.用图 16.2 来说明,那么(1,k)相当 于 H1,在 H1 中需要不断生成的那 100 个 1,这 100 个 1 被放置在单元格 A11:A110 中以便 进行模拟运算.为节省篇幅,图 16.2 对第 13 行至第 107 行单元格作了隐藏.在具体模拟运 算时,将单元格 H1 要引用列的单元格区域取为 A11:A110,就可以按照第一次得到 B10:N10 中 0-1 随机数的方法独立重复地模拟小球下落过程 100 次,这 100 次模拟试验的结果分别由 单元格区域 B11:N110 中的 100 行 0-1 随机数所给出. 最后,我们可用求和函数 SUM 把这 100 行 0-1 随机数按列相加,统计出 100 次模拟试 验后最后一行每一个格子内所容纳的小球的个数,即图 16.3 中 Q2,R2,…,W2 中的数值, 并绘出柱形图. 2 模拟试验结果 下面是在 Excel 中模拟 100 次高尔顿钉板试验的结果,这里钉子设计了 7 行. 图 16.2 100 个小球,7 行钉子的高尔顿钉板试验结果 86 对图 16.2 表格中的部分公式的说明: G2“= RANDBETWEEN(0,1)”; I2“=1-G2”; …… G4“=IF(F3=1,IF(E4=1,0,1),IF(H3=1,RANDBETWEEN(0,1),0))”; …… B10“=B7”, D10“=D7”,…, N10“=N7”; 图 16.3 频数表和柱形图 对图 16.3 中表格部分公式的说明: Q2“=SUM(B11:B110)”;R2“=SUM(D11:D110)”;……;W2“=SUM(W11:W110)”. 其中 B,D,…,W 列中的数据见图 16.2. 3 试验结果分析 从上面所做的 100 次模拟高尔顿钉板试验的结果中,可以看出,所得到的频率柱形图的 轮廓线大致为中间高、两边低且呈轴对称的倒钟形曲线,接近于正态分布的密度函数曲线.注 意这里只设计了 7 行钉子,试验次数也只有 100 次,就大致呈现出了高尔顿钉板试验的效果, 即独立同分布随机变量序列和的极限分布近似于正态分布,这就通过模拟试验验证了独立同 分布情形下的中心极限定理.当钉子的行数和试验次数越来越大,这样的近似程度就会越来 越高,高尔顿钉板试验的效果将会更加明显. 上面的模拟试验还有一个突出的优点:可以进行动态演示.由于上述试验基于随机数的 产生,在 Excel 中,只要按一次 F9 键,试验就会重复进行一次,同时图 16.1 及图 16.2 中的那 一列“1”还会随之变化,显示出小球下落的轨迹,另外,频率柱形图也会同时随之发生相 应变化.这一点对中心极限定理的教学和动态展示特别有用. 二、二维高尔顿钉板试验的设计原理及模拟 1 设计原理 上面讨论的高尔顿钉板试验可以理解为一维高尔顿钉板试验,它实际上是一个小球在一 个垂直平面上下落的情形.我们现在把这个实验推广到二维情形:让一个小球从空间某处自 由下落,穿过其下方一层一层的网格线,每一层网格线由距离相等并且垂直相交两组平行线 构成.下一层网格线的交叉点总是对准其上一层网格线的格子中点,小球每次碰到上一层格 子交叉点后,等可能地从以这个交叉点为公共点的四个格子之一落下,再碰到下一层网格中 对准此格子中点的一个交叉点,再随机地从以这个交叉点为公共点的四个格子之一落下,如 87 此等等.最后小球落入置放在最后一层网格线下方的 n×n 个小正方形格子中的某一个格 内.独立重复地做 N 次这样的实验,可以想象小球将在最下方的 n×n 个小正方形格子中堆 积成一个“小山头”,中间高,周围逐渐变低,呈中心轴对称,形状接近于一个倒扣的古钟. 我们可以类似于上述一维高尔顿钉板试验,在 Excel 中对这样的二维高尔顿钉板试验进 行模拟.实验的原理是类似的,但相应的 Excel 命令较为复杂.为简单起见,我们只考虑独 立情形下的二维高尔顿钉板实验. 设二维随机变量(X, Y )的两个分量 X, Y 相互独立,取值为( xi, yj ), (i, j=1,2,…,n),由 X, Y 的独立性可得(X, Y )的联合分布列为 pij=P(X=xi,Y=yj)=P(X=xi)P(Y=yj), i, j=1,2,…,n. 把上面的一维高尔顿钉板试验独立重复地做两个,用 xi 表示前一个高尔顿钉板试验中最 后落入第 i 个格子中的小球个数,i=1,2,…,n;同理,用 yj 来表示后一个高尔顿钉板试验中最 后落入第 j 个格子中的小球个数,j=1,2,…,n. 见图 16.4 中单元格区域 B112:N112 和 Q112:AC112 中的数值.由大数定律知 P(X=xi)≈ni /N(即频数/试验次数);同理,P(Y=yj)≈ nj /N. 故有 pij≈( ni /N)( nj /N) ≈( ni nj /N 2). 图 16.6 中给出了 ni nj /N 2 的模拟值,并据此做出了 三维柱形图,其中 N=100. 2 模拟试验结果 下面模拟的是 100 个小球落在一个 77× 的平面单元格区域内的二维高尔顿钉板试验. 首先按上述方法,独立地做两个一维高尔顿模拟试验,并计算分别落入 7 个格子中小球 的频数,见图 16.4 中单元格区域 B112:N112 和 Q112:AC112 中的数值. 图 16.4 两个一维的高尔顿钉板实验 其次,利用两个一维高尔顿钉板实验数据绘出频数柱形图,见图 16.5. 图 16.5 柱状图 88 最后,利用独立性,给出二维联合分布列,并绘出三维立体柱形图,见图 16.6. 图 16.6 二维高尔顿钉板试验的联合概率图 对图 16.6 表格中部分公式的说明(参见图 16.4): AX2“=B112”; AY2“=D112”;… ;BD2“=N112”;AW3“=Q112”;AW4“=S112”;… ; AW9=AC112;… ;AY4 ≈AW4*AY2/10000=0.0048…等等. 3 试验结果分析 在独立情形下我们将一维高尔顿试验推广到了二维情形.通过对二维高尔顿钉板试验数 据模拟结果以及柱形图的观察,可发现该三维柱形图的轮廓曲面接近于二维正态分布的联合 密度的图形.理论上也可证明,在独立同分布条件下,二维情形下的中心极限定理仍然成 立.另外,二维情形也同样可以进行动态演示,在 Excel 中,逐次按 F9 键,即可不断地实现 二维高尔顿钉板试验的动态演示,图 16.6 中的三维柱形图也会随之不断变化,但其轮廓面总 是近似于中间高,周围逐渐变低,呈中心轴对称的倒扣的古钟形状. 16.4 讨论 我们知道,Excel 几乎是每一台计算机的必装软件,接触过计算机的人对它都有或多或 少的了解,目前不管是大学概率统计教学还是中学的新课改教学内容,都会涉及到高尔顿钉 板试验的模拟.对高尔顿钉板试验的模拟并不是那么容易的,要么需要使用复杂的软件,要 么要进行复杂的编程.而使用大众化的统计软件 Excel 来完成对高尔顿钉板试验的模拟就显 得很有理论意义和实用价值,本实验所介绍的方法并不复杂,人人都可以做,还可以进行推 广.对二维高尔顿钉板试验的模拟从独立情形推广到非独立情形,可以想象,这时情况会比 较复杂,感兴趣的读者可进行更深一步的讨论. 89 实验十七 样本均值的抽样分布模拟 17.1 实验原理 样本均值 1 1 n i i X Xn = = ∑ 和修正样本方差 () 2 2 1 1 1 n i i SXXn = =−− ∑ 为是最常用的统计 量.关于 X 的抽样分布,有如下重要结论: 定理 17.1 设 12,,,nX XX 是来自总体 X 的简单随机样本, X 为样本均值. (1) 若总体 X 的分布为 2(, )N μ σ ,则 X 的精确分布为 2(, )Nnμσ ; (2) 若总体 X 分布未知或者不是正态分布,但 ()EX μ= , 2()Var X σ= ,则当 n 充分 大时 X 的渐进分布也近似为 2(, )Nnμσ . 该定理利用样本的独立性和中心极限定理容易证明.定理结论告诉我们,不管样本是否 来自正态总体,只要样本容量 n 充分大,则 X 的分布就为(或近似为) 2(, )Nnμσ .说明 样本均值的期望等于总体均值,而样本均值的方差 ()Var X 仅为总体方差 2σ 的1 n .或者说, 总体标准差σ 是样本均值标准差 ()Var X 的 n 倍,即 ()nVar Xσ = ,故均值运算能有效降 低数据的波动误差.再注意 2() ()XE SVarX= ,于是近似有 XnSσ ≈ . 17.2 实验目的及要求 实验目的是通过实验模拟验证上述定理的结论.要求从正态总体中抽取多份简单随机样 本,计算每份样本的均值,考察样本均值的分布特征并用柱形图和折线图进行比较,对非正 态总体情形可作类似分析,留作讨论题. 17.3 实验过程 例 17.1 产生 1000 份大小为 9 的 2(6, )N σ 随机样本,计算每份样本的均值 ( 1,2, ,1000)iXi= .再计算这 1000 个样本均值的均值 X 和样本均值的标准差 XS .在总体 方差 =1σ 和 =4σ 两种情形下比较总体分布和样本均值 X 分布的差异,并用柱形图和折线图 进行比较. 解 先在单元格 B3 内输入动态随机数命令“=NORMINV(RAND( ), $C$1, $E$1)”,然 后将此公式拖放填充至 J1002,就可以在单元格区域 B3:J1002 中生成 1000 份容量均为 9 的 (6,16)N 随机样本,每一行(如 B2:J2)就是一份随机样本,如图 17.1. 接着在单元格 K3 内 90 输入函数命令“=AVERAGE(B2:J2)”计算出该份样本的均值,再将其拖放填充至 K1002, 就计算出了这 1000 份样本的均值,它们所在的单元格区域为 K3:K1002.在单元格 L3 内输 入函数命令“=AVERAGE($K$3:$K$1002)”计算出这 1000 个样本均值的均值 X =5.999,这 个值非常接近于总体均值 6;在单元格 M3 内输入函数命令“=STDEV($K$3:$K$1002)”计 算出这 1000 个样本均值的标准差 XS =1.313,这个值近似等于总体标准差σ =4 的 1=19=13n ,即有 9 =3.938 4=XS σ≈ .这就验证了定理 17.1 的结论.为了从图形上获 得直观认识,我们在单元格区域 L11:L35 中以间隔为 0.5 输入区间分点:0,0.5,1,…,12. 然后选取单元格区域 M11:M35,再在编辑栏中输入公式“=FREQENCY($K$3:$K$1002, L11:L35)”,按 Ctrl+Shift+Enter 组合键(这是数组命令,不能仅按 Enter 键),得到单元格 区域 M11:M35 中的频数计算结果,如图 17.1. 再以这组频数(即 M11:M35 中频数)为数据 图 17.1 1000 份容量为 9 的正态总体 N(6,16)样本均值的均值、标准差、频数和频率 系列,以 L11:L35 为水平(分类)轴标签作柱形图,其轮廓线显然很近似于正态分布密度曲 线 左图.若将每个小区间内频数均除以 1000 得到每个小区间内的频率(即 N11:N35 中频率).再以 L11:L35 和 N11:N35 中的数据作频率散点图,并与 N(6,16)密度曲 线进行比较,可以看出,该频率散点图比 N(6,16)密度曲线的分布要更集中、陡峭.两者关 系是:样本均值的方差 2() XVar X S≈ 仅为总体方差 22=4 =16σ 的19,或者说,样本均值的标 准差 () XVar X S≈ 仅为总体标准差 =4σ 的13.见图 17.1 中单元格 N3 和 E1 中的值,它们 91 近似有关系式 9=31.313=3.9384=XS σ×≈.参见图 17.2 右图. 图 17.2 正态总体 N(6,16)样本均值柱形图及频率折线)密度曲线的比较 类似考虑总体为 (6,1)N 的情形.将总体标准差σ (即 E1 中的值)改为 1(注:图 17.3 和图 17.1 中均设置了滚动条对 E1 进行控制).同上讨论,用动态随机数命令 “=NORMINV(RAND( ), $C$1, $E$1)”在单元格区域 B3:J1002 中生成 1000 份容量均为 9 的 (6,1)N 随机样本,在单元格区域 K3:K1002 内计算出相应的 1000 个样本均值.在单元格 L3 内用均值函数计算出这 1000 个样本均值的均值 X =5.996,这个值非常接近于总体均值 6; 在单元格 M3 内计算出这 1000 个样本均值的标准差 XS =0.336,这个值近似等于总体标准差 σ =1 的1=19=13n ,即有 3 =1.009 1=XS σ≈ .这也验证了σ =1 时定理 17.1 的结论. 图 17.3 1000 份容量为 9 的正态总体 N(6,1)样本均值的均值、标准差、频数和频率 92 图 17.4 左图给出了总体标准差σ =1 时模拟数据的柱形图,而图 17.4 右图给出了频率散 点图与 N(6,1)密度曲线的比较,它们同样验证了定理 17.1 的结论. 图 17.4 正态总体 N(6,1)样本均值柱形图及频率折线)密度曲线 讨论 上述实验过程可以从两方面进行动态演示: (1) 由于使用了动态随机数命令,因此每按一次 F9 键,所有随机数和与其相关的均值、 标准差、频数、频率、柱形图和折线图都会发生相应的变化,并且变化过程中,总是近似保 持关系式 XnSσ ≈ 成立; (2) 由于对将总体标准差σ (即单元格 E1 中的值)设置了滚动条 (在单元格 F1 位 置)对其进行滚动控制,只要单击该滚动条的左右小黑三角就会发现 E1 中的值会在 1~ 4 之 间发生大小变化,表中相关的数值和图形也会发生相应变化,并且变化过程中,也总是近似 保持关系式 XnSσ ≈ 成立. 此外,对非正态总体,也可以进行类似讨论,如下面的讨论题: 讨论题:产生 1000 份来自均匀分布总体 ~1,9XU()的简单随机样本,计算每份样本 的均值,再计算这 1000 个样本均值的均值、方差和标准差,并和总体的均值和标准差进行 比较. 93 实验十八 时间序列分析实验 18.1 实验原理 时间序列分析(Time Series Analysis)就是利用按时间先后顺序排列起来所形成的数列, 应用数理统计方法加以处理,以预测未来事物的发展.时间序列分析常用在国民经济宏观控 制、企业经营管理、市场潜量预测、气象预报、地震前兆预报、农作物病虫灾害预报等方面. 对时间序列进行预测,很重要的工作是趋势预测.移动平均和指数平滑是时间序列趋势 预测的两种常见平滑方法,以下将详细介绍移动平均和指数平滑的原理. 18.1.1 移动平均原理 移动平均又称滑动平均(Moving average,MA)是用一组最近的实际数据值来预测未 来一期或几期数据的一种常用方法.其基本思想是:根据时间序列逐项推移,依次计算包含 一定项数的序时平均值,以反映长期趋势.所平均的范围可以是整个序列(整体平均数), 也可以是序列中的一部分(局部平均数).当原始数据既不快速增长也不快速下降,且不存 在季节性因素时,移动平均法能有效地消除预测中的随机波动.具体计算公式为: 设某一时间序列为 12,,,tyy y ,则 1t + 时刻的预测值为 1 11 1 0 1ˆ . n tt tn ttj j yy yyynn − −−+ +− = +++==∑ 式中: n 称为进行移动平均计算的过去期间的个数. 18.1.2 指数平滑原理 指数平滑是加权移动平均的进一步发展和完善,最初由美国经济学家布朗(Robert G..Brown)于 1959 年首次提出.所谓平滑就是通过某种平均方式消除历史统计序列中的随机 波动.指数平滑法是通过计算指数平滑值,配合一定的时间序列预测模型对现象的未来进行 预测.其原理是任一期的指数平滑值都是本期实际观察值与前一期指数平滑值的加权平均. 根据平滑次数不同,指数平滑法一般可分为:一次指数平滑、二次指数平滑和三次指数 平滑等. 一、 次指数平滑 当时间数列无明显的趋势变化,可用一次指数平滑进行预测.一次指数平滑以预测目标 的本期实际值和本期预测值为基数,分别给予不同的权数求出指数平滑值作为未来的预测 值.其指数平滑计算公式可表示为 )1( 1 )1( )1( −−+= ttt SYS αα . 其中 (1) tS 和 (1) 1tS − 分别为 t 期和 1t − 期的一次指数平滑值, (0 1)α α

  , 1ρ ≤ . 由二维正态分布的性质知,当 X 与Y 独立等价于 X 与Y 的相关系数 0ρ = .或者说, 当 0ρ ≠ 时 X 与 Y 不独立. 19.2 实验目的及要求 掌握在 Excel 中产生分量相互独立的二维正态分布随机数;理解产生分量不独立的二维 正态分布随机数的原理和方法. 19.3 实验过程 例 19.1 产生服从二维正态分布 (7,1;6,1;0)N 的随机向量(X, Y). 解 若随机变量 X 与 Y 相互独立(即 0ρ = )且 2 11~(,)XNμ σ , 2 22~(, )YNμ σ ,则 22 11 2 2(,)~(, ; , ;0)XY Nμσ μσ , 于是只要在 Excel 工作表中 A 列上产生 2 11(, )N μ σ 随机数 X,在 B 列上产生 2 22(, )N μ σ 随机 数 Y,则由 X 和 Y 组成的二维随机向量(X, Y)就服从二维正态分布 22 11 2 2(, ; , ;0)N μσ μσ .用 这种方法,只要由随机数发生器分别产生 (7,1)N 随机数 X 和 (6,1)N 随机数 Y,则由 X 和 Y 组成的二维随机向量(X, Y)就服从二维正态分布 (7,1;6,1;0)N . 101 例 19.2 产生服从二维正态分布 (7,1;6,1;0.6)N 的随机向量(X, Y). 解 此时情形要复杂得多.先讨论一般情形,设 n 维正态分布随机向量 1(,, )nXX X′= … 的联合密度函数为 11 2 1 22 ()() 1 1() ( , , ) (2 ) n xx nfx fx x e μ μ π −′−Σ −== Σ , 其中 1(, , )nx xx′= ; 1(, , )nμ μμ′= 是 X 的均值向量, ()ij n nσ ×Σ = 是 X 的协方差阵,这 里 0Σ

  为正定阵, [( )( )]ij i i j jEX Xσ μμ=− −.由于 Σ 为正定阵,故存在下三角阵 C, 使得 CC′Σ= ;若设 1(,, )nUU U′= ,U 的各个分量相互独立均服从 (0,1)N 分布,那么 可以证明 X CUμ= + 服从以 1(, , )nμ μμ′= 为均值向量,以 CC′Σ = 为协方差阵的 n 维正态分布. 由上述讨论可得产生分量相关的 n 维正态分布随机数的方法如下: 首先产生独立同 (0,1)N 分布的随机变量 1,,nUU ;接着对正定阵 Σ 作 Cholesky 分解 (即下三角×上三角分解): CC′Σ= ;其中 11 21 22 12 00 0 nn nn c ccC cc c ⎛⎞ ⎜⎟ ⎜⎟= ⎜⎟ ⎜⎟ ⎝⎠ , 令 1 n kk kii i X cUμ = =+∑ ,(1,,)kn= ,这样就得到服从 (,)nN μ Σ 的 n 维正态分布随机向量 1(,, )nXX X′= . 由于二维正态分布 (7,1;6,1;0.6)N 的协方差阵 Σ 可作如下分解: 2 112 2 12 2 1 0.6 1 0 1 0.6 0.6 1 0.6 0.8 0 0.8 σρσσ ρσ σ σ ⎛⎞⎛⎞⎛⎞⎛⎞Σ= = =⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝⎠⎝⎠ , 故得到变换公式 710 7 6 0.6 0.8 6 0.6 0.8 XUU YVUV +⎛⎞⎛⎞⎛ ⎞⎛⎞⎛ ⎞=+ =⎜⎟⎜⎟⎜ ⎟⎜⎟⎜ ⎟++⎝⎠⎝⎠⎝ ⎠⎝⎠⎝ ⎠ . 操作过程如图 19.1. 102 图 19.1 产生相关二维正态分布随机数 例 19.3 设随机变量 X 与Y 相互独立,且 11~(,)XNμ σ , 22~(,)YNμ σ ,则 ()()22 12 22 12 1 2 12 1(, ) () () , ( , )2 xy XYfxy f xf y e xy μμ σσ πσ σ ⎡⎤−−⎢⎥−+ ⎢⎥⎣⎦== −∞

  0,COUNTIF(E:E,E5)

  1),E5)”,按住填充柄向下拖动到 K33,条件区域定 义完毕. 依次点击菜单栏中【数据】/【筛选】/【高级】按钮调出【高级筛选】窗口,选择“将筛 选结果复制到其他位置”,勾选“选择不重复的记录”. 在“列表区域”文本框中填入“$E$4:$E$33”. 在“条件区域”文本框中填入“$H$4:$H$33”. 在“复制到”文本框中填入 “$I$4”. 按下确定按钮即可完成筛选工作,E 列中所有重复两次以上且大于 0 的 的不同数值 显示在 I5:I9. 在 J5 输入公式“=COUNTIF(E:E,=&I5)”,按住填充柄下拉至 J9,完成对上述 的不同 数值的频数统计. (5)在单元格 M6 中键入公式“=SUMIF(F:F,

  0)”,计算正 的秩和(简称正秩和). 108 (7)在单元格 M8 中填入公式“=COUNT(E:E)-COUNTIF(E:E,=0)”,计算非 0 的 的 个数 n. (8)在单元格 I9 中键入公式 “=(ABS(M7-M8*(M8+1)/4)-0.5)/SQRT((M8*(M8+1)*(2*M8+1)/24)-SUM((J5:J9)^3-(J5:J9))/48 )”计算 W 统计量的值. (9)在单元格 M10 中键入公式“=(1-NORMSDIST(M9))*2”计算 p 值. 由于 p 值非常小,从而我们有很大把握判断方法 2 优于方法 1. 检验完成后的工作表见图 20.2. 图 20.2 阅读能力测验数 据的 Wilcoxon 符号-秩检验 20.4 讨论 统计实践中的数据,根据其所含信息的丰富程度从低到高可分为四个等级,高级数据 一定可以视为低级数据,反之则不然. (1)分类数据(nominal scale data)——如果数据仅仅表示类别,而无其他含义.如用 1 表示男性,0 表示女性,这里的 1 和 0 就是类别数据,对其进行运算或比较大小均无意义. (2)等级数据(ordinal data)——如果数据仅仅表示一定的顺序或等级,而无指定数 值或其数值并无数量上的含义.如我们用 3、2、1、0 表示学生的某种能力的优、良、中、差 四个等级,这里的 3、2、1、0 仅仅表示能力的高低,从而对其比较大小是有意义的,但是 其数值没有数量上的意义,0 并非意味没有能力,3 与 2 的差和 1 与 0 的差虽然都是 1,但是 这两个 1 显然含义是不同的. (3)等距数据(interval scale data)——如果数据除了表示等级,其差距还有明确的数 量含义,且相同差距表示相同的含义.如温度数据就是典型的等距数据,35 摄氏度与 30 摄氏 109 度的差距和 5 摄氏度与 0 摄氏度的差距含义是一样的,但是要注意的是说“30 度是 15 度的 两倍”是没有意义的,即两个数据的比是没有意义的,因为温度测量中的 0 度并非意味着没 有温度,即在等距数据中没有绝对的 0. (4)比率数据(ratio data)——如果数据既是等距数据又有绝对(或固定)的 0 点.如 体重、身高等数据就是典型的比率数据.对比率数据而言,任何两个数据的比是有意义的, 如 30 公斤是 15 公斤的两倍. 对于比率和等距数据,均值和标准差都是有意义的. 如果教育测验设计得足够精细(这通常是难以做到的),所得分数经过适当变换仍然可 以看成比例甚至比率数据,从而常用的针对比率数据的检验方法如 t 检验等仍然可以使用. 本实验中进行高级筛选时,可直接在在“列表区域”中填入“$E$4:$E$33”,也可先单击文本 框后面的折叠按钮,然后在工作表中选取相应区域完成.特别注意:这里选取的是 E 列第 4 至 33 行,E 的第 4 行内容是“d”,而非数据,但是由于是该行的开头字段,必须包含在“列表 区域”中,否则将出现错误“提取区域中字段名丢失或非法”.事实上,如果“列表区域”的第一 个单元格的值和结果将要复制到的区域的第一个单元格的值不相同就会出现这个错误.另外 两个文本框的填写同“列表区域”. 请读者仿例 20.1 自己搜集数据,用 Excel 完成 Mann-Whitney U 检验. 110 实验二十一 含有虚拟变量的多元线 实验原理 在实际建模过程中,因变量(被解释变量)不仅受到定量变量影响,通常还受定性变量 影响.例如需要考虑性别、民族、文化程度、季节差异、不同时间段等因素的影响.这些因 素也应该包括在模型中. 由于定性变量通常表示的是某种特征的有和无,所以量化方法可采用取值为 1 或 0.这 种变量称作虚拟变量(dummy variable),用 D 表示.例如中国成年人体重 y(kg)与身高 x(cm)的回归关系如下: 105, 1 ( ),100 5 100, 0 ( ). xDyxDxD −=⎧=− + − =⎨ −=⎩ 男 女 这里 D 取 1 或 0 就分别表示男性和女性.虚拟变量可以看成在模型中新增加的自变量,对 其回归系数的估计与检验与普通的定量变量相同. 注意: ① 若定性变量含有 m 个类别,应引入 m-1 个虚拟变量,否则会导致多重共线性,称作 虚拟变量陷阱(dummy variable trap). ② 关于定性变量中的哪个类别取 0,哪个类别取 1,是任意的,不影响检验结果. ③ 定性变量中取值为 0 所对应的类别称作基础类别(base category). ④ 对于多于两个类别的定性变量可采用设一个虚拟变量而对不同类别采取赋值不同的 方法处理.比如可用虚拟变量 D 赋值 1、0 和-1 分别表示文化程度为大学、中学和小学. 21.2 实验的目的和要求 实验目的是学习并初步掌握在 Excel 中建立含有虚拟变量的多元线性回归模型的方 法.具体要求为:对含有属性变量的实际数据,在 Excel 中引入哑变量建立包含虚拟变量的 多元线性回归模型;对模型进行改进优化,并把建立的模型与不含虚拟变量的回归模型进行 比较;利用所建立的模型进行预测. 21.3 实验过程 例 21.1 图 21.1 给出了 1996 年 1 季度至 2000 年 4 季度我国季节 GDP 数据(数据来源 于《中国统计年鉴》1998-2001).要求: (1) 用软件 Excel 建立 GDP 关于时间 t 以及季度虚拟变量 D1,D2 和 D3 的多元线性回 归模型,其中 D1,D2 和 D3 是表征 1 季度,2 季度和 3 季度的虚拟变量,即 1 (1 季度数据), 1 (2 季度数据) 1 (3 季度数据) D1 = D2 = D3 = 0 (2, 3,4 季度数据). 0 (1, 3, 4 季度数据) 0 (1, 2, 4 季度数据). (2) 若只是建立 GDP 关于时间 t 的一元线性回归模型,效果如何? 111 (3) 利用上面建立的模型预测 2001 年 1 季度至 4 季度 GDP 的值. 图 21.1 1996 年 1 季度至 2000 年 4 季度我国季节 GDP 数据 解 (1) 利用 Excel 中的回归分析工具可建立的含有季度虚拟变量 D1,D2 和 D3 的多 元线 含有季度虚拟变量 D1,D2 和 D3 的多元线 所以 GDP 关于时间 t,D1,D2 和 D3 的的多元线GDP t D D D= +− − − . (21.1) 标准误差为 0.05,比较小.可决系数 2R =0.9863,接近于 1,拟合优度较高.诸回归系数均 通过了 t 检验,都很显著,这由它们对应的 p 值均很小可以看出.方差分析表中 F 统计量的 p 值也非常小:8.801E-14. 说明多元回归模型整体是非常显著的. (2) 若只是建立 GDP 关于时间 t 的一元线),可以使用数据分析工具, 也可以利用 C 列中的时间 t 数据和 B 列中的 GDP 数据画带数据点的散点图,再向该散点图 添加趋势线的方法,并且选中“显示公式”和“显示 R 平方值”(参见实验十二).得到的 一元线yx=+,如图可决系数仅为 2 0.399R = ,比上面的 2R =0.9863 小得多,说明拟合优度不高. 图 21.3 向散点图添加趋势线) 利用上面建立的含有季度虚拟变量 D1,D2 和 D3 的多元线 季度 GDP 的值. 在单元格 B22 中输入命令“=$B$44+$B$45*C22+B46*D22+$B$47*E22+$B$48*F22”, 再拖放填充至 B25,则可得到 2001 年 1 季度至 4 季度 GDP 的预测值值. 图 21.4 利用含有季度虚拟变量 D1,D2 和 D3 的多元线 年 1 季度至 4 季度 GDP 的预测值的图形展。

http://hostgladjens.net/shuzhouyufa/212.html
锟斤拷锟斤拷锟斤拷QQ微锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷微锟斤拷
关于我们|联系我们|版权声明|网站地图|
Copyright © 2002-2019 现金彩票 版权所有