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

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

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

3.1 随机变量分布与数字特征实例及MA TLAB求解 (1)

3.1.1 MATLAB实现 (1)

3.1.2 相关实例求解 (2)

3.2 数理统计实例分析及MATLAB求解 (4)

3.1.1 MATLAB实现 (4)

3.1.2 相关实例求解 (4)

3.3参数估计与假设检验实例分析及MATLAB求解 (5)

3.1.1 MATLAB实现 (5)

3.1.2 相关实例求解 (5)

3.4 方差分析实例求解 (10)

3.1.1 MATLAB实现 (10)

3.1.2 相关实例求解 (10)

3.5 判别分析应用实例及求解 (14)

3.1.1 MATLAB实现 (14)

3.1.2 相关实例求解 (14)

3.6 聚类分析应用实例及MATLAB求解 (16)

3.1.1 MATLAB实现 (16)

3.1.2 相关实例求解 (16)

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) % 生成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,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 设X 的概率密度为????

?????

<<-≤≤=其他。0;30001500,1500

3000;15000,

1500)(2

2x x

x x x f ，求)(X E 。

求解程序：

syms x

f1=x/1500^2;

f2=(3000-x)/1500^2;

Ex=int(x*f1,0,1500)+int(x*f2,1500,3000)

运行结果：

Ex =1500

Ex =1/3

例2-3：绘制 =0.5,1,3,5,10 时Poisson 分布的概率密度函数与概率分布函数曲线。

代码如下：

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))];

end

plot(x,y1), figure; plot(x,y2)

图2 泊松分布概率密度函数图

图3 泊松分布概率分布函数

3.2 数理统计实例分析及MATLAB 求解

3.1.1 MATLAB 实现

在MATLAB 中各种随机数可以认为是独立同分布的，即简单随机样本。

常用分布的随机数产生方法，可用分布英文名称缩写加上rnd ，例如：

x=betarnd(a,b,m,n) 参数为a,b 的beta 分布； x=binornd(N,p,m,n) 参数为N,p 的二项分布；

3.1.2 相关实例求解

例2-4：设总体密度函数

cos ,,2

2

2()0,.

x

x f x π

π

?-

<<

?

=???

其他 试从该总体中抽取容量为1000的简单随机样本。

解 利用MATLAB 编辑窗口保存以下程序，保存为ex11.m

n=1000; x=zeros(1,n); k=0; while k

a=rand*pi-pi/2; b=rand/2;

if b<(cos(a)/2)

k=k+1;

x(k)=a;

end

end

hist(x,-pi/2:0.2:pi/2)

保存完成之后，在命令窗口执行ex11，则x被赋值，且可以得到这个容量为1000的样本的直方图。

图7 直方图

3.3参数估计与假设检验实例分析及MATLAB求解

3.1.1 MATLAB实现

3.1.2 相关实例求解

例3-5:对某型号的20辆汽车记录其5L汽油的行驶里程（公里），观测数据如下：

29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 28.7

28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.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];

x=[x1 x2]';

muhat=mean(x)

sigma2hat=moment(x,2)%样本二阶中心矩

var(x,1);

运行结果：

muhat = 28.6950

sigma2hat =0.9185

例3-6：对某型号的20辆汽车记录其5L汽油的行驶里程（公里），观测数据如下：

29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 28.7

28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.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];

x=[x1 x2]';

p=mle('norm',x);

muhatmle=p(1)

sigma2hatmle=p(2)^2

运行结果：

muhatmle =28.6950

sigma2hatmle =0.9185

例3-7设两台车床加工同一零件，各加工8件，长度的误差为：

A：-0.12 -0.80 -0.05 -0.04 -0.01 0.05 0.07 0.21 B：-1.50 -0.80 -0.40 -0.10 0.20 0.61 0.82 1.24 求方差比的置信区间。

解：用Matlab计算如下：

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.9 10.9 11.1 9.6 10.2 10.3 9.6 9.9 11.2 10.6 9.8 10.5 10.1 10.5 9.7

设装配时间的总体服从正态分布，标准差为0.4，是否可以认为装配时间的均值在0.05的水平下不小于10。

解：

0H ：10<μ vs 1H ：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) 运行结果： h =1

sig =0.01267365933873 muci =

10.05287981908398 Inf

因此，在0.05的水平下，可以认为装配时间的均值不小于10。

例3-9 某种电子元件的寿命X (以小时计)服从正态分布，μ和2σ均未知。现测得16只元件的寿命如下：

159 280 101 212 224 379 179 264 222 362 168 250 149 260 485 170 问是否有理由认为元件的平均寿命大于225(小时)？ 解：

0H ：225≤μ vs 1H ：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=0 sig=0.2570

198.2321 Inf

由于sig=0.257，因此没有充分的理由认为元件的平均寿命大于225小时。 而对于

0H ：225≥μ vs 1H ：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=0 sig=0.7430 muci =

-Inf 284.7679

由于sig=0.743，因此更没有充分的理由认为元件的平均寿命小于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

根据专业经验，硬度服从正态分布，且方差保持不变，试在显著性水平05.0=α下判断镍合金的硬度是否有明显提高。

解：

0H ：铜镍μμ≤ vs 1H ：铜镍μμ>

程序：

%正态总体的方差齐但求知时的均值检验

x=[76.43 76.21 73.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=1 sig=0.0142

1.4148 Inf

因此，在显著性水平05.0=α下，可以判断镍合金的硬度有明显提高。 例3-11 q-q 图 程序：

%分布拟合检验(q-q 图法)

x=normrnd(0,1,100,1);%生成服从正态分布的随机数 y=normrnd(0.5,2,50,1);

z=weibrnd(2,0.5,100,1);%生成服从威布尔分布的随机数 subplot(2,2,1)%说明生成子图的位置 qqplot(x);hold on

subplot(2,2,2);qqplot(x,y);hold on subplot(2,2,3);qqplot(z);hold on subplot(2,2,4);qqplot(x,z); hold off

运行结果见下图：

图12 q-q 分布图

结论：上面的q-q 图中第1个子图用x 的数据绘图，因为服从正态分布，图中数据点呈直线分布；第2个子图用x 数据和y 数据（均服从正态分布），数据点的主体部分呈直线；第3个子图用z 数据绘图，由于它服从威布尔分布，所以数据点不在一条直线上；第4个子图是用x 数据和z 数据绘制的，因为它们不是

同分布的，图中数据点不呈直线分布。

例3-12：诸X 为 1.1，3.3，5.5，7.7，诸Y 为2.2，4.4，6.6，以下列表给出混合样本及秩

混合样本

1.1

3.3

5.5

7.7

2.2

4.4

6.6

秩 1 3 5 7 2 4 6

则

24612W =++=.

若H0成立，则W 的值应该适中。注意到每个秩序的平均值为(1)2n +，故H0成立时，E(W)=n(N+1)2，W 的值在此值附近应该是正常的。若W 的值异常偏大，说明第二个总体确有增加效应。利用MATLAB 自身的函数

p = ranksum(X,Y)

可以进行双侧的秩和检验。返回的p 值小于给定的α则拒绝原假设，认为H1：

0?≠成立。

H0成立时，可以证明W 关于E(W)=n(N+1)2对称，要检验H1：0?>，只要判定(1)2n N W +>，并且p = ranksum(X,Y)2α<即可。

3.4 方差分析实例求解

3.1.1 MATLAB 实现 3.1.2 相关实例求解

例3-13 某化工产品的产量是衡量经济效益的重要指标。为了考察反应温度（因子A ）对该化工产品产量y 是否有显著影响，我们选取因子A 的5个水平为

1A ：60℃，2A ：65℃，3A ：70℃，4A ：75℃，5A ：80℃。每个水平下各作5次重复试验，试验结果如下表所示：

表8 反应温度观测值

设该试验的线性统计模型为：

2

,

1,,5;5...(0,).ij i ij ij

y i j i i d N μεεσ=+?==?? 诸 程序：

%等重复的单因子试验的方差分析

y1=[90 87 92 91 88]';y2=[97 91 93 95 92]'; 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 水箱平衡图

结论：

可见，反应温度（因子A）对该化工产品产量y有显著影响。

例3-14 为了考察某种电池的最大输出电压受板极材料与使用电池的环境温度的影响，材料类型(因子A)取3个水平(即3种不同的材料)，温度(因子B)也取3个水平，每个水平组合下重复4次试验，数据如下：

表 9 实验数据

解：程序

1).%两因子等重复试验（方差分析）

y1=[130 155 174 180 34 40 80 75 20 70 82 58]';

y2=[150 188 159 126 136 122 106 115 25 70 58 45]';

y3=[138 110 168 160 174 120 150 139 96 104 82 60]';

y=[y1 y2 y3];

p=anova2(y,4)

运行结果：

图15 运行结果图

结论：列因子A、行因子B以及其交互作用均显著。

2).再作列因子的多重比较，程序如下：

%两因子等重复试验（列因子的多重比较）

y1=[130 155 174 180 34 40 80 75 20 70 82 58]';

y2=[150 188 159 126 136 122 106 115 25 70 58 45]';

y3=[138 110 168 160 174 120 150 139 96 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 significant at the level you specified. Testing main effects under these conditions is questionable.

图16 分布示意图

可见，列因子A、行因子B以及其交互作用均显著。

3.5 判别分析应用实例及求解

3.1.1 MATLAB实现

3.1.2 相关实例求解

例3-15我国山区某大型化工厂，在厂区及邻近地区挑选有代表性的15个大气取样点，每日4次同时抽取大气样品，测定其中含有的6种气体的浓度，前后共4天，每个取样点每种气体实测16次。计算每个取样点每种气体的平均浓度，数据如下表所示，气体数据对应的污染地区分类如下表中最后一列所示。现有两个取自该地区的4个气体样本，气体指标如表中后4行所示，试判别这4个样品的污染分类。

表10 指标分布表

解:首先将数据保存在文件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','H2:H16');%样本的分组信息数据

obs=[1:19]';样本编号

%距离判别，判别函数类型为mahalanobis，返回判别结果向量C和判误概率err [C,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 2

err =

由此可知，气体样品1和样品2污染分类属于第一类，而气体样品3和样品4则属于污染分类的第二类，且判误概率极小，约等于0，这说明该分类结果是可行的。

3.6 聚类分析应用实例及MATLAB求解

3.1.1 MATLAB实现

3.1.2 相关实例求解

例3-16为了研究世界各国森林、草原资源的分布规律，共抽取了21个国家的数据，每个国家4项指标，原始数据见下表。试用该数据对国别进行聚类分析。

表 16 世界绿化情况表

MATLAB提供了两种方法进行聚类分析。

一种是一次聚类。它的优点是可利用clusterdata函数对样本数据进行一次聚类。其缺点是可供用户选择的面比较窄，不能更改距离的计算方法。

另一种是分步聚类。可以分以下步骤进行分步聚类：

⑴找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变

量之间的距离；

⑵用linkage函数定义变量之间的连接；

⑶用Cophenetic函数评价聚类信息；

⑷cluster函数创建聚类。

1、一次聚类

X=xlsread('d.xls','B2:E22');

T=clusterdata(X,0.9)

通过分析T =[10，10，3，8，2，9，8，10，10，7，1，8，1，8，8，1，4，9，9，5，6]，可知最后21个国家的绿化情况可分为5类：{中国，美国，法国，澳大利亚，印度尼西亚，尼日尼亚}，{德国，意大利，波兰，南斯拉斯，罗马尼亚}，{墨西哥，巴西}，{印度，日本}，{英国，捷克，匈牙利，保加利亚}。

2、分步聚类

⑴．寻找变量之间的相似性

用pdist函数计算相似矩阵，有多种方法可以计算距离，进行计算之前最好先将数据用zscore函数进行标准化。

X2=zscore(X); %标准化数据

Y2=pdist(X2); %计算距离

⑵．定义变量之间的连接

Z2=linkage(Y2);

⑶．评价聚类信息

C2=cophenet(Z2,Y2);

运行结果C2=0.9470，可见，其聚类效果较佳。

⑷．创建聚类，并作出谱系图

T=cluster(Z2,6);

H=dendrogram(Z2);

分类结果：{中国，美国，法国，澳大利亚，印度尼西亚，尼日尼亚}，{德国，意大利，波兰，南斯拉斯，罗马尼亚}，{墨西哥，巴西}，{印度，日本}，{英国，捷克，匈牙利，保加利亚}。

4

71512 414111613 5 61819 31720 1 9 2 82110

图17 聚类图