1 大学数学实验 Experiments in Mathematics 实验 10 数据的统计与分析 清华大学数学科学系 数据的统计与分析的两类方法 第一类：一般意义的统计（普查） 对生产的全部 1000件产品逐一检验，发现 18件次品 对全区居民逐一调查，得到月平均支出为 828元 结果分析： 次品率： 1.8%；月平均支出为 828元 优点： 结果完全确定，可信 缺点：调查、收集的数据量可能很大，经费投入 大，有些产品不允许全部检验，如灯泡、电器的寿 命等 缺点：结果是随机的，是否可信？ 第二类：数理统计（抽查） 全部产品中随机抽取 100件，发现 2件次品 随机调查了 200位居民，得到月平均支出为 788元 结果分析： 次品率： 2%；月平均支出 788元 优点：调查、收集的数据量小，经费投入小，适合 不允许全部检验的产品，如灯泡、电器的寿命等 任务：怎样用它来估计整体的状况（全部产品的 次品率，全体居民的月平均支出） 数据的统计与分析基本内容 2.数据的整理和描述 3.随机变量的概率分布及数字特征 1.实例及其分析 4. 用随机模拟计算数值积分 5. 实例的建模和求解 1. 实例及其分析 实例1： 报童的利润 报童每天从发行商处购进报纸零售，晚上将没有卖掉的报纸退 回。如果每份报纸的购进价为 a，每份报纸的零售价为 b，每份 报纸的退回价（发行商返回报童的钱）为 c, 且满足 b≥ a≥ c。每 天报纸的需求量是随机的。为了获得最大的利润，该报童每天 应购进多少份报纸？ 假定 a=0.8元， b=1元， c=0.75元，为报童提供最佳决策。 2815203532221393天数 280 — 260 — 240 — 220 — 200 — 180 — 160 — 140 — 120 — 100 — 需求 量 159天报纸需求量的情况 2 分析： 每天报纸的需求量随机，报童每天的利润也是随机 的。只能以长期售报过程中每天的平均利润最大为目标，确定 最佳决策。 数学模型近似： 决策变量： 报童每天购进报纸的份数 n 可以通过历史数据得到每天需求量为 r的天数所占的百分比， 记做 f(r) ，如需要 200份所占的百分比为 35/159=22% 平均利润： V(n) ∑∑ ∞ = ? = ?+????= nr n r rfnabrfrncarabnV )(])[()()])(()[()( 1 0 实例1： 报童的利润（续） 实例2： 路灯更换策略 管理部门： 不亮灯泡，折合计时进行罚款。 路政部门 : 路灯维护 条件 : 需要专用云梯车进行线路检测和更换灯泡 向相应的管理部门提出电力使用和道路管制申请 向雇用的各类人员支付报酬等 更换策略 : 整批更换 路政部门的问题： 多长时间进行一次灯泡的全部更 换，换早了，很多灯泡还没有坏；换晚了，要承受 太多的罚款。 2. 数据的整理和描述 ? 数据的收集和样本的概念 ? 数据的整理、频数表和直方图 ? 统计量 ? MATLAB命令 数据的收集 某银行为使顾客感到亲切以吸引更多的资金，计划对柜台的高 度进行调整。银行随机选了 50名顾客进行调查，测量每个顾客 感觉舒适时的柜台高度，表 2为得到的数据。银行怎样依据它 确定柜台高度呢？ 50顾客感觉舒适高度（单位：厘米） 110115128124134119116130112115 12010813112212114010911795105 119125119127119110106114117118 12212112210211693108115113126 991191209510010497136110100 基本概念 ? 总体 --研究对象的全体。如所有顾客感觉舒适的高度 ? 个体 --总体中一个基本单位。如一位顾客的舒适高度 ? 样本 --若干个体的集合。如 50位顾客的舒适高度 ? 样本容量 --样本中个体数。如 50 顾客群体的舒适高度 ~随机变量 X，概率分布 F(x) n位顾客的舒适高度 { x i , i= 1,…n} (样本 )~相互独 立的、分布均为 F(x)的一组随机变量。 样本 ：随机取值的一组数据； 一组相互独立的、同分布的随机变量。 样本 ——统计研究的主要对象 数据的整理 比较直观，比较清晰的结论 21—50岁的中青年患者大约占总发病人数的 3/4，提醒民众中青年是易感人群。 100%17.77%17.50%20.14%35.69%7.64%1.27%比例 189733733238267714524人数 总数51岁以 上 41-50岁31-40岁21-30岁11-20岁10岁以 下年龄 北京地区SARS患者的统计数据（截至2003年5月5日） 3 频数表和直方图 将数据的取值范围划分为若干个区间，统计这组数据在 每个区间中出现的次数，称为 频数 ，得到一个 频数表 。 柜台高度频数表 22451286344 频数 137.65132.95128.25123.55118.85114.15109.45104.7 5 100.0595.35 中点 推测出总体的某些简单性质。如表 6表明选择柜台高度 在 107.10至 125.90的有 31人，占总人数的 62%，柜台高 度设计在这个范围内，会得到大多数顾客的满意。 直方图 (histogram)，或频数分布图 90 95 100 105 110 115 120 125 130 135 140 0 2 4 6 8 10 12 柜台高度直方图 统计量 平均值 (mean，简称样本均值 )定义为 频数表和直方图给出某个范围的状况，无法直接给 出具体值，如例 1关于确定柜台高度的问题 ∑ = = n i i x n x 1 1 26.115=x 可作为设计柜台高度的参考值 两个班的一次考试成绩 787783868082828177907385乙班 81995390857993558788687978888669甲班 32313029282726252423222120191817序号 79848578837578849082818285828384乙班 80948788739365888784799592858892甲班 16151413121110987654321序号 现象：甲班的平均值： 82.75分，乙班的平均值： 81.75分 结论：大致表明甲班的平均成绩稍高于乙班 现象：甲班中 90分以上的有 7人，但有 2人不及格，分数比较分 散。乙班全在 73分到 90分之间，分数相对集中 40 60 80 100 0 2 4 6 8 10 12 14 40 60 80 100 0 2 4 6 8 10 12 14 16 18 为了描述数据的这种分散程度（统计上称为变异）， 统计上引入标准差的概念。 样本 x=(x 1 , x 2 , … , x n )的 标准差 (Standard deviation)为： 2/12 1 ])( 1 1 [ xx n s n i i ? ? = ∑ = 甲班的标准差为 10.98分，乙班的标准差为 3.98分， 表明甲班成绩的分散程度远大于乙班。 统计量 ：由样本加工出来的、集中反映样本数量特 征的函数。 三类统计量：表示 位置的 ，表示 变异程度的 ，表示 分布形状的 。 4 3 1 3 1 )( 1 xx ns g n i i ?= ∑ = 4 1 42 )( 1 xx ns g n i i ?= ∑ = 偏度 (skewness)：分布对称性 峰度 (kurtosis )：分布形状 表示位置的还有： 中位数 (median)：将数据由小到大排序后处于中间位 置的那个数值。 当样本容量 n为奇数时，中位数唯一确定；当 n 为偶数时，定义为中间两个数的平均值。 表示变异程度的还有： 极差 (range)： x 1 , x 2 , … , x n 的最大值与最小值之差。 方差 (variance)：标准差的平方 s 2 。 表示分布形状的： MATLAB数据描述的常用命令 峰度 g 2 同上峰度kurtosis(x) 偏度 g 1 同上偏度skewness(x) var(x,1)：同上方差 s 2 同上方差var(x） std(x,1): (3)式 中 n-1改成 n 标准差 s同上标准差std(x） 极差同上极差range(x) 中位数同上中位数median(x) x: 原始数据行向量均值mean(x) 同上直方图同上直方图hist(x,k) [n,y]=hist(x)中 k 取缺省值 10 n: 频数行向量 y: 区间中点行向量 x: 原始数据行向量 k：等分区间数 频数表[n,y]=hist(x,k) 注意事项输出输入名称 命令 求银行柜台高度的频数表、直方图及均值等统计量： X =[100 110 136 97 104 100 95 120 119 99 ... % 输入表 2数据， ...为延续符号 126 113 115 108 93 116 102 122 121 122 ... 118 117 114 106 110 119 127 119 125 119 ... 105 95 117 109 140 121 122 131 108 120 ... 115 112 130 116 119 134 124 128 115 110]; [N,Y]=hist(X)， % 频数表 hist(X), % 直方图 x1=mean(X)， x2=median(X) % 各个统计量 x3=range(X)， x4=std(X) x5=skewness(X)， x6=kurtosis(X) 示例 exam1001a.m 输出图和下列结果： N = 4 4 3 6 8 12 5 4 2 2 Y= 95.3500 100.0500 104.7500 109.4500 114.1500 118.8500 123.5500 128.2500 132.9500 137.6500 x1 = 115.2600， x2 =116.5000 x3 =47， x4 =10.9690 x5 = -0.0971， x6 =2.6216 exam1001b.m 3. 随机变量的概率分布及数字特征 ? 频率与概率 ? 概率密度与分布函数 ? 期望和方差 ? 常用的概率分布 ? MATLAB命令 频率与概率 在保证抽取样本的随机性和独立性，当样本容量无 限增大时，频率会趋向一个确定值，这个值称为随 机变量 X落入区间（ a,b]的 概率 （ Probability），记 作 频率 : 样本数据在一个确定区间（ a,b]的频数 k与样本 容量 n的比值 n k bXaf =≤< )( )( bXaP ≤< 90 95 100 105 110 115 120 125 130 135 140 0.24 p(x) 5 概率密度与分布函数 dxxpbXaP b a ∫ =≤< )()( 概率密度函数 (Probability density function，简称 概率 密度 ) ： 0)( ≥xp ∫ ∞ ∞? = 1)( dxxp 概率分布函数 （ Cumulative distribution function，简称 分布函数 ） ∫ ∞? =≤= x dxxpxXPxF )()()( 1)(,0)( =∞=?∞ FF )()(}{ aFbFbXaP ?=≤< dx dF xp =)( 对于连续随机变量 期望和方差 随机变量 X的期望就是平均值的意思，记作 EX或 μ ∫ ∞ ∞? = dxxxpEX )( ∫ ∞ ∞? = dxxxpEX )( ∫ ∞ ∞? ?= dxxpEXxDX )()( 2 ExEx n xE n i i == ∑ =1 1 n Dx Dx n xD n i i == ∑ =1 2 1 常用的概率分布 均匀分布 (Uniform distribution) ： X~U (a,b) ? ? ? ? ? ∈ ? = 其他。 ,0 ],,[ , 1 )( bax ab xp 12 )( , 2 2 ab DX ba EX ? = + = 指数分布 (Exponential distribution)： X~Exp (λ) ? ? ? ? ? ≥ = ? 其他 ,0 0 , 1 )( xe xp x λ λ 2 , λλ == DXEX -1 0 1 2 3 4 5 6 0 0.1 0.2 0.3 0.4 0.5 U(0,2) U(1,5) 图5 均匀分布概率密度函数图形 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 Exp(2) Exp(4) 正态分布 (Norm distribution)： ) 2 )( exp( 2 1 )( 2 2 σ μ σπ ? ?= x xp 2 , σμ == DXEX -6 -4 -2 0 2 4 6 0 0.1 0.2 0.3 0.4 N(0,1) N(0,2 2 ) ),(~ 2 σμNX 2 χ 分布 （ Chi square）： ∑ = = n i i XY 1 2 )(~ 2 nY χ ， n称 自由度 0 5 10 15 20 0 0.05 0.1 0.15 0.2 Chi2(5) Chi2(10) 定义： 服从标准正态分布的随机变量 n XXXL,, 21 其中 相互独立、 6 )(~ ntT n为自由度t分布 （ Student分布） -6 -4 -2 0 2 4 6 0 0.1 0.2 0.3 0.4 t(2) t(20) nY X T / = )1,0(~ NX )(~ 2 nY χ 其中， X,Y相互独立定义： F分布： ),(~ 21 nnFF ),( 21 nn 称自由度 0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 F(10,5) F(10,50) 2 1 / / nY nX F = )(~ 1 2 nX χ )(~ 2 2 nY χ 相互独立定义： 一次实验只有两种结果：成功和失败，记成功的概率为 p， q=1-p， n次独立实验中成功的次数是随机变量 X nkqp k n kXP knk ,,1,0,)(L= ? ? ? ? ? ? ? ? == ? npqDXnpEX == , 当二项分布的 n→∞， np→λ（常数）时 L,2,1,0, ! )( === ? ke k kXP k λ λ λλ == DXEX , 二项分布 （ Binomial distribution） X~B (n, p) 泊松分布 （ Poisson distribution） X~Poiss(λ), 背景问题： 产品检验中的废品个数 背景问题：服务系统在一定 时间内接到的呼唤数 MATLAB命令 2 χ poissbinoft chi2 normexpunif字符 泊松分布二项分布F分布t分布分布正态分布指数分布均匀分布分布 rndstatinvcdfpdf字符 随机数生成均值与方差逆概率分布分布函数概率密度功能 y=normpdf(1.5,1,2) 正态分布 (μ=1, σ=2) x=1.5处的概率密 度（标准正态分布的 μ, σ 可省略） y=normcdf([-1 0 1.5],0,2) 在 x= -1, 0, 1.5处的分布函数)2,0( 2 N [m,v]=fstat(3,5) 计算 F(3,5)的期望和方差 4. 用随机模拟计算数值积分 ? 4.1 定积分的计算 ? 4.2 重积分的计算 ? 4.3 MATLAB实现 方法的直观解释 ——随机投石 y 1 0 1 x · 向单位正方形里随机投 n块小石头 * * * * * * * 若有k块小石头落在1/4单位圆内，当n很大时 1/4单位圆的面积 n k ≈ 4 π （ 计算 π的一种方法 ） 1）随机投点法 目的：计算1/4单位圆的面积 7 大数定律（贝努利定理） 设 k是 n次独立重复试验 中事件 A发生的次数。 p是事件 A在每次试验中发生 的概率，则对任意的正数 ε, 有 1)|(|lim =<? ∞→ εp n k P n 随机变量(X,Y)在单位 正方形内均匀分布 1,0,1),( ≤≤= yxyxp n k p ≈= 4 π 点( x i , y i )落在1/4单位圆内概率 2 1 ii xy ?≤即满足 y 1 0 1 x · 2 1 ii xy ?≤ 一般地 投点坐标 (x i , y i ),i=1,2,…n, x i , y i 是相互独立、 (0,1)内均匀分布的随机变量（ (0,1)随机数） ∫∫ ? =?∈ dxdyyxpYXP ),()),(( 1)(0 ,10: ≤≤≤ ≤≤? xfy x 产生 n组 (0,1)随机数 (x i ,y i ),其中 k 组满足 )( ii xfy ≤ n k YXP ≈?∈ )),(( 1)(0,)( 1 0 ≤≤≈ ∫ xf n k dxxf 随机投点法 y 1 0 1 x · )(xfy= ? ∫∫∫ =?= 1 0 1 0 )( 0 )(1 dxxfdydx xf 大数定律（辛钦定理） 设随机变量 n YYY ,,, 21 L 相互独立，服从同一个分布，且具有数学期望 ),,2,1( niEY i L== μ 则对任意的正数 ε有 1)| 1 (|lim 1 =<? ∑ = ∞→ εμ n i i n Y n P 随机变量 X 的概率密度为 bxaxp ≤≤),( )(XfY= 的期望为 ∫ = b a dxxpxfXfE )()())(( 2) 均值估计法 产生 (0,1)随机数 x i (i=1,2,…n), n很大 ∑ ∫ = ≈ n i i xf n dxxf 1 1 0 )( 1 )( 用随机模拟方法计算任意区间上的积分 = ∫ b a dxxf )( ∫ ?+? 1 0 ))(()( duuabafab uabax )( ?+= ))(( 1 ∑ = ?+ ? ≈ n i i uabaf n ab 其中 u i 为(0,1)随机数 均值估计法 ? 不要产生 y i , 不用比较 )( ii xfy ≤ 1)(0 ≤≤ xf 限制；? 没有 均值估计法 的优点 用随机模拟法 计算重积分 ∫∫ ? dxdyyxf ),( 1)( )(0 ,10: 2 1 ≤≤ ≤≤ ≤≤? xg yxg x x y 1 0 1 g 2 (x) g 1 (x) ? ),( 1 ),( 1 k m k k yxf n dxdyyxf ∑ ∫∫ = ? ? 产生相互独立的(0,1)随机数 x i ,y i, ,i =1, …n, 落在 ?内的 m个点记作( x k ,y k ),k= 1, …m ? 可用于任意的 f, ?, 且可推广至高维 ? 结果的精度和收敛速度与维数无关 ? 计算量大，精度低，结果具有随机性 重积分的计算 dxgyxgcbxadxdyyxf ≤≤≤≤≤≤? ∫∫ ? )()(,:,),( 21 ),1(, niyx ii L= 分别为 [a, b]和 [c, d] 区间上的 均匀分布随机数，判断每个点是否落在 Ω 域内， 将落在 Ω 域内的 m个点记作 mkyx kk L,1),,( = 则 ),( ))(( ),( 1 k m k k yxf n cdab dxdyyxf ∑ ∫∫ = ? ?? ? 8 MATLAB实现 随机数的产生 ： unifrnd(a,b,m,n) 产生 m行 n列[ a,b]区间上的均匀分布随机数。当 a=0, b=1时，可用 rand(m,n) 随机投点法计算 π n=10000; x=rand(2,n); k=0; for i=1:n if x(1,i)^2+x(2,i)^2<=1 k=k+1; end end p=4*k/n Cal_pi.m 5.实例的建模和求解 ? 报童的利润 ? 路灯更换策略 报童的利润 1） 每份报纸的购进价为 a，每份报纸的零售价为 b， 每份报纸的退回价（发行商返回报童的钱）为 c。 假设： 2）需求为连续随机变量 x，大致服从正态分布 3）将历史的统计表看作需求量的频率，由此可以计算需求 量的均值和标准差 报童每天的平均利润 V(n) {} ∫∫ ∞ ?+????= n n dxxnpabdxxpxncaxabnV )()()())(()()( 0 ) 2 )( exp( 2 1 )( 2 2 σ μ σπ ? ?= x xp 其中 μ和 σ由3）的假设 计算得到 0)()()()( )()()()()()()()()( 0 0 =?+??= ?+?????=′ ∫∫ ∫∫ ∞ ∞ n n n n dxxpabdxxpca dxxpabnnpabdxxpcannpabnV ca ab dxxp dxxp n n ? ? = ∫ ∫ ∞ )( )( 0 cb ab dxxp n ? ? = ∫ ∞? )( 简化为 定性分析 在b≥a ≥c的条件下讨论a、b、c的变化 对最佳决策 n的影响。 1) 当b>a=c时，即购进价与退回价相同且零售价高于 购进价，报童不承担任何卖不出去的风险，他将从发 行商处购进尽可能多（无穷多）的报纸。这样必然造 成发行商的损失。 2) 当b=a>c时，即零售价与购进价相同且高于退回 价，报童无利润可得，他不从发行商处购进任何报纸 (n<<0)。这样发行商也无法获得任何利益。 cb ab dxxp n ? ? = ∫ ∞? )( 3) 当b>a>c时，发行商和报童才能获得利益，报童将 根据需求量的随机规律制定自己的应对策略。b-a越 大，购进量 n越大；b-c越大，购进量 n越小。这些都 符合直观的理解。 定量求解 根据 159天报纸需求量的分布情况 2815203532221393天数 280 — 260 — 240 — 220 — 200 — 180 — 160 — 140 — 120 — 100 — 需求量 a=0.8元， b=1元， c=0.75元 a=0.8;b=1;c=0.75; q=(b-a)/(b-c); r=[3 9 13 22 32 35 20 15 8 2]; rr=sum(r); x=110:20:290; % 需求量取表中小区间的中点 rbar=r*x'/rr % 计算均值 s=sqrt(r*(x.^2)'/rr-rbar^2) % 计算标准差 n=norminv(q,rbar,s) % 用逆概率分布计算n 计算结果：rbar =199.4340 s =38.7095 n =232.0127 报童应购进 2 3 2份报纸 Newsboy.m 9 路灯的更换问题 假设 1) 每个灯泡更换价格： a ——灯泡的成本和安装时分摊到每个灯泡的费用 2) 不亮灯泡单位时间罚款： b 3) 假定灯泡寿命服从 4) 更换周期 : T 5）灯泡总数： K ),( 2 σμN 模型： T dxxpxTKbKa TF T ∫ ∞? ?+ = )()( )( 0= dT dF b a dxxxp T = ∫ ∞? )( 计算： 可得： 定性结果分析 结论 1： a/b越大，更换价格与惩罚费用之比越大， 更换周期 T应越长 结论 2：若以灯泡的平均寿命为更换周期， 惩罚费用为： ∫ ∞? = μ dxxxp a b )( 考虑灯泡寿命服从 π σμ 22 ? = a b 具体示例 某品牌灯泡服从平均寿命为 4000小时，标准差为 100小时的正态分布，每个灯泡的安装价格为 80元， 管理部门对每个不亮的灯泡制定的惩罚费用为 0.02元 ／小时，求最佳更换周期。 则： 计算结果： 最佳更换周期为 4459（小时） 不同惩罚费用对更换周期的影响 37153797382839183977T 1010.50.10.05b ),( 2 σμN a=80; b=0.02; aoverb=a/b mu=4000; s=100; t=mu; % 设定 T的初值 step=0.1; % T增加或减少的步长 var=0.01; %（ 50）左右端比较的误差限 vp= mu*normcdf(t,mu,s)-s^2*normpdf(t,mu,s); % 计算（ 50）左端 if vp>aoverb % （ 50）左端大于右端， T减少 while (vp-aoverb)>var t=t-step; vp= mu*normcdf(t,mu,s)-s^2*normpdf(t,mu,s); end end if vp<aoverb % （ 50）右端大于左端， T增加 while (aoverb-vp)>var t=t+step; vp= mu*normcdf(t,mu,s)-s^2*normpdf(t,mu,s); end end vp， t 计算灯泡更换周期的程序 Lamp.m 布置实验 目的 1)掌握数理统计的基本概念； 内容课上布置，或参见网络学堂 2) 掌握用随机的方法计算积分； 3) 掌握对实际问题的建立概率模计和进行计算