《 数学实验》7  概率论 样本描述 参数估计 假设检验  方差分析

7.1 概率论 利用MATLAB统计工具箱，可以进行基本概论和数理统计分析，以及进行比较复杂的多元统计分析。 7.1.1 分布率和概率密度函数（P132表7-1） 以正态分布为例，用normpdf函数计算其概率密度函数，调用格式为： Y=normpdf(X,MU,SIGMA) 计算数据X中各值处参数为MU和SIGMA的正态概率密度函数的值。其中参数SIGMA必须为正。 2/25

7.1.2 分布函数 [例7-1] 计算参数为mu和1的正态分布概率密度函数在1.5处的值，其中mu为1到2之间以0.2为间隔的小数。 y=normpdf(1.5,mu,1) y=0.1295 0.1714 0.2179 0.2661 0.3123 0.3521 0.3814 0.3970 0.3970 0.3814 0.3521 7.1.2 分布函数 若X为随机变量，x为任意实数，则函数 为X的分布函数。如果知道X的分布函数，就可以知道落在任一区间(x1,x2)上的概率。 3/25

p=normcdf(X,MU,SIGMA) 用normcdf函数计算正态分布的分布函数，调用格式为： 计算参数为MU和SIGMA的正态分布分布函数在数据X中每个值处的值。其中参数SIGMA必须为正。 [例7-2]求标准正态分布的一个观察量落在区间[-1 1]中的值。 p=normcdf([-1,1]); p(2)-p(1) ans= 0.6872 4/25

7.1.3 随机变量的数字特征 返回相应分布的数学期望和方差 [M,V]=binostat(N,P) [M,V]=expstat(MU) 7.1.3 随机变量的数字特征 [M,V]=binostat(N,P) [M,V]=expstat(MU) [M,V]=normstat(MU,SIGMA) C=cov(X) 返回X的协方差或协方差矩阵 C=cov(X,Y) 返回X与Y的协方差矩阵 R=corrcoef(X) 返回源于矩阵的相关系数矩阵 M=moment(X,order) 返回X的order阶中心矩 [例7-3]生成一个6行5列的随机矩阵，然后计算每列数据的3阶中心矩。 返回每一列的3阶中心矩 X=randn([6,5]); M=moment(X,3) 5/25

7.2 样本描述 7.2.1 集中趋势 描述样本数据集中趋势的统计量有算术平均值、中位数、众数、几何均值、调和均值和截尾均值等。 7.2 样本描述 7.2.1 集中趋势 描述样本数据集中趋势的统计量有算术平均值、中位数、众数、几何均值、调和均值和截尾均值等。 样本数据x1,x2,…,xn的 （1）几何均值 m=geomean(X) （2）调和均值 m=harmmean(X) 6/25

对样本数据进行排序后，去掉两端的部分极值，然后对剩下的数据求算术平均值，得到截尾均值。 （3）算术平均值 m=mean(X) m=median(X) （4）中值 （5）截尾均值 对样本数据进行排序后，去掉两端的部分极值，然后对剩下的数据求算术平均值，得到截尾均值。 m=trimmean(X,percent) 剔除测量值中最大和最小percent%的数据后，计算样本X的均值。 7/25

描述离散趋势的统计量包括均值绝对差、极差、方差和标准差。 7.2.2 离中趋势 描述离散趋势的统计量包括均值绝对差、极差、方差和标准差。 （1）均值绝对差 y=mad(X) 若X为矢量，则y用mead(abs(X-mean(X)))计算；若X 为矩阵，则y为包含X中每列数据均值绝对差的行矢量。 mad(X,0): 与mad(X)相同，使用均值。 mad(X,1): 基于中值计算y,即:median(abs(X-median(X))). （2）极差 y=range(X) 返回X中数据的最小值与最大值之间的差值。 8/25

7.3 参数估计 （3）方差 y=var(X) （4）标准差 y=std(X) 7.3.1 点估计：用单个数值作为参数的估计 7.3 参数估计 7.3.1 点估计：用单个数值作为参数的估计 （1）矩法：用总体的样本矩来估计总体的同阶矩。 [例7-13]随机取8个活塞环，测得它们的直径为(以mm计):74.001 74.005 74.003 74.001 74.000 73.998 74.006 74.002,设环直径的测量值服从正态分布，现估计总体的方差。 解：因为样本的2阶中心矩是总体方差的矩估计量，所以可以用moment函数进行估计。 X=[74.001 74.005 74.003 74.001 74.000 73.998 74.006 74.002]； moment(X,2) 9/25

使用data矢量中的样本数据，返回dist指定的分布的最大似然估计。 （2）最大似然法 p=mle(‘dist’,data) 使用data矢量中的样本数据，返回dist指定的分布的最大似然估计。 [例7-14]用最大似然估计法解例7-3。 X=[74.001 74.005 74.003 74.001 74.000 73.998 74.006 74.002]； p=mle(‘norm’,X); p(2)*p(2) 7.3.2 区间估计： 区间估计不仅仅给出了参数的近似取值，还给出了取该值的误差范围。 10/25

[phat,pci]=mle(‘dist’,data):返回最大似然估计和95%置信区间。 [phat,pci]=mle(‘dist’,data,alpha):返回指定分布的最大似然估计值和100(1-alpha)置信区间。 [phat,pci]=mle(‘dist’,data,alpha,p1):该形式仅用于二项分布，其中p1为试验次数。 [例7-15]从一批灯泡中随机地取5只作寿命试验,测得寿命(以小时计)为:1050 1100 1120 1250 1280,设灯泡寿命服从正态分布,求灯泡寿命平均值的95%置信区间. X=[1050 1100 1120 1250 1280]; [p,ci]=mle(‘norm’,X,0.05) p=1.0e+003* ci=1.0e+003* 1.1600 0.0892 1.0818 0.0339 1.2382 0.1445 11/25

7.3.3 常见分布的参数估计(P142表7-4) Matlab统计工具箱还提供了具体函数的参数估计函数。如用normfit函数对正态分布总体进行参数估计: [muhat,sigmahat,muci,sigmaci]=normfit(X):对于给定的服从正态分布的数据矩阵X,返回参数 和 的估计值nuhat和sigmahat。muci和sigmaci为 和 的95%置信区间。 [muhat,sigmahat,muci,sigmaci]=normfit(X,alpha):进行参数估计并计算100(1-alpha)置信区间。 [例7-16]用normfit函数求解例7-6。 X=[1050 1100 1120 1250 1280]; [muhat,sigmahat,muci,sigmaci]=normfit(X,0.05) muhat=1160 sigmahat=99.7497 muci=1.0e+003* sigmaci=59.7633 1.0361 289.6364 1.2839 12/25

7.4 假设检验 7.4.1 单个正态总体均值的假设检验 检验关于分布或参数未知的总体的假设是否合理 7.4 假设检验 　检验关于分布或参数未知的总体的假设是否合理 7.4.1 单个正态总体均值的假设检验 在均方差未知时，用t统计量检验样本均值的显著性 函数：ttest；调用格式： • ttest(x,m) 在0.05的显著水平上检验矢量样本x的均值为m的假设（零假设），返回结果为0，表示接受零假设，为1，则拒绝零假设。 • h=ttest(x,m,alpha)自定义显著水平alpha,其余同上 • [h,sig,ci]=ttest(x,m,alpha,tail)　tail取0, 1, -1分别表 示备择假设为均值不等于，大于，小于m。（注意此时的零假设）sig为与t统计量有关的p值，ci为均值真值的1-alpha置信区间。 13/25

返回结果h为0，接受零假设；h为1，则拒绝零假设。 [例7-17] 某电子元件寿命x（以小时计）服从正态分布,μ和σ2均未知。现测得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 170]; >> [h,p,ci]=ttest(x,225,0.05,1) h = 　 0 p = 　0.2570 ci = 　 198.2321 Inf 在0.05的显著性水平上接受μ≤ 225的零假设 不能认为元件的平均寿命大于225小时 14/25

7.4.1 两个正态总体均值差的检验 对两个独立同方差（方差未知）正态总体的样本均 值差异进行t检验 函数：ttest2； 调用格式： • h=ttest2(x,y) • [h,significance,ci]=ttest2(x,y,alpha) • ttest2(x,y,alpha,tail)　 tail取0, 1, -1分别表示备择假设为μx≠ μy , μx> μy , μx<μy 。 [例7-18] 在平炉上进行一项试验以确定改变操作方法的建议是否会增加钢的得率，试验是在同一个平炉上进行的。每炼一炉钢，除操作方法外其它条件都尽可能做到相同。先用标准方法炼一炉，然后 15/25

μx≤ μy 用建议的新方法炼一炉，以后交替进行，各炼10炉，其钢的得率分别为： 标准方法　78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3 新方法　　79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1 设这两个样本相互独立，且分别来自正态总体，均值和方差都未知。问建议的新操作方法是否能 提高钢的得率？ >> x=[78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3]; >> y=[79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1]; >> [h,sig,ci]=ttest2(x,y,0.05,1) h = sig = 0.9998 ci = 　-4.4917 Inf μx≤ μy 16/25

q-q图：用qqplot函数生成两个样本的q-q（quan-tile分位数）图。若两样本来自同一分布，图中数 >> h=ttest2(x,y,0.05,0) h = 1 μx≠μy 综合以上得μx<μy ，新方法得钢率比标准方法高。 7.4.3 基于成对数据的检验 还是用ttest函数进行t检验，具体见P145例7-19 7.4.4 分布拟合检验 q-q图：用qqplot函数生成两个样本的q-q（quan-tile分位数）图。若两样本来自同一分布，图中数 据点呈直线关系，否则为曲线关系。调用格式为： • qqplot(X)：显示X的样本值与服从正态分布的理 论数据之间的q-q图。 17/25

生成一个服从威布尔分布的样本 [例7-18] • qqplot(X, Y)：显示X和Y两个样本的q-q图。 • h=qqplot(X, Y, pvec)：返回直线的句柄到h中。 [例7-18] 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) subplot(2,2,3) qqplot(z) subplot(2,2,4) qqplot(x,z) hold off 生成两个不同均值，方差的正态样本 生成一个服从威布尔分布的样本 18/25

由第一个子图 看出X服从正态 分布。 由第二个子图 看出X和Y可看作 同分布的。 由第三个子图 看出Z不服从正 态分布。 由第四个子图 　由第一个子图 看出X服从正态 分布。 　由第二个子图 看出X和Y可看作 同分布的。 　由第三个子图 看出Z不服从正 态分布。 　由第四个子图 看出X和Z不是同 分布的。 19/25

峰度-偏度检验(Jarque-Bera检验) 基于数据样本的偏度和峰度，评价给定数据服从未 知均值和方差正态分布的假设是否成立。 函数：jbtest,　调用格式为： • H=jbtest(X): 以0.05的显著水平对数据矢量X进行 Jarque-Bera检验，返回值H=0，接受X服从正态分 布的假设，H=1，则拒绝该假设。 • H=jbtest(X，alpha):　指定显著水平alpha • [H, P, JBSTAT, CV]=jbtest(X，alpha):　返回H值， 检验的p值P，检验统计量值JBSTAT和临界值CV。 例：见P147例7-21 注：Jarque-Bera检验不能用于小样本检验。 20/25

7.5 方差分析 7.5.1 单因子方差分析 研究不同因素及因素的不同水平对事件发生的影响程度 函数：anova1 7.5 方差分析 　　研究不同因素及因素的不同水平对事件发生的影响程度 7.5.1 单因子方差分析 函数：anova1　 调用格式：p = anova1(X, group, ‘displayopt’) 　　当X为矩阵时，每一列看成一组，定义group变量（字符数组或单元数组，长度等于X 的列数）做为各组标签。返回零假设（各组均值相等）成立的概率p。一般p小于0.05或0.01时，拒绝零假设，否则，接受零假设。（适用于平衡数据） 　　‘displayopt’参数为’on’(默认设置)时，显示anova表和箱型图，为’off’时不显示。 21/25

>> group = ｛‘a1’ , ‘a2’, ‘a3’, ‘a4’｝ >> p=anova1(x,group) 例 >> x =[73 44 18 31; 45 62 49 60; 35 65 40 17; 　　 14 54 46 62; 56 24 61 24; 66 35 17 58] >> group = ｛‘a1’ , ‘a2’, ‘a3’, ‘a4’｝ >> p=anova1(x,group) p = 0.7918　　　 认为列均值是相同的 22/25

23/25

7.5.1 双因子方差分析 当X为矢量时，通过输入变量group （字符数组或单元数组，长度与X 相同）进行标示分组。（适用于不平衡数据） 其它调用格式：p = anova2(X); p = anova2(X, group); [p, table] = anova1(…); [p, table,stats] = anova1(…) 7.5.1 双因子方差分析 函数：anova2　 调用格式： p = anova2(X, reps) 样本Ｘ中不同列中的数据代表因子Ａ的变化，不同行中数据代表因子Ｂ的变化。变量reps表示每个单元中观测值的个数。当reps=1（默认值）时，返回两个p值到p矢量中： 24/25

零假设H0A（源于因子Ａ的所有样本（X中列样本） 取自相同总体）的p值。 零假设H0B（源于因子B的所有样本（X中行样本） 当reps >1时，返回第三个p值到p矢量中： 零假设H0AB（因子Ａ和因子B间无交互效应）的p值。 若任一p值接近于0（一般判断时为小于0.05或0.01）, 则认为与之相关的零假设不成立。 其它调用格式： • p = anova2(X, group, ‘displayopt’) • [p, table] = anova2(…) • [p, table,stats] = anova2(…) 例　见P151例7-23, P152例7-24 25/25