第3章 概率统计实例分析及MatlAb求解

1、第3章 概率统计实例分析及MatlAb求解3.1 随机变量分布与数字特征实例及MATLAB求解3.1.1 MATLAB实现用mvnpdf和mvncdf函数可以计算二维正态分布随机变量在指定位置处的概率和累积分布函数值。利用MATLAB统计工具箱提供函数，可以比较方便地计算随机变量的分布律（概率密度函数）、分布函数及其逆累加分布函数，见附录2-1，2-2，2-3。MATLAB中矩阵元素求期望和方差的函数分别为mean和var，若要求整个矩阵所有元素的均方差，则要使用std2函数。随机数生成函数：rand( )和randn( )两个函数伪随机数生成函数： A=gamrnd(a,lambda,n,m

2、) % 生成n*m的分布的伪随机矩阵 B=raylrnd(b,n,m) %生成rayleigh的伪随机数3.1.2 相关实例求解例2-1计算服从二维正态分布的随机变量在指定范围内的累积分布函数值并绘图。程序：%二维正态分布的随机变量在指定范围内的累积分布函数图形mu=0 0;sigma=0.25 0.3;0.3 1;%协方差阵x=-3:0.1:3;y=-3:0.2:3;x1,y1=meshgrid(x,y);%将平面区域网格化取值f=mvncdf(x1(:) y1(:),mu,sigma);%计算累积分布函数值F=reshape(f,numel(y),numel(x);%矩阵重塑surf(x,

3、y,F);caxis(min(F(:)-0.5*range(F(:),max(F(:);%range(x)表示最大值与最小值的差，即极差。axis(-3 3 -3 3 0 0.5);xlabel(x);ylabel(y);zlabel(Probability Density);图1 二维正太分布累积分布函数值图例2-2设的概率密度为，求。求解程序：syms xf1=x/15002;f2=(3000-x)/15002;Ex=int(x*f1,0,1500)+int(x*f2,1500,3000)运行结果：Ex =1500Ex =1/3例2-3：绘制=0.5,1,3,5,10 时 Poisson

4、分布的概率密度函数与概率分布函数曲线。代码如下：x=0:15; y1=; y2=; lam1=0.5,1,3,5,10;for i=1:length(lam1)y1=y1,poisspdf(x,lam1(i); y2=y2,poisscdf(x,lam1(i);endplot(x,y1), figure; plot(x,y2)图2 泊松分布概率密度函数图图3 泊松分布概率分布函数3.2 数理统计实例分析及MATLAB求解3.1.1 MATLAB实现在MATLAB中各种随机数可以认为是独立同分布的，即简单随机样本。常用分布的随机数产生方法，可用分布英文名称缩写加上rnd，例如：x=betarnd

5、(a,b,m,n) 参数为a,b的beta分布；x=binornd(N,p,m,n) 参数为N,p的二项分布；3.1.2 相关实例求解例2-4：设总体密度函数试从该总体中抽取容量为1000的简单随机样本。解 利用MATLAB编辑窗口保存以下程序，保存为ex11.mn=1000;x=zeros(1,n);k=0;while kn a=rand*pi-pi/2; b=rand/2; if b(cos(a)/2) k=k+1; x(k)=a; endendhist(x,-pi/2:0.2:pi/2)保存完成之后，在命令窗口执行ex11，则x被赋值，且可以得到这个容量为1000的样本的直方图。图7 直

6、方图3.3参数估计与假设检验实例分析及MATLAB求解3.1.1 MATLAB实现3.1.2 相关实例求解例3-5:对某型号的20辆汽车记录其5L汽油的行驶里程（公里），观测数据如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.9试估计总体的均值和方差。求解程序：%矩法估计x1=29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 28.7;x2=28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.9;

7、x=x1 x2;muhat=mean(x)sigma2hat=moment(x,2)%样本二阶中心矩var(x,1);运行结果：muhat = 28.6950sigma2hat =0.9185例3-6：对某型号的20辆汽车记录其5L汽油的行驶里程（公里），观测数据如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.9设行驶里程服从正态分布，试用最大似然估计法估计总体的均值和方差。求解程序：%最大似然估计x1=29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.

8、0 27.9 28.7;x2=28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.9;x=x1 x2;p=mle(norm,x);muhatmle=p(1)sigma2hatmle=p(2)2运行结果：muhatmle =28.6950sigma2hatmle =0.9185例3-7 设两台车床加工同一零件，各加工8件，长度的误差为：A：-0.12 -0.80 -0.05 -0.04 -0.01 0.05 0.07 0.21B：-1.50 -0.80 -0.40 -0.10 0.20 0.61 0.82 1.24求方差比的置信区间。解：用Matlab

9、计算如下：x=-0.12,-0.80,-0.05,-0.04,-0.01,0.05,0.07,0.21;y=-1.50,-0.80,-0.40,-0.10,0.20,0.61, 0.82,1.24;v1=var(x); v2=var(y);c1=finv(0.025,7,7); c2=finv(0.975,7,7);a=(v1/v2)/c2; b=(v1/v2)/c1; a,b计算结果为： 0.0229 0.5720结论：方差比小于1的概率至少达到了95%，说明车床A的精度明显高。例3-8 下面列出的是某工厂随机选取的20只零部件的装配时间（分）：9.8 10.4 10.6 9.6 9.7 9

10、.9 10.9 11.1 9.6 10.210.3 9.6 9.9 11.2 10.6 9.8 10.5 10.1 10.5 9.7设装配时间的总体服从正态分布，标准差为0.4，是否可以认为装配时间的均值在0.05的水平下不小于10。解：程序：%正态总体的方差已知时的均值检验x1=9.8 10.4 10.6 9.6 9.7 9.9 10.9 11.1 9.6 10.2;x2=10.3 9.6 9.9 11.2 10.6 9.8 10.5 10.1 10.5 9.7;x=x1 x2;m=10;sigma=0.4;a=0.05;h,sig,muci=ztest(x,m,sigma,a,1)运行结果

11、：h =1sig =0.01267365933873muci = 10.05287981908398 Inf因此，在0.05的水平下，可以认为装配时间的均值不小于10。例3-9某种电子元件的寿命(以小时计)服从正态分布，和均未知。现测得16只元件的寿命如下：159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170问是否有理由认为元件的平均寿命大于225(小时)？解：程序：%正态总体的方差未知时的均值检验x=159 280 101 212 224 379 179 264 222 362 168 250 149 260 485

12、170;m=225;a=0.05;h,sig,muci=ttest(x,m,a,1)运行结果：h=0sig=0.2570muci= 198.2321 Inf由于sig=0.257，因此没有充分的理由认为元件的平均寿命大于225小时。而对于：程序：%正态总体的方差未知时的均值检验x=159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170;m=225;a=0.05;h,sig,muci=ttest(x,m,a,-1)运行结果：h=0sig=0.7430muci = -Inf 284.7679由于sig=0.743，因此更没有充

13、分的理由认为元件的平均寿命小于225小时。例3-10某厂铸造车间为提高铸件的耐磨性而试制了一种镍合金铸件以取代铜合金铸件，为此，从两种铸件中各独立地抽取一个容量分别为8和9的样本，测得其硬度（一种耐磨性指标）为：镍合金 76.43 76.21 73.58 69.69 65.29 70.83 82.75 72.34铜合金 73.66 64.27 69.34 71.37 69.77 68.12 67.27 68.07 62.61根据专业经验，硬度服从正态分布，且方差保持不变，试在显著性水平下判断镍合金的硬度是否有明显提高。解：程序：%正态总体的方差齐但求知时的均值检验x=76.43 76.21 7

14、3.58 69.69 65.29 70.83 82.75 72.34;y=73.66 64.27 69.34 71.37 69.77 68.12 67.27 68.07 62.61;a=0.05;h,sig,ci=ttest2(x,y,a,1)运行结果：h=1sig=0.0142ci = 1.4148 Inf因此，在显著性水平下，可以判断镍合金的硬度有明显提高。例3-11q-q图程序：%分布拟合检验(q-q图法)x=normrnd(0,1,100,1);%生成服从正态分布的随机数y=normrnd(0.5,2,50,1);z=weibrnd(2,0.5,100,1);%生成服从威布尔分布的随机

15、数subplot(2,2,1)%说明生成子图的位置qqplot(x);hold onsubplot(2,2,2);qqplot(x,y);hold onsubplot(2,2,3);qqplot(z);hold onsubplot(2,2,4);qqplot(x,z);hold off运行结果见下图：图12 q-q分布图结论：上面的q-q图中第1个子图用x的数据绘图，因为服从正态分布，图中数据点呈直线分布；第2个子图用x数据和y数据（均服从正态分布），数据点的主体部分呈直线；第3个子图用z数据绘图，由于它服从威布尔分布，所以数据点不在一条直线上；第4个子图是用x数据和z数据绘制的，因为它们不是

16、同分布的，图中数据点不呈直线分布。例3-12：诸为 1.1，3.3，5.5，7.7，诸为2.2，4.4，6.6，以下列表给出混合样本及秩混合样本1.13.35.57.72.24.46.6秩1357246则.若H0成立，则的值应该适中。注意到每个秩序的平均值为，故H0成立时，的值在此值附近应该是正常的。若的值异常偏大，说明第二个总体确有增加效应。利用MATLAB自身的函数p = ranksum(X,Y)可以进行双侧的秩和检验。返回的p值小于给定的则拒绝原假设，认为H1：成立。H0成立时，可以证明关于对称，要检验H1：，只要判定，并且p = ranksum(X,Y)即可。3.4 方差分析实例求解3

17、.1.1 MATLAB实现3.1.2 相关实例求解例3-13某化工产品的产量是衡量经济效益的重要指标。为了考察反应温度（因子）对该化工产品产量是否有显著影响，我们选取因子的5个水平为：60，：65，：70，：75，：80。每个水平下各作5次重复试验，试验结果如下表所示：表8 反应温度观测值反映温度观察值（）606570758090 87 92 91 8897 91 93 95 9296 92 93 96 9584 82 86 83 8884 81 85 86 82设该试验的线性统计模型为： 程序：%等重复的单因子试验的方差分析y1=90 87 92 91 88;y2=97 91 93 95 9

18、2;y3=96 92 93 96 95;y4=84 82 86 83 88;y5=84 81 85 86 82;y=y1 y2 y3 y4 y5;p=anova1(y)%y的列表示重复观测值。运行结果：图13 运行结果图14 水箱平衡图结论：可见，反应温度（因子）对该化工产品产量有显著影响。例3-14为了考察某种电池的最大输出电压受板极材料与使用电池的环境温度的影响，材料类型(因子)取3个水平(即3种不同的材料)，温度(因子)也取3个水平，每个水平组合下重复4次试验，数据如下：表 9 实验数据温度152535材料类型1130 155 174 18034 40 80 7520 70 82 582

19、150 188 159 126136 122 106 11525 70 58 453138 110 168 160174 120 150 13996 104 82 60 解：程序1).%两因子等重复试验（方差分析）y1=130 155 174 18034 40 80 7520 70 82 58;y2=150 188 159 126136 122 106 11525 70 58 45;y3=138 110 168 160174 120 150 13996 104 82 60;y=y1 y2 y3;p=anova2(y,4)运行结果：图15 运行结果图结论：列因子A、行因子B以及其交互作用均显著。

20、2).再作列因子的多重比较，程序如下：%两因子等重复试验（列因子的多重比较）y1=130 155 174 18034 40 80 7520 70 82 58;y2=150 188 159 126136 122 106 11525 70 58 45;y3=138 110 168 160174 120 150 13996 104 82 60;y=y1 y2 y3;p,table,stats=anova2(y,4)multcompare(stats,0.05)运行结果：Note: Your model includes an interaction term that is significanta

21、t the level you specified. Testing main effects under theseconditions is questionable.图16 分布示意图可见，列因子A、行因子B以及其交互作用均显著。3.5 判别分析应用实例及求解3.1.1 MATLAB实现3.1.2 相关实例求解例3-15我国山区某大型化工厂，在厂区及邻近地区挑选有代表性的15个大气取样点，每日4次同时抽取大气样品，测定其中含有的6种气体的浓度，前后共4天，每个取样点每种气体实测16次。计算每个取样点每种气体的平均浓度，数据如下表所示，气体数据对应的污染地区分类如下表中最后一列所示。现有两

22、个取自该地区的4个气体样本，气体指标如表中后4行所示，试判别这4个样品的污染分类。表10 指标分布表气体氯硫化氢二氧化碳碳4环氧氯丙烷环己烷污染分类10.0560.0840.0310.0380.00810.022120.040.0550.10.110.0220.007130.050.0740.0410.0480.00710.02140.0450.050.110.10.0250.006150.0380.130.0790.170.0580.043260.030.110.070.160.050.046270.0340.0950.0580.160.20.029180.030.090.0680.180.

23、220.039190.0840.0660.0290.320.0120.0412100.0850.0760.0190.30.010.042110.0640.0720.020.250.0280.0382120.0540.0650.0220.280.0210.042130.0480.0890.0620.260.0380.0362140.0450.0920.0720.20.0350.0322150.0690.0870.0270.050.0890.0211样品10.0520.0840.0210.0370.00710.022样品20.0410.0550.110.110.0210.007样品30.030.1

24、120.0720.160.0560.021样品40.0740.0830.1050.190.021解:首先将数据保存在文件example12_4.mat中，然后再编写程序。注意：对于样本数据的每一种分类，要求其对应的观测案例数（即training的行数）必须大于观测的变量个数（即training的列数）。例如本例中，变量个数为6，则对应于第1类和第2类的观测案例数必须大于6，否则将出错。代码如下：sample=xlsread(b.xls,B2:G20);% 读取文件全部数据training=xlsread(b.xls,B2:G16);% 已知组别的样本数据group=xlsread(b.xls,

25、H2:H16);%样本的分组信息数据obs=1:19;样本编号%距离判别，判别函数类型为mahalanobis，返回判别结果向量C和判误概率errC,err = classify(sample,training,group,mahalanobis);obs, C % 查看判别结果err % 查看误判概率运行结果：ans =1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 1 1 1 1 2 2 1 1 2 2 2 2 2 2 1 1 1 2 2err = 0由此可知，气体样品1和样品2污染分类属于第一类，而气体样品3和样品4则属于污染分类的第二类，

26、且判误概率极小，约等于0，这说明该分类结果是可行的。 3.6 聚类分析应用实例及MATLAB求解3.1.1 MATLAB实现3.1.2 相关实例求解例3-16为了研究世界各国森林、草原资源的分布规律，共抽取了21个国家的数据，每个国家4项指标，原始数据见下表。试用该数据对国别进行聚类分析。表 16 世界绿化情况表国别森林面积（万公顷）森林覆盖率（%）林木蓄积量（万公顷）草原面积（万公顷）中国1197812.593.531908美国2844630.420223754日本250167.224.858德国102828.414599英国2108.61.51147法国145826.7161288意大利63521.13.6514加拿大3261332.7192.82585澳大利亚1070013.910.545190前苏联9200041.1841.537370捷克45835.88.9168波兰86827.811.440