E N D
第九章 概率统计计算 北京交通大学
9.1 概率统计软件包 • Mathematica可以处理概率统计方面的计算,有关的命令都在Mathematica自带的统计软件包中, 这些软件包存放在Mathematica系统自己带有程序包,存放在C:\wnmath22\Packges\Statisti目录中,用户可以在Mathematica的工作窗口键入Ctrl+ O,调出Open窗口, 将该窗口左下脚的文件类型选为Packages (*.m), 并用鼠标双击文件夹packages打开其中的子文件夹,然后任意双击Statisti文件夹, 就可以在窗口左上部分看到很多以.m为扩展名的Mathematica所有自带的概率统计软件包文件: (见图) 下一页
Mathematica中的部分概率统计软件包文件名,调用名称及涉及的问题Mathematica中的部分概率统计软件包文件名,调用名称及涉及的问题
9.2 Mathematica概率统计软件包中最常用的命令 • 为了使用的方便,下面写出一些概率统计软件包中最常用的内容及其调用文件名 需调用Statistics`DescriptiveStatistics`软件包才能使用的函数: Mean[data] 计算样本数据data的均值 Median[data] 计算样本数据data的中值 Variance[data] 计算样本数据data的方差 StandardDeviation[data] 计算样本数据data的标准差 • 注意:data是由离散数据组成的表
例1:1) 已知样本数据为dat={3.2,5.1,1,4,2},试计算dat的均值、中值、方差、标准差。 2) 产生[0,1]上的20个随机实数,并计算它们的均值、中值、方差、标准差。 解: In[1]:= <<Statistics`DescriptiveStatistics` *调用统计软件包 In[2]:=dat={3.2, 5.1, 1, 4, 2}; In[3]:=Mean[dat] Out[3]:=3.06 In[4]:=Median[dat] Out[4]:=3.2 In[5]:= Variance[dat] Out[5]:= 2.608 In[6]:= StandardDeviation[dat] Out[6]:= 1.61493 In[7]:=dat1=Table[Random[],{20}] Out[7]:= {0.93234, 0.439331, 0.407442, 0.469035, 0.741679, 0.884562, 0.111029, 0.696056, 0.0591917, 0.622276, 0.825287, 0.540449, 0.594691, 0.597846, 0.490196, 0.463414, 0.404672, 0.19069, 0.105273, 0.942455} In[8]:=Mean[dat1] Out[8]:= 0.525896 In[9]:=Median[dat1] Out[9]:= 0.515323 In[10]:= Variance[dat1] Out[10]:= 0.0724088 In[11]:= StandardDeviation[dat1] Out[11]:= 0.269089
需调用Statistics`DiscreteDistributions`软件包才能使用的概率分布和函数:需调用Statistics`DiscreteDistributions`软件包才能使用的概率分布和函数: • BernoulliDistribution[p] 表示均值为p的离散伯努力分布 • BinomialDistribution[n, p] 表示参数为n,p的二项分布b(n,p) • GeometricDistribution[p] 表示参数为p的几何分布 • HypergeometricDistribution[n, nsucc, ntot] 表示参数为n, nsucc, ntot 的超几何分布 • PoissonDistribution[mu] 表示参数为mu的F泊松分布 • PDF[distribution, k] 离散分布distribution的分布律P{=k} • CDF[distribution, x] 概率分布为distribution且随机变量小于值x的概率P{<x} • Mean[distribution] 计算离散分布distribution的均值 • Variance[distribution] 计算离散分布distribution的方差 • StandardDeviation[distribution] 计算离散分布distribution的标准差 • Random[distribution] 产生具有概率分布为 distribution一个伪随机数
例2: 设随机变量服从参数为0.8的泊松分布 (1)求随机变量的均值、中值、方差、标准差和分布律。 (2)求随机变量 4的概率 解:泊松分布是离散分布,故需调用处理离散概率问题的软件包,执行命令为 In[1]:= <<Statistics`DiscreteDistributions` *调用统计软件包 In[2]:=s=PoissonDistribution[0.8] Out[2]:= PoissonDistribution[0.8] In[3]:= {Mean[s], Variance[s], StandardDeviation[s] } Out[3]:={0.8, 0.8, 0.894427} In[4]:= PDF[s, k] Exp[-1*0.8] 0.8 k Out[4]:=If[!Negative[k], If[IntegerQ[k], -------------------- , 0 ], 0] k! In[5]:= 1-CDF[s,3] *因为概率P( 4)=1- P( < 4) Out[5]:= 0.00907986
例5:假设投掷一个均匀硬币只能出现正面和反面两种情况, 用Mathematica命令来验证投掷出现正面的概率为0.5。 解:设X表示投掷一个均匀硬币出现正面和反面的随机变量,它只取两个值0和1, 采用具有概率分布均值为0.5的离散伯努力分布BernoulliDistribution[0.5] 产生的伪随机数Random[BernoulliDistribution[0.5]] 来模拟实际投掷一个均匀硬币的情况,规定出现随机数是1表示投掷硬币出现正面;0 表示投掷硬币出现反面。命令中分别用产生的100个伪随机数、500个伪随机数和1000个伪随机数出现数1的频率来验证投掷出现正面的概率为0.5的结论,命令为: In[1]:= <<Statistics`DiscreteDistributions` *调用统计软件包 In[2]:= sy[n_]:=Module[{face,s}, *定义模拟函数 s=BernoulliDistribution[0.5]; For[face=0;i=1, i<=n, i=i+1, If[Random[s]==1, face=face+1] ]; N[face/n] ] In[3] = { sy[100], sy[500], sy[1000] } Out[3]={ 0.53, 0.514, 0.472 } • 从模拟试验结果可以看到投掷出现正面的概率在0.5附近波动。
需调用Statistics`ContinuousDistributions`软件包才能使用的概率分布和函数需调用Statistics`ContinuousDistributions`软件包才能使用的概率分布和函数 • BetaDistribution[, ] 表示参数为 和的Beta连续分布 • CauchyDistribution[, ] 表示参数和的柯西连续分布 • ChiSquareDistribution[n] 表示有 n个自由度的2 连续分布 • ExponentialDistribution[lambda] 表示参数为 的指数连续分布 • "FRatioDistribution[n1, n2] 表示分子参数为n1和分母参数为n2的F连续分布 • NormalDistribution[, ] 表示均值为标准差为的正态分布N (, 2) • RayleighDistribution[] 表示参数为的瑞利连续分布 • "StudentTDistribution[n] 表示有 n个自由度的t 连续分布 • UniformDistribution[min, max] 表示[min, max] 区间上的均匀分布 • PDF[distribution, x] 概率分布为distribution的分布密度函数f(x) • CDF[distribution, x] 概率分布为distribution且随机变量小于值x的概 率P{<x} • Mean[distribution] 计算概率分布为distribution均值 • Variance[distribution] 计算概率分布为distribution方差 • StandardDeviation[distribution] 计算概率分布为distribution标准差 • Random[distribution] 产生具有概率分布为 distribution一个伪随机数
例3:设随机变量服从正态分布N(0,32), (1)求出对应的分布密度函数,并画出对应的分布密度函数图形 (2)求随机变量<2的概率 解:Mathematica命令为: In[1]:= <<Statistics`ContinuousDistributions` In[2]:= dis=NormalDistribution[0,3] Out[2]:= NormalDistribution[0, 3] In[3]: =PDF[dis,x] 1 Out[3]= ------------------- x2 /18 3 E Sqrt[2 Pi] In[4]:= Plot[PDF[dis,x], {x,-10,10}, PlotRange->All ] Out[4]:=-Graphics- In[5]:= CDF[dis,2] *求随机变量<2的概率 Out[5]=0.747507
实验 1 袋内有6个白球4个黑球,从中任取两个球,求取出的两个球都是白球的概率。 分析:基本事件总数C102,有利的事件数C62,故所求概率P= C62 / C102 Mathematica 命令 In[1]:= Binomial[6,2]/ Binomial[10,2] Out[1]= 1/3 故 取出两个球都是白球的概率为1/3 2 已知在1000个灯泡中坏灯泡的个数从0到5均等可能,求从中任取100个都是好灯泡的概率。 Mathematica 命令 In[1]:= pbi =Table[1/6,{6}] ; In[2]:= pabi =Table[Binomial[1000-i,100],{i,0,5}]/ Binomial[1000,100]; In[3]: = pa =Sum[pbi[[i]]*pabi[[i]],{i,1,6}] Out[3]= In[4]:= N[pa] (*将精确结果转化为有6位有效数字的近似数*) Out[4]= 0.780693
实验 3 生成自由度为12的t分布的连续型随机变量及其概率密度函数,分布函数,并用图形显示。 Mathematica 命令 In[1]:= << Statistics`ContinuousDistributions` In[2]:= rv =StudentTDistribution[12]; In[3]:= f =PDF[rv,x] Out[3]:= (*t(12)的概率密度函数*) In[4]:=Plot[f, {x, -5, 5 } ] In[5]:= g =CDF[rv, x]; Out[5]:= (*t(12)的分布函数*) In[6]:=Plot[g,{x,-4,4}]
实验 4 某地区18岁女青年的血压(收缩压,以mm-Hg计)服从N(110,122)。在该地区任选一个18岁的女青年,测量她的血压X。求P(X≤105)和P(100<X≤120),画出血压X概率密度函数的图像。 Mathematica 命令 In[1]:= <<Statistics`ContinuousDistributions` rv =NormalDistribution[110,12]; f[x_]:=PDF[ rv, x]; F[x_]:=CDF[rv,x]; N[F[105]] Out[5]= 0.338461 (*P(X105)= 0.338461*) In[6]:= N[F[120]-F[100]] Out[6]= 0.595343 (*P(100<X100)= 0.595343*) In[7]:= Plot[ f[x], {x, 110-12*3, 110+12*3} ]
实验 5 设随机变量X-b(20-0.4),计算 (1)P{X=0} (2)P{X=1} (3)P{X<2} (4)P{X≤6} (5)P{X>10} (6)P{X≥15} Mathematica 命令 In[1]:= <<Statistics`DiscreteDistributions` brv = BinomialDistribution[20, 0.4]; f[x_]: = PDF[ brv, x]; df[x_]:= CDF[brv, x]; f[0] Out[5]= 0.0000365616 (*得P(=0) = 0.0000365616*) In[6]:= f[1] Out[6]= 0.000487488 (*得P(=1) = 0.000487488*) In[7]:= df[2]- f[2] Out[7]= 0.000524049 (*得P(<2) = 0.000524049*) In[8]:= df[6] Out[8]= 0.250011 (*得P(6 ) = 0.250011*) In[9]:= 1-df[10] Out[9]= 0.127521 (*得P( >10 ) = 0.127521*) In[10]:= 1-df[15]+f[15] Out[10]= 0.00161152 (*得P(15 ) = 0.00161152*)
实验 6 在以原点为圆心的单位圆周上任取一点,求随机地取到点的横坐标的概率密度函数。 Mathematica 命令 In[1]:= <<Statistics`ContinuousDistributions` In[2]:= rv=UniformDistribution[ -Pi, Pi ]; f =PDF[ rv, t] Out[3] = In[4]:= F=CDF[rv, t ] Out[4]= In[5]:= G[x_]:=1-ArcCos[x]/Pi; D[G[x], x] Out[5]=
实验 7 给出20个服从均值为0、标准差为3的正态分布N(0,32)随机数组成的表 Mathematica 命令 In[1]:= <<Statistics`ContinuousDistributions` In[2]:= rv =NormalDistribution[0,3 ]; RandomArray[rv, 20] Out[3] = {0.636589,-4.25557,2.04924,1.58478,0.0244065,0.371864,-0.933664, 3.54688,-0.888601,-0.650029,-2.49356,-3.07764,-2.44536,-0.512286, -1.68181,3.8912,-4.28302,-2.01939,-0.294215,2.13797}
实验 8 n个人每人携带一件礼物参加联欢会。联欢会开始后,先把所有的礼物编号,然后每人任意抽取一个号码,按号码领取礼物。请分别就参加联欢会的人数n=1到20人求所有人都得到别人赠送礼物的概率,并从这些概率值推断随着参加联欢会的人数增加是否会出现所有人都得到别人赠送礼物的概率会不断变小的情况? Mathematica 命令 In[1]:= p[n_]:=Sum[(-1)^k*1/k!,{k,2,n}] In[2]:= Table[N[p[k],18],{k,1,20}] Out[2]= {0,0.500000000000000000, 0.333333333333333333, 0.375000000000000000, 0.366666666666666667, 0.368055555555555556, 0.367857142857142857, 0.367881944444444444, 0.367879188712522046, 0.367879464285714286, 0.367879439233605900, 0.367879441321281599, 0.367879441160691161, 0.367879441172161906, 0.367879441171397190, 0.367879441171444985, 0.367879441171442173, 0.367879441171442329, 0.367879441171442321, 0.367879441171442322} 从计算结果可以看到,随着参会人数的增加,所有人都得到别人赠送礼物的概率不会不断变小,而是会收敛到一个约为0.367879,也就是e-1。
实验 9 在某纺织厂中,一个工人要照顾800个纱锭。每个纱锭旋转时,由于偶然的原因,纱会被扯断。假设在某一段时间内,每个纱锭的纱被扯断的概率为0.005,求在这段时间内,纱被扯断次数不大于10的概率。 分析:相当于进行800次独立试验,用X表示纱被扯断次数,则有X服从b(800,0.005)的二项分布,而所求概率为P{X≤10}可以用求b(800,0.005)的分布函数得到。 Mathematica 命令 In[1]:= <<Statistics`DiscreteDistributions` In[2]:= rvb=BinomialDistribution[800, 0.005]; f[k_]:=CDF[rvb, k] f[10] Out[4]= 0.997239 所以在这段时间内,纱被扯断次数不大于10的概率为0.997 239。
练习 1 某种检验方法对癌症的准确率时95%,一个人接受了检测并且结果呈阳性,假定这个人来自一个有100 000人口的地区,该地区2 000人得到这种癌症,推断接受检测者患这种癌症的概率是多少? 2 福利彩票摇奖的大转盘上的圆周等分成100份,第i份对应奖金i千元(i=1~100)。转动一次大转盘,求奖金金额为5万元到6万元的概率。 3 某厂产品中有4%废品,而在100件合格品中有80件一等品,求任取一件是一等品的概率。 4 生成服从二项分布b(25,0.3)的随机变量,显示其分布率,并画出其分布率图形。
实验 10 设样本数据为{ 110.1,25.2,50.5,50.5,55.7,30.2,35.4,30.2,4.9,32.3,50.5,30.5,32.3,74.2,60.8} 求该样本的均值、方差、标准差、中位数、众数。 Mathematica 命令 In[1]:= <<Statistics`DescriptiveStatistics` In[2]:= d1={ 110.1,25.2,50.5,50.5,55.7,30.2,35.4,30.2,4.9,32.3,50.5,30.5,32.3,74.2,60.8}; In[3]:= Mean[d1] Out[3]= 44.8867 (*均值为44.8867*) In[4]:= var=Variance[d1] Out[4]=614.89 (*方差为614.89*) In[5]:= Sqrt[var] Out[5]=24.797 (*标准差为24.797*) In[6]:=Median[d1] Out[6]:=35.4 (*中位数为 35.4*) In[7]:= Mode[d1] Out[7]:= 50.5 (*众数为50.5*)
实验 11 设样本数据为{16.5,13.8, 16.6, 15.7, 16.0, 16.4, 15.3},求该样本的均值、几何均值和调和均值。 Mathematica 命令 In[1]:= <<Statistics`DescriptiveStatistics` In[2]:= d2={16.5,13.8, 16.6, 15.7, 16.0, 16.4, 15.3}; In[3]:= Mean[d2] Out[3]= 15.7571 (*均值为15.7571*) In[4]:= GeometricMean[d2] Out[4]= 15.7296 (*几何均值为15.7296*) In[5]:= HarmonicMean [d2] Out[5]=15.7007 (*调和均值为15.7007*)
实验 12 设样本数据为{6.5, 3.8, 6.6, 5.7, 6.0, 6.4, 5.3} ,画出该样本的条形图和饼形图。 Mathematica 命令 In[1]:= << Graphics`Graphics` In[2]:= d3={6.5, 3.8, 6.6, 5.7, 6.0, 6.4, 5.3}; In[3]:= BarChart[d3] (*画样本条形图10-1*) In[4]:= PieChat[d3] (*画样本饼形图10-2*)
实验 13 某企业在1995到2001年间,每年的生产总值分别是上一年的60%、80%、90%、100%、105%、110%,试计算该企业生产总值的年平均发展速度。 Mathematica 命令 In[1]:= <<Statistics`DescriptiveStatistics` In[2]:= d6={ 0.60,0.80,0.90,1.00,1.05,1.10}; In[3]:= GeometricMean[d6] Out[3]= 0.89059
实验 14 设有如下30个样本数据 0.192, -1.382, 0.508, -0.813, 0.531, -0.536, 0.826, 1.404, -1.372, -0.349, 1.054, 1.372, 1.624, 0.709, 1.034, 1.670, -0.205, -0.017, -0.204, 0.056, -1.179,-0.645,1.201,0.453,0.304,-1.832,0.058,1.870,0.912,-1.769 (1)画出具有10个等距子区间的直方图; (2)画出具有16个等距子区间的直方图。 (1)Mathematica 命令 In[1]:= d ={0.192, -1.382, 0.508, -0.813, 0.531, -0.536, 0.826, 1.404, -1.372, -0.349, 1.054, 1.372, 1.624, 0.709, 1.034, 1.670, -0.205, -0.017, -0.204, 0.056, -1.179,-0.645,1.201,0.453,0.304,-1.832,0.058,1.870,0.912,-1.769 } In[2]:= h1 = -2; (*左端点值为-2*) h = ( 2 - (-2) )/10; (*小区间的长度为0.4*) num=10; (*小区间的个数为10*) dnum = 30; (*样本容量为30*) fpn=Table[m=0; Do[If[d[[i]]>=h1&&d[[i]]<h1+h,m=m+1],{i,1,dnum}]; h1=h1+h;m, {k,1,num}] Out[2]={2,2,2,2,4,4,4,4,3,3} In[3]:= << Graphics`Graphics` In[4]:= h1=-2; GeneralizedBarChart[ Table[{h1+h*i , fpn[[i+1]],h}, {i,0,num-1}], PlotRange->{0,5}]
(2)Mathematica 命令 In[1]:= d={0.192, -1.382, 0.508, -0.813, 0.531, -0.536, 0.826, 1.404, -1.372, -0.349, 1.054, 1.372, 1.624, 0.709, 1.034, 1.670, -0.205, -0.017, -0.204, 0.056, -1.179,-0.645,1.201,0.453,0.304,-1.832,0.058,1.870,0.912,-1.769} In[2]:= h1 = -2; (*左端点值为-2*) h = ( 2 - (-2) ) /16; (*小区间的长度为0.25*) num=16; (*小区间的个数为16*) dnum = 30; (*样本容量为30*) fpn=Table[m=0; Do[If[d[[i]]>=h1&&d[[i]]<h1+h,m=m+1],{i,1,dnum}]; In[3]:= h1=-2; GeneralizedBarChart[ Table[{h1+I*h, fpn[[i+1]],h}, {i,0,num-1}], PlotRange->{0,4}]
实验 15 某厂用自动包装机进行产品装包作业,现从某日生产的产品中随机抽取80包,测得数据如下 101.1, 100.6, 101.1, 101.7, 102.4, 102.7, 103.2, 103.7, 99.6, 99.1, 98.6, 98.1, 97.6, 96.8, 97.7, 98.2, 98.4,103.1,102.8,102.0,102.5,102.3,101.9,101.2, 101.1, 99.6, 99.9, 99.1, 98.1, 102.2, 102.3, 101.8, 101.7,102.0,101.8, 101.8,102.0, 101.5,101.3,101.4, 100.9,100.6,98.6,100.2,100.8,101.4,101.5,101.3,99.4, 99.5, 99.1, 101.0, 100.3, 100.5, 100.0,99.9, 99.7,99.6, 100.4, 100.3, 100.2, 100.0, 100.1, 100.5,99.8, 99.6, 100.0, 100.3,100.5,100.2, 99.0, 98.6, 99.4, 99.3, 99.1, 100.1, 100.2, 101.4, 100.9, 101.0 求该自动包装机装包重量的均值、方差、中位数、众数,并画出该样本数据的直方图。 Mathematica 命令 In[1]:= d={101.1, 100.6, 101.1, 101.7, 102.4, 102.7, 103.2, 103.7, 99.6, 99.1, 98.6, 98.1, 97.6, 96.8, 97.7, 98.2, 98.4,103.1,102.8,102.0,102.5,102.3,101.9,101.2, 101.1, 99.6, 99.9, 99.1, 98.1, 102.2, 102.3, 101.8, 101.7,102.0,101.8, 101.8,102.0, 101.5,101.3,101.4, 100.9,100.6,98.6,100.2,100.8,101.4,101.5,101.3,99.4, 99.5, 99.1, 101.0, 100.3, 100.5, 100.0,99.9, 99.7,99.6, 100.4, 100.3, 100.2, 100.0, 100.1, 100.5,99.8, 99.6, 100.0, 100.3,100.5,100.2, 99.0, 98.6, 99.4, 99.3, 99.1, 100.1, 100.2, 101.4, 100.9, 101.0} ;
In[2]:= Mean[d] Out[2]= 100.49 In[3]:= Variance[d] Out[3]=2.02344 In[4]:=Median[d] Out[4]:= 100.45 In[5]:= Mode[d] Out[5]:={99.1,99.6,100.2} In[6]:=<<Graphics`Graphics` In[7]:= Histogram[d]
实验 16 某矿脉中10个相邻样本点处一种伴生金属的含量数据如下表所示。 1)分析 本题没有事先指定回归关系,因此应改通过散点图形状自己尝试,从中找出较为合适的回归形式。 2)实验操作 In[1]:=dat = {{2,106.42},{3,108.2},{4,109.58},{5,109.5},{7,110.},{8,109.93}, {10,110.49},{11,110.59},{14,110.6},{15,110.9}} In[2]:=ListPlot[d] Out[2]=
粗劣观察此散点图,觉得似乎回归方程形式可取线性形式:y= a+bx. 但仔细观察,发现散点图有一些上凸,而上凸的图形具有y= a+b形式。下面用Mathematica 数学软件尝试这两种形式回归的优劣。 In[3]:=<<Statistics`LinearRegression` In[4]:= r=Regress[dat,{1,x},x] Out[4]= In[5]:= r = Regress[dat,{1,Sqrt[x]},x] Out[5]= 实验结果发现,两个回归方程的线性关系都是高度显著的,但回归方程y= a+b无论在变量系数b的检验概率还是对应回归方程检验概率都比回归方程y= a+bx要小近一个数量级,因此回归方程y= a+b要比y= a+bx形式好。我们最后选定本题的回归方程为: y= 105.767+1.43079
练习 1 设样本数据为{0.2,0.1,0.3,0.4},画出该样本的条形图和饼形图。 样本{1050,1100,1120,1250,1280}来自正态总体,求该总体均值、方差的置信度为0.95置信区间。 如果一组样本数据的散点图具有y=aebx形状,怎样求其回归方程?
第9章结束 • 谢谢