(精品)matlab在科学计算中的应用8.ppt

《(精品)matlab在科学计算中的应用8.ppt》由会员分享，可在线阅读，更多相关《(精品)matlab在科学计算中的应用8.ppt（113页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第8章 概率论与数理统计问题的求解概率分布与伪随机数生成统计量分析数理统计分析方法及计算机实现统计假设检验方差分析及计算机求解 8.1概率分布与伪随机数生成 8.1.1 概率密度函数与分布函数概述图4-9通用函数计算概率密度函数值 函数 pdf格式 P=pdf(name，K，A)P=pdf(name，K，A，B)P=pdf(name，K，A，B，C)说明 返回在X=K处、参数为A、B、C的概率密度值，对于不同的分布，参数个数是不同；name为分布函数名。例如二项分布：设一次试验，事件Y发生的概率为p，那么，在n次独立重复试验中，事件Y恰好发生K次的概率P_K为：P_K=PX=K=pdf(bin

2、o，K，n，p)例:计算正态分布N（0，1）的随机变量X在点0.6578的密度函数值。解：pdf(norm,0.6578,0,1)ans=0.3213例:自由度为8的卡方分布，在点2.18处的密度函数值。解：pdf(chi2,2.18,8)ans=0.0363 随机变量的累积概率值(分布函数值)通用函数cdf用来计算随机变量的概率之和（累积概率值）函数 cdf格式 cdf(name，K，A)cdf(name，K，A，B)cdf(name，K，A，B，C)说明 返回以name为分布、随机变量XK的概率之和的累积概率值，name为分布函数名.例:求标准正态分布随机变量X落在区间(-，0.4)内的概

3、率。解：cdf(norm,0.4,0,1)ans=0.6554例：求自由度为16的卡方分布随机变量落在0，6.91内的概率。解：cdf(chi2,6.91,16)ans=0.0250随机变量的逆累积分布函数 MATLAB中的逆累积分布函数是已知，求x。命令 icdficdf 计算逆累积分布函数格式 icdf(name，F，A)icdf(name，F，A，B)icdf(name，F，A，B，C)说明 返回分布为name，参数为a1,a2,a3，累积概率值为F的临界值，这里name与前面相同。如果F=cdf(name，X，A，B，C)，则 X=icdf(name，F，A，B，C)例:在标准正态分布

4、表中，若已知F=0.6554,求X解：icdf(norm,0.6554,0,1)ans=0.3999例：公共汽车门的高度是按成年男子与车门顶碰头的机会不超过1%设计的。设男子身高X（单位：cm）服从正态分布N（175，6），求车门的最低高度。解：设h为车门高度，X为身高。求满足条件 FX=0.99，即 FXh h=icdf(norm,0.99,175,6)h=188.95818.1.2 常见分布的概率密度函数与分布函数 8.1.2.1 Poisson分布其要求x是正整数。其中：x为选定的一组横坐标向量，y为x各点处的概率密度函数值。例：绘制 l l=1,2,5,10 时 Poisson 分布的

5、概率密度函数与概率分布函数曲线。x=0:15;y1=;y2=;lam1=1,2,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)8.1.2.2 正态分布正态分布的概率密度函数为:例：x=-5:.02:5;y1=;y2=;mu1=-1,0,0,0,1;sig1=1,0.1,1,10,1;sig1=sqrt(sig1);for i=1:length(mu1)y1=y1,normpdf(x,mu1(i),sig1(i);y2=y2,no

6、rmcdf(x,mu1(i),sig1(i);end plot(x,y1),figure;plot(x,y2)8.1.2.3 分布例：绘制 为（1,1),(1,0.5),(2,1),(1,2),(3,1)时 x=-0.5:.02:5;x=-eps:-0.02:-0.5,0:0.02:5;x=sort(x);替代 y1=;y2=;a1=1,1,2,1,3;lam1=1,0.5,1,2,1;for i=1:length(a1)y1=y1,gampdf(x,a1(i),lam1(i);y2=y2,gamcdf(x,a1(i),lam1(i);end plot(x,y1),figure;plot(x,

7、y2)8.1.2.4 分布（卡方分布）其为一特殊的 分布，a=k/2,l l=1/2。例：x=-eps:-0.02:-0.5,0:0.02:2;x=sort(x);k1=1,2,3,4,5;y1=;y2=;for i=1:length(k1)y1=y1,chi2pdf(x,k1(i);y2=y2,chi2cdf(x,k1(i);end plot(x,y1),figure;plot(x,y2)8.1.2.5 分布概率密度函数为:其为参数k的函数，且k为正整数。例：x=-5:0.02:5;k1=1,2,5,10;y1=;y2=;for i=1:length(k1)y1=y1,tpdf(x,k1(i

8、);y2=y2,tcdf(x,k1(i);end plot(x,y1),figure;plot(x,y2)8.1.2.6 Rayleigh分布例：x=-eps:-0.02:-0.5,0:0.02:5;x=sort(x);b1=.5,1,3,5;y1=;y2=;for i=1:length(b1)y1=y1,raylpdf(x,b1(i);y2=y2,raylcdf(x,b1(i);end plot(x,y1),figure;plot(x,y2)8.1.2.7 F 分布其为参数p,q的函数，且p,q均为正整数。例：分别绘制（p,q）为（1,1),(2,1),(3,1)(3,2),(4,1)时F分

9、布的概率密度函数与分布函数曲线。x=-eps:-0.02:-0.5,0:0.02:1;x=sort(x);p1=1 2 3 3 4;q1=1 1 1 2 1;y1=;y2=;for i=1:length(p1)y1=y1,fpdf(x,p1(i),q1(i);y2=y2,fcdf(x,p1(i),q1(i);end plot(x,y1),figure;plot(x,y2)8.1.3 概率问题的求解图4-9例：b=1;p1=raylcdf(0.2,b);p2=raylcdf(2,b);P1=p2-p1P1=0.8449 p1=raylcdf(1,b);P2=1-p1P2=0.6065例：syms

10、 x y;f=x2+x*y/3;P=int(int(f,x,0,1/2),y,0,1/2)P=5/192 syms x y;f=x2+x*y/3;P=int(int(f,x,0,1),y,0,2)P=18.1.4 随机数与伪随机数例：b=1;p=raylrnd(1,30000,1);xx=0:.1:4;yy=hist(p,xx);hist()找出随机数落入各个子区间的点个数，并由之拟合出生成数据的概率密度。yy=yy/(30000*0.1);bar(xx,yy),y=raylpdf(xx,1);line(xx,y)8.2 统计量分析 8.2.1 随机变量的均值与方差例：均值 syms x;sy

11、ms a lam positive p=lama*x(a-1)/gamma(a)*exp(-lam*x);m=int(x*p,x,0,inf)m=1/lam*a 方差 s=simple(int(x-1/lam*a)2*p,x,0,inf)s=a/lam2已知一组随机变量样本数据构成的向量:求该向量各个元素的均值、方差和标准差、中位数medianmedian例：生成一组 30000 个正态分布随机数，使其均值为 0.5，标准差为1.5，分析数据实际的均值、方差和标准差，如果减小随机变量个数，会有什么结果？p=normrnd(0.5,1.5,30000,1);mean(p),var(p),std(

12、p)ans=0.4879 2.2748 1.5083300个随机数 p=normrnd(0.5,1.5,300,1);mean(p),var(p),std(p)ans=0.4745 1.9118 1.3827可见在进行较精确的统计分析时不能选择太小的样本点。例：m,s=raylstat(0.45)m=0.5640s=0.08698.2.2 随机变量的矩例：求解原点矩 syms x;syms a lam positive;p=lama*x(a-1)/gamma(a)*exp(-lam*x);for n=1:5,m=int(xn*p,x,0,inf),endm=1/lam*a m=1/lam2*a

13、*(a+1)m=1/lam3*a*(a+1)*(a+2)m=1/lam4*a*(a+1)*(a+2)*(a+3)m=1/lam5*a*(a+1)*(a+2)*(a+3)*(a+4)有规律 syms n;m=simple(int(x)n*p,x,0,inf)直接求出m=lam(-n)*gamma(n+a)/gamma(a)for n=1:6,s=simple(int(x-1/lam*a)n*p,x,0,inf),end 中心距s=0s=a/lam2 s=2*a/lam3s=3*a*(a+2)/lam4s=4*a*(5*a+6)/lam5s=5*a*(3*a2+26*a+24)/lam6 好像无规

14、律例：考虑前面的随机数，可以用下面的语句得出随机数的各阶矩。A=;B=;p=normrnd(0.5,1.5,30000,1);n=1:5;for r=n,A=A,sum(p.r)/length(p);B=B,moment(p,r);end A,BA=0.5066 2.4972 3.5562 18.7530 41.5506B=0 2.2405 0.0212 15.1944 0.0643求各阶距的理论值：syms x;A1=;B1=;p=1/(sqrt(2*pi)*1.5)*exp(-(x-0.5)2/(2*1.52);for i=1:5 A1=A1,vpa(int(xi*p,x,-inf,inf

15、),12);B1=B1,vpa(int(x-0.5)i*p,x,-inf,inf),12);end A1,B1A1=.500000000001,2.50000000000,3.50000000001,18.6250000000,40.8125000000 B1=0,2.25000000000,0,15.1875000000,08.2.3 多变量随机数的协方差分析例：p=randn(30000,4);cov(p)ans=1.0033 0.0131 0.0036 0.0020 0.0131 1.0110 0.0061 -0.0154 0.0036 0.0061 1.0055 -0.0004 0.0

16、020 -0.0154 -0.0004 0.98818.2.4 多变量正态分布的联合概率密度即分布函数例：mu1=-1,2;Sigma2=1 1;1 3;%输入均值向量和协方差矩阵 X,Y=meshgrid(-3:0.1:1,-2:0.1:4);xy=X(:)Y(:);%产生网格数据并处理(两列2501*2）p=mvnpdf(xy,mu1,Sigma2);%求取联合概率密度 P=reshape(p,size(X);Change size(2501*161*41)surf(X,Y,P)对协方差矩阵进行处理，可计算出新的联合概率密度函数。Sigma2=diag(diag(Sigma2);%消除协方

17、差矩阵的非对角元素 p=mvnpdf(xy,mu1,Sigma2);P=reshape(p,size(X);surf(X,Y,P)R为m行n列。例：mu1=-1,2;Sigma2=1 1;1 3;R1=mvnrnd(mu1,Sigma2,2000);plot(R1(:,1),R1(:,2),o)Sigma2=diag(diag(Sigma2);figure;R2=mvnrnd(mu1,Sigma2,2000);plot(R2(:,1),R2(:,2),o)8.3数理统计分析方法及计算机实现 8.3.1 参数估计与区间估计 无论总体X的分布函数 的类型已知或未知，我们总是需要去估计某些未知参数或

18、数字特征，这就是参数估计问题.即参数估计就是从样本 出发，构造一些统计量 （i=1，2，k）去估计总体X中的某些参数（或数字特征）（i=1，2，k）.这样的统计量称为估计量估计量.1、点估计、点估计：构造 的函数 作为参数 的点估计量，称统计量 为总体X参数 的点估计量.2.区间估计区间估计：构造两个函数 (X1，X2，Xn）和 (X1，X2，Xn）做成区间，把这 作为参数 的区间估计.区间估计的求法区间估计的求法 设总体X的分布中含有未知参数 ，若对于给定的概率 ，存在两个统计量 (X1，X2，Xn）和 (X1，X2，Xn）,使得 则称随机区间 为参数 的置信水平为 的置信区间置信区间,称

19、为置信下限置信下限,称 为置信上限置信上限.由极大似然法估计出该分布的均值、方差 及其置信区间。置信度越大，得出的置信区间越小，即得出的结果越接近于真值。还有gamfit(),raylfit(),poissfit()，unifit()（均匀分布）等参数估计函数例：p=gamrnd(1.5,3,30000,1);Pv=0.9,0.92,0.95,0.98;A=;for i=1:length(Pv)a,b=gamfit(p,Pv(i);A=A;Pv(i),a(1),b(:,1),a(2),b(:,2)end AA=0.9000 1.5137 1.5123 1.5152 2.9825 2.9791

20、2.9858 0.9200 1.5137 1.5126 1.5149 2.9825 2.9798 2.9851 0.9500 1.5137 1.5130 1.5144 2.9825 2.9808 2.9841 0.9800 1.5137 1.5135 1.5140 2.9825 2.9818 2.9831 num=300,3000,30000,300000,3000000;A=;for i=1:length(num)p=gamrnd(1.5,3,num(i),1);a,b=gamfit(p,0.95);A=A;num(i),a(1),b(:,1),a(2),b(:,2);end A(:,2,3

21、,4,5,6,7)ans=1.4795 1.4725 1.4865 2.9129 2.8960 2.9299 1.4218 1.4198 1.4238 3.1676 3.1623 3.1729 1.4898 1.4891 1.4904 3.0425 3.0409 3.0442 1.4998 1.4996 1.5000 3.0054 3.0049 3.0059 1.5006 1.5005 1.5007 2.9968 2.9966 2.9969 要达到参数估计效果良好，随机数不能选得太少，也不能选得太多，此例中为30000为好。8.3.2 多元线性回归与区间估计例：a=1-1.232 2.23 2

22、 4 3.792;X=randn(120,6);y=X*a;a1=inv(X*X)*X*y;a1ans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920 a,aint=regress(y,X,0.02);a,aintans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920ans=1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920 1.0000 -1.2320 2.2300 2.0000 4.0000 3.7920 yhat=y+sqrt(0.5)*randn(120,1);a,aint=r

23、egress(yhat,X,0.02);a,aint a=1-1.232 2.23 2 4 3.792ans=1.0576 -1.3280 2.1832 2.0151 4.0531 3.7749ans=0.8800 -1.5107 2.0284 1.8544 3.8788 3.6221 1.2353 -1.1453 2.3379 2.1757 4.2274 3.9276 errorbar(1:6,a,aint(:,1)-a,aint(:,2)-a)errorbar()用图形绘制参数估计的置信区间。yhat=y+sqrt(0.1)*randn(120,1);a,aint=regress(yhat

24、,X,0.02);errorbar(1:6,a,aint(:,1)-a,aint(:,2)-a)8.3.3 非线性函数的最小二乘参数估计与区间估计r为参数下的残差构成的向量。J为各个Jacobi行向量构成的矩阵。例：f=inline(a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x),a,x);x=0:0.1:10;y=f(0.12,0.213,0.54,0.17,1.23,x);a,r,j=nlinfit(x,y,f,1;1;1;1;1);aans=0.11999999763418 0.21299999458274 0.5400000019681

25、8 0.17000000068705 1.22999999996315 ci=nlparci(a,r,j)0.12,0.213,0.54,0.17,1.23ci=0.11999999712512 0.11999999814323 0.21299999340801 0.21299999575747 0.54000000124534 0.54000000269101 0.17000000036077 0.17000000101332 1.22999999978603 1.23000000014028 y=f(0.12,0.213,0.54,0.17,1.23,x)+0.02*rand(size(x

26、);a,r,j=nlinfit(x,y,f,1;1;1;1;1);aans=0.12655784086874 0.17576593556541 0.54363873794463 0.17129712329146 1.23139632101927 ci=nlparci(a,r,j)ci=0.12240417108574 0.130711510651740.16754837168468 0.183983499446140.53737093469422 0.549906541195040.16845014477426 0.174144101808661.22983289563708 1.232959

27、74640145 errorbar(1:5,a,ci(:,1)-a,ci(:,2)-a)例：a=1;1;1;1;1;1;f=inline(a(1)*x(:,1).3+a(2).*sin(a(3)*x(:,2)，.*x(:,3)+(a(4)*x(:,3).3+a(5)*x(:,3)+a(6),a,x);X=randn(120,3);y=f(a,X)+sqrt(0.2)*randn(120,1);ahat,r,j=nlinfit(X,y,f,0;2;3;2;1;2);ahatahat=0.99166464884539 1.04776526972943 0.97668595800756 1.0202

28、2345889541 0.88639528713563 1.09317291667891 ci=nlparci(ahat,r,j);ci 置信区间ci=0.89133624667624 1.09199305101455 0.86664749663205 1.22888304282680 0.83628948119418 1.11708243482094 0.98466523279168 1.05578168499914 0.73055684224143 1.04223373202984 0.99932407018303 1.18702176317478 errorbar(1:6,ahat,ci

29、(:,1)-ahat,ci(:,2)-ahat)y1=f(ahat,X);plot(y y1)绘制曲线8.4 统计假设检验8.4.1 正态分布的均值假设检验 H为假设检验的结论，当H0时表示不拒绝H0假设，否则表示拒绝该假设。s为接受假设的概率值，为其均值的置信区间。若未知正态分布的标准差时，可用此函数。例：设某车间用一台包装机包装葡萄糖，包得的袋装糖重量是一个随机数，它服从正态分布。当机器正常时，其均值为0.5公斤，标准差为0.015。某日开工后检验包装机是否正常，随机地抽取它所包装的的糖9袋，称得净重为（公斤）0.497,0.506,0.518,0.524,0.498,0.511,0.52

30、,0.515,0.512,问机器是否正常？解：（分析）总体均值、标准差已知，则可设样本的标准差为0.015，于是 问题就化为根据样本值来判断 还是 。为此提出假设：x=0.497,0.506,0.518,0.524,0.498,0.511,0.52,0.515,0.512;H,p,ci=ztest(x,0.5,0.015,0.05)H=1p=0.0248%样本观察值的概率 ci=0.5014 0.5210%包含均值95置信区间，均值0.5在此区间之外 结果H1，说明在0.05的水平下，拒绝原假设，即认为这天包装机工作不正常。例：某种电子元件的寿命X（以小时计）服从正态分布，均值、方差均未知。现

31、测得16只元件的寿命如下：159 280 101 212 224 379 179 264 222 262 168 250 149 260 485 170,问是否有理由认为元件的平均寿命大于225（小时）：解：按题意需做如下假设：取 x=159 280 101 212 224 379 179 264 222 262 168 250 149 260 485 170;H,p,ci=ttest(x,225,0.05)H=0p=0.6677ci=185.3622 285.1378%均值225在该置信区间内 结果表明，H0，即在显著水平为0.05的情况下，不能拒绝原假设。即认为元件的平均寿命不大于225小

32、时。8.4.2 正态分布假设检验 由随机样本判定分布是否为正态分布，可用下面两个假设算法的函数。s为接受假设的概率值，s越接近于0，则可以拒绝是正态分布的原假设.?例：X=216,203,197,208,206,209,206,208,202,203,206,213,218,207,208,.202,194,203,213,211,193,213,208,208,204,206,204,206,208,209,.213,203,206,207,196,201,208,207,213,208,210,208,211,211,214,.220,211,203,216,224,211,209,218

33、,214,219,211,208,221,211,218,.218,190,219,211,208,199,214,207,207,214,206,217,214,201,212,.213,211,212,216,206,210,216,204,221,208,209,214,214,199,204,.211,201,216,211,209,208,209,202,211,207,202,205,206,216,206,.213,206,207,200,198,200,202,203,208,216,206,222,213,209,219;H,p=jbtest(X,0.05)%P为接受假设的概

34、率值，P越接近于0，则可以拒绝是正态分布的原假设；H=0 p=0.7281 mu1,sig1,mu_ci,sig_ci=normfit(X,0.05);mu=mu1,mu_cimu=208.8167 207.6737 209.9596该分布的均值及置信区间 sig=sig1,sig_cisig=6.3232 5.6118 7.2428该分布的方差及置信区间例：r=gamrnd(1,3,400,1);H,p,c,d=jbtest(r,0.05)H=1p=0c=504.2641d=5.9915%P为接受假设的概率值，P越接近于0，则可以拒绝是正态分布的原假设；c为测试统计量的值，d为是否拒绝原假设

35、的临界值,cd,故拒绝。8.4.3 其它分布的Kolmogorov-Smirnov 检验 此函数（Kolmogorov-Smirnov 算法）可对任意已知分布函数进行有效的假设检验。其中cdffun为两列的值，第一列为自变量，第二列为对应的分布函数的值。例：r=gamrnd(1,3,400,1);alam=gamfit(r)alam=0.9708 3.1513检验：r=sort(r);H0,p=kstest(r,r gamcdf(r,alam(1),alam(2),0.05)H0=0p=0.6067 8.5方差分析及计算机求解 8.5.1 单因子方差分析 对一些观察来说，只有一个外界因素可能对

36、观测的现象产生影响。单因素方差分析是比较两组或多组数据的均值，它返回原假设均值相等的概率,若p值接近于0，则原假设受到怀疑，说明至少有一列均值与其余列均值有明显不同。X为需要分析的数据，每一列对应于随机分配的一个组的测试数据，这样会返回概率p，tab为方差分析表。stats为统计结果量，为结构变量，包括每组均值等。单因子方差分析表例：建立A矩阵，并求各列的均值。A=5,4,6,7,9;8,6,4,4,3;7,6,4,6,5;7,3,5,6,7;10,5,4,3,7;8,6,3,5,6;mean(A)ans=7.5000 5.0000 4.3333 5.1667 6.1667 p,tbl,sta

37、ts=anova1(A)%单因子方差分析p=0.0136%F Columns 36.4667 4 9.1167 3.8960 0.0136 Error 58.5000 25 2.3400 Total 94.9667 29 stats=gnames:5x1 char n:6 6 6 6 6 source:anova1 means:7.5000 5 4.3333 5.1667 6.1667 df:25 s:1.5297有显著影响。盒式图 单因子方差表A=5,4,6,7,9;8,6,4,4,3;7,6,4,6,5;7,3,5,6,7;10,5,4,3,7;8,6,3,5,6;例：设有3台机器，用来生

38、产规格相同的铝合金薄板。取样测量薄板的厚度，精确至厘米。得结果如下：机器1：0.236 0.238 0.248 0.245 0.243机器2：0.257 0.253 0.255 0.254 0.261机器3：0.258 0.264 0.259 0.267 0.262检验各台机器所生产的薄板的厚度有无显著的差异？X=0.236 0.238 0.248 0.245 0.243;0.257 0.253 0.255 0.254 0.261；0.258 0.264 0.259 0.267 0.262;P=anova1(X)P=1.3431e-0058.5.2 双因子方差分析 如果有两种因子可能影响到某现

39、象的统计规律，则应该引入双因子方差分析的概念。这时观察值可表示为一个三维数组。根据双因子的特点，可以引入3个假设双因素方差表表中记号的定义求解双因子方差分析问题:例：比较 3 种松树在4 个不同地区的生长情况有无差别,在每个地区对每种松树随机地选择 5 株，测量它们的胸径，对它们进行双因子方差分析。B=23,15,26,13,21,25,20,21,16,18,21,17,16,24,27,14,17,19,20,24;28,22,25,19,26,30,26,26,20,28,19,24,19,25,29,17,21,18,26,23;18,10,12,22,13,15,21,22,14,1

40、2,23,25,19,13,22,16,12,23,22,19;anova2(B,5);5表示每一单元观察点的数目 小，树种对观测树的胸径有显著影响；很大，所以没有理由拒绝另外两个假设。故得出结论：地区对树的胸径无显著影响，不同区域对不同树种的胸径观测结果也无显著影响。计算均值：C=;for i=1:3 for j=1:4 C(i,j)=mean(B(i,1:5+(j-1)*5);end,end C=C;mean(C);C=C mean(C)C=19.6000 20.0000 21.0000 18.8000 19.8500 24.0000 26.0000 23.2000 21.0000 23.

41、5500 15.0000 16.8000 20.4000 18.4000 17.6500 19.5333 20.9333 21.5333 19.4000 20.3500例、下表记录了3 位操作工分别在4 台不同机器上操作的日产量，试检验 操作工之间的差异是否显著?机器之间的差异是否显著?交互作用是否显著()?机 操作工 机 操作工器 1 2 3 器 1 2 3M1 15 15 17 19 19 16 16 18 21 M3 15 17 16 18 17 16 18 18 18M2 17 17 17 15 15 15 19 22 22 M4 18 20 22 15 16 17 17 17 17解

42、：A=15,15,17,19,19,16,16,18,21;17,17,17,15,15,15,19,22,22;15,17,16,18,17,16,18,18,18;18,20,22,15,16,17,17,17,17;p,tbl=anova2(A,3)p=0.6645 0.0023 0.0002tbl=Source SS df MS F ProbF Columns 2.7500 3 0.9167 0.5323 0.6645 Rows 27.1667 2 13.5833 7.8871 0.0023 Interaction 73.5000 6 12.2500 7.1129 1.9217e-00

43、4 Error 41.3333 24 1.7222 Total 144.7500 35 由上表可知：机器自由度（df）是（41）3，p值为0.6645,故机器之间的差异不显著。操作工自由度（df）是（31）2，p值为0.0023,故操作工之间的差异显著。交互作用自由度（df）是（41）（31）6，p值为1.9217e-004，故交互作用差异非常显著。8.5.3 多因子方差分析8.6统计作图统计作图8.6.1正整数的频率表正整数的频率表命令 正整数的频率表函数 tabulate格式 table=tabulate(X)%X为正整数构成的向量，返回3列：第1列中包含X的值,第2列为这些值的个数，第3

44、列为这些值的频率。例 A=1 2 2 5 6 3 8A=1 2 2 5 6 3 8 tabulate(A)Value Count Percent 1 1 14.29%2 2 28.57%3 1 14.29%4 0 0.00%5 1 14.29%6 1 14.29%7 0 0.00%8 1 14.29%8.6.2经验累积分布函数图形经验累积分布函数图形函数 cdfplot格式 cdfplot(X)%作样本X（向量）的累积分布函数图形h=cdfplot(X)%h表示曲线的环柄图h,stats=cdfplot(X)%stats表示样本的一些特征例 X=normrnd(0,1,50,1);h,stat

45、s=cdfplot(X)h=151.0077stats=min:-2.3646%样本最小值 max:2.3093%最大值 mean:-0.0298%平均值 median:-0.0084%中间值 std:1.0221%样本标准差 8.6.3最小二乘拟合直线最小二乘拟合直线函数 lsline格式 lsline%最小二乘拟合直线 h=lsline%h为直线的句柄例 X=2 3.4 5.6 8 11 12.3 13.8 16 18.8 19.9;plot(X,+)lsline8.6.4绘制正态分布概率图形绘制正态分布概率图形函数 normplot格式 normplot(X)%若X为向量，则显示正态分布

46、概率图形，若X为矩阵，则显示每一列的正态分布概率图形。h=normplot(X)%返回绘图直线的句柄说明 样本数据在图中用“+”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可能在图中产生弯曲。例:X=normrnd(0,1,50,1);normplot(X)8.6.5绘制威布尔绘制威布尔(Weibull)概率图形概率图形函数 weibplot格式 weibplot(X)%若X为向量，则显示威布尔(Weibull)概率图形，若X为矩阵，则显示每一列的威布尔概率图形。h=weibplot(X)%返回绘图直线的柄说明 绘制威布尔(Weibull)概率图形的目的是用图解法估计来自威布尔分布

47、的数据X，如果X是威布尔分布数据，其图形是直线的，否则图形中可能产生弯曲。例:r=weibrnd(1.2,1.5,50,1);weibplot(r)8.6.6样本数据的盒图样本数据的盒图函数 boxplot格式 boxplot(X)%产生矩阵X的每一列的盒图和“须”图，“须”是从盒的尾部延伸出来，并表示盒外数据长度的线，如果“须”的外面没有数据，则在“须”的底部有一个点。boxplot(X,notch)%当notch=1时，产生一凹盒图，notch=0时产生一矩箱图。boxplot(X,notch,sym)%sym表示图形符号，默认值为“+”。boxplot(X,notch,sym,vert)

48、%当vert=0时，生成水平盒图，vert=1时，生成竖直盒图（默认值vert=1）。boxplot(X,notch,sym,vert,whis)%whis定义“须”图的长度，默认值为1.5，若whis=0则boxplot函数通过绘制sym符号图来显示盒外的所有数据值。例:x1=normrnd(5,1,100,1);x2=normrnd(6,1,100,1);x=x1 x2;boxplot(x,1,g+,1,0)8.6.7 给当前图形加一条参考线 函数 refline格式 refline(slope,intercept)%slope表示直线斜率，intercept表示截距refline(slo

49、pe)slope=a b，图中加一条直线：y=b+ax。例:y=3.2 2.6 3.1 3.4 2.4 2.9 3.0 3.3 3.2 2.1 2.6;plot(y,+)refline(0,3)8.6.8在当前图形中加入一条多项式曲线在当前图形中加入一条多项式曲线函数 refcurve格式格式h=refcurve(p)%在图中加入一条多项式曲线，在图中加入一条多项式曲线，h为曲线的环柄，为曲线的环柄，p为多项式系数向量，为多项式系数向量，p=p1,p2,p3,pn，其中其中p1为最高幂项系数。为最高幂项系数。例:火箭的高度与时间图形，加入一条理论高度曲线，火箭初速为100m/秒。h=85 16

50、2 230 289 339 381 413 437 452 458 456 440 400 356;plot(h,+)refcurve(-4.9 100 0)8.6.9样本的概率图形样本的概率图形函数 capaplot格式 p=capaplot(data,specs)%data为所给样本数据，specs指定范围，p表示在指定范围内的概率。说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率例:data=normrnd(0,1,30,1);p=capaplot(data,-2,2)p=0.97938.6.10附加有正态密度曲线的直方图附加有正态密度曲线的直方图函数 histfit格式 h