1、Matlab在统计分析中的应用,MATLAB概述与运算基础,MATLAB语言：广泛应用于工程计算及数值分析领域的新型高级语言 1984年由美国 MathWorks 公司推向市场 已成为国际公认的最优秀的工程应用开发环境 功能强大、简单易学、编程效率高,在欧美各高等院校，MATLAB已成为线性代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等课程的基本教学工具，成为大学生、硕士生以及博士生必须掌握的基本技能。,MATLAB：MATrix LABorotory,一. MATLAB特点:,一. MATLAB特点:,1. 数值计算和符号计算功能 数值计算功能包括：矩阵运算、多项式
2、和有理分式运算、数据统计分析、数值积分、优化处理等。 符号计算将得到问题的解析解。,2.MATLAB语言 除了命令行的交互式操作以外，还可以程序方式工作。使用MATLAB很容易地实现C或FORTRAN语言的几乎全部功能。,3.图形功能 两个层次的图形命令：一种是对图形句柄进行的低级图形命令，另一种是建立在低级图形命令之上的高级图形命令。利用MATLAB的高级图形命令可以轻而易举地绘制二维、三维乃至四维图形，并可进行图形和坐标的标识、视角和光照设计、色彩精细控制等等。,4.应用工具箱 包括：基本部分和各种可选的工具箱。 基本部分中有数百个内部函数。 工具箱分为两大类：功能性工具箱和学科性工具箱。
3、 a. 功能性工具箱：主要用来扩充符号计算功能、可视建模仿真功能及文字处理功能等。 b. 学科性工具箱：专业性较强，如控制系统工具箱、信号处理工具箱、神经网络工具箱、最优化工具箱、金融工具箱等，用户可以直接利用这些工具箱进行相关领域的科学研究。,命令窗口的命令编辑区用于输入命令和显示计算结果。例：键入 A=1 2 3;4 5 6;7 8 9输出 A = 1 2 34 5 67 8 9,二. MATLAB集成环境,矩阵运算,【例】 求解线性方程组AX=B 1 1.5 2 9 70 3.6 0.5 -4 4其中A= 7 10 -3 22 333 7 8.5 21 63 8 0 90 -203-4B
4、= 20516,在MATLAB命令窗口输入命令： a=1,1.5,2,9,7; 0,3.6,0.5,-4,4; 7,10,-3,22,33; 3,7,8.5,21,6; 3,8,0,90,-20;b=3; -4; 20; 5; 16;x=ab,得到的结果是：x =3.5653-0.9255-0.26950.14350.0101,【例】 求方程 x4+7x3 +9x-20=0的全部根。 命令窗口输入：p=1,7,0,9,-20; %建立多项式系数向量x=roots(p) %求根,得到的结果是：x =-7.2254-0.4286 + 1.5405i-0.4286 - 1.5405i1.0826,二
5、维图形,一、 plot函数 函数格式：plot(x,y) 其中x和y为坐标向量 函数功能：以向量x、y为轴，绘制曲线。【例】 在区间0X2内，绘制正弦曲线Y=SIN（X）x=0:pi/100:2*pi; y=sin(x); plot(x,y),【例】同时绘制正、余弦两条曲线Y1=SIN(X)和Y2=COS(X) x=0:pi/100:2*pi; y1=sin(x); y2=cos(x); plot(x,y1,x,y2),（一）线型与颜色 格式：plot(x,y1,cs,.) 其中c表示颜色， s表示线型,一、 plot函数,【例】 用不同线型和颜色重新绘制上例图形x=0:pi/100:2*pi
6、; y1=sin(x); y2=cos(x); plot(x,y1,go,x,y2,b-.) 其中：参数go和b-.表示图形的颜色和线型。 g：绿色 o：线型为圆圈 b：蓝色 -.：表示图形线型为点划线,（二）图形标记title(加图形标题); xlabel(加X轴标记); ylabel(加Y轴标记); text(X,Y,添加文本);,一、 plot函数,三维图形,一、 plot3函数 函数格式：plot3(x1,y1,z1,c1,x2,y2,z2,c2,) 其中x1,y1,z1表示三维坐标向量；c1,c2表示线形或颜色。 函数功能：以向量x，y，z为坐标，绘制三维曲线。,【例】 绘制三维螺旋
7、曲线t=0:pi/50:10*pi; x=sin(t),y=cos(t); plot3(x,y,t); title(helix),text(0,0,0,origin); xlabel(sin(t),ylabel(cos(t),zlabel(t); grid;,二、mesh函数 绘制三维网格图。函数格式：mesh(x,y,z,c) 其中：x，y控制X和Y轴坐标矩阵z是由(x，y)求得Z轴坐标(x,y,z)组成三维空间的网格点c用于控制网格点颜色,【例】 绘制三维网格曲面图 x=0:0.15:2*pi; y=0:0.15:2*pi; z=sin(y)*cos(x); %矩阵相乘 mesh(x,y,
8、z);,【例】画出由函数形成的立体网状图：,x=linspace(-2, 2, 25); % 在x轴上取25点 y=linspace(-2, 2, 25); % 在y轴上取25点 xx,yy=meshgrid(x,y); % xx和yy都是21x21的矩阵 zz=xx.*exp(-xx.2-yy.2); % 计算函数值，zz也是21x21的矩阵 mesh(xx, yy, zz); % 画出立体网状图,三、surf函数 绘制三维曲面图，各线条之间的补面用颜色填充。surf函数和mesh函数的调用格式一致。 函数格式: surf (x,y,z) 其中x，y控制X和Y轴坐标，矩阵z是由x，y求得的曲
9、面上Z轴坐标。,【例】 绘制三维曲面图 x=0:0.15:2*pi; y=0:0.15:2*pi; z=sin(y)*cos(x); %矩阵相乘 surf(x,y,z); xlabel(x-axis),ylabel(y-axis),zlabel(z-axis); title(3-D surf);,【例】剔透玲珑球 X0,Y0,Z0=sphere(30); %产生单位球面的三维坐标 X=2*X0;Y=2*Y0;Z=2*Z0; %产生半径为2的球面的三维坐标 surf(X0,Y0,Z0); %画单位球面 shading interp %采用插补明暗处理 hold on; mesh(X,Y,Z);h
10、old off %画外球面 hidden off %产生透视效果 axis off %不显示坐标轴,【例】卫星返回地球的运动轨线示意。 R0=1; %以地球半径为一个单位 a=12*R0;b=9*R0;T0=2*pi; %T0是轨道周期 T=5*T0;dt=pi/100;t=0:dt:T;f=sqrt(a2-b2); %地球与另一焦点的距离 th=12.5*pi/180; %卫星轨道与x-y平面的倾角 E=exp(-t/20); %轨道收缩率 x=E.*(a*cos(t)-f);y=E.*(b*cos(th)*sin(t);z=E.*(b*sin(th)*sin(t); plot3(x,y,z
11、,g) %画全程轨线 X,Y,Z=sphere(30);X=R0*X;Y=R0*Y;Z=R0*Z; %获得单位球坐标 grid on,hold on,surf(X,Y,Z),shading interp %画地球 x1=-18*R0;x2=6*R0;y1=-12*R0;y2=12*R0;z1=-6*R0;z2=6*R0; axis(x1 x2 y1 y2 z1 z2) %确定坐标范围 view(117 37),comet3(x,y,z,0.02),hold off %设视角、画运动轨线,五、等高线图 函数contour3 【例】 多峰函数peaks的等高线图 x,y,z=peaks(30);
12、contour3(x,y,z,16); xlabel(x-axis),ylabel(y-axis),zlabel(z-axis); title(contour3 of peaks),4.6 动画设计 动画功能函数：getframe、moviein和movie,【例】 播放一个不断变化的眼球程序。 m=moviein(20); %建立一个20个列向量组成的矩阵 for j=1:20plot(fft(eye(j+10) %绘制出每一幅眼球图并保存到m矩阵中m(:,j)=getframe; end movie(m,10);%以每秒10幅的速度播放画面,【例 】求下列三阶线性代数方程组的近似解,MAT
13、LAB程序为： A=2 -5 4;1 5 -2;-1 2 4； b=5;6;5； x=Ab,解法1：分别解方程组 (1)Ax=b1；(2)Ay=b2 A=1 -1 1;5 -4 3;2 1 1; b1=2;-3;1; b2=3;4;-5; x=Ab1 x =-3.80001.40007.2000,y=Ab2-3.6000-2.20004.4000,得两个线性代数方程组的解：(1) x1= -3.8， x2= 1.4， x3= 7.2； (2) y1= -3.6， y2=2.2， y3= 4.4,解法2：将两个方程组连在一起求解：Az=b b=2 3;-3 4;1 -5 z=Ab z =-3.8
14、000 -3.60001.4000 -2.20007.2000 4.4000,一、 基本统计处理,1、查取最大值 MAX函数的命令格式： Y,I= max (X)：将max(X)返回矩阵X的各列中的最大元素值及其该元素的位置赋予行向量Y与I；当X为向量时，则Y与I为单变量。 Y,I=max(X, ,DIM)：按数组X的第DIM维的方向查取其最大的元素值及其该元素的位置赋予向量Y与I。,【例】查找下面数列x的最大值。 x=3 5 9 6 1 8 % 产生数列x x = 3 5 9 6 1 8 y=max(x) % 查出数列x中的最大值赋予y y = 9 y,l=max(x) % 查出数列x中的最
15、大值及其该元素的位置赋予y,l y = 9 l = 3,【例】分别查找下面34的二维数组x中各列和各行元素中的最大值。 x=1 8 4 2;9 6 2 5;3 6 7 1 % 产生二维数组x x = 1 8 4 29 6 2 53 6 7 1 y=max(x) % 查出二维数组x中各列元素的最大值产生赋予行向量y y = 9 8 7 5,2、查取最小值MIN函数 3、求中值 Y=median(X)：将median(X)返回矩阵X各列元素的中值赋予行向量Y。若X为向量，则Y为单变量。 Y=median(X,DIM)：按数组X的第DIM维方向的元素求其中值赋予向量Y。若DIM=1，为按列操作；若D
16、IM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。,【例】分别求下面数列x1与x2的中值。 x1=9 -2 5 7 12; % 奇数个元素 y1=median(x) y1 =7 x2=9 -2 5 6 7 12; % 偶数个元素 y2=median(x) y2 =6.5000,4、求和 命令格式： Y=sum(X)：将sum(X)返回矩阵X各列元素之和赋予行向量Y；若X为向量，则Y为单变量。 Y=sum(X,DIM)：按数组X的第DIM维的方向的元素求其和赋予Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则
17、Y为单变量。,例： x=4 5 6;1 4 8 x =4 5 61 4 8 y=sum(x,1) y =5 9 14 y=sum(x,2) y =1513,5、求平均值 MEAN函数调用的命令格式： Y= mean(X)：将mean (X)返回矩阵X各列元素之的平均值赋予行向量Y。若X为向量，则Y为单变量。 Y= mean(X,DIM)：按数组X的第DIM维的方向的元素求其平均值赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。,例： x=4 5 6;1 4 8; y1= mean(x,1) y1 =2.5000 4.5
18、000 7.0000 y2= mean(x,2) y2 =5.00004.3333,6、求积 命令格式： Y= prod(X)：将prod(X)返回矩阵X各列元素之积赋予行向量Y。若X为向量，则Y为单变量。 Y= prod(X,DIM)：按数组X的第DIM维的方向的元素求其积赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。,例： x=4 5 6;1 4 8; y1= prod(x,1) y1 =4 20 48 y2= prod(x,2) y2 =12032,7、 求累计和、累积积、标准方差与升序排序,累计和函数CUMS
19、UM 累积积CUMPROD 标准方差STD,微积分,求极限syms x; %定义符号变量 f=(x*(exp(sin(x)+1)-2*(exp(tan(x)-1)/sin(x)3; %确定符号表达式 w=limit(f) %求函数的极限,解方程,解方程函数的格式： solve(expr1,expr2,.,exprN,var1,var2,.varN) 或 solve(expr1,expr2,.,exprN) 功能：求解代数方程组expr1,expr2,.,exprN的根，未知数为var1,var2,.varN。 说明：若不指明符号表达式expr1,expr2,.,exprN的值，系统默认为0。如
20、给出一个表达式x2-3*x-8,则系统将按x2-3*x-8=0进行运算；,【例】解代数方程：a*x2-b*x-6=0 syms a b x solve(a*x2-b*x-6) ans = 1/2/a*(b+(b2+24*a)(1/2) 1/2/a*(b-(b2+24*a)(1/2) 即该方程有两个根： x1=1/2/a*(b+(b2+24*a)(1/2)； x2=1/2/a*(b-(b2+24*a)(1/2),例 研究矮壮素使玉米矮化的效果，在抽穗期测定喷矮壮素小区8株、对照区玉米9株，其观察值如下表：,Matlab实现t检验,从理论上判断，喷施矮壮素只可能矮化无效而不可能促进植物长高，因此假
21、设H0：喷施矮壮素的株高与未喷的相同或更高，即H0： 12对HA： 12，即喷施矮壮素的株高较未喷的为矮。显著水平=0.05。 测验计算：,按=7+8=15，查t 表得一尾t0.05=1.753(一尾测验t0.05等于两尾测验的t0.10)，现实得t=-3.05- t0.05=-1.753，故P0.05。 推断：否定H0： 12，接受HA： 12，即认为玉米喷施矮壮素后，其株高显著地矮于对照。,x=160 160 200 160 200 170 150 210; y=170 270 180 250 270 290 270 230 170; h,p,ci,stats=ttest2(x,y,0.0
22、5,left),例 选生长期、发育进度、植株大小和其它方面皆比较一致的两株番茄构成一组，共得7组，每组中一株接种A处理病毒，另一株接种B处理病毒，以研究不同处理方法的纯化的病毒效果。,假设：两种处理对纯化病毒无不同效果，即： H0：d=0 ；对HA：d0。显著水平=0.01。,测验计算：,查附表4， =7-1=6时，t0.01=3.707。实得|t| t0.01,故P0.01。 推断：否定H0：d=0，接受HA：d0，即A、B两法对纯化病毒的效应有极显著差异。,y1=10 13 8 3 20 20 6; y2=25 12 14 15 27 20 18; h,p,ci,stat=ttest(y1
23、,y2,0.05,both),例 以A、B、C、D4种药剂处理水稻种子，其中A为对照，每处理各得4个苗高观察值(cm)，试做方差分析。,水稻药剂处理苗高方差分析表,x=18 20 10 28; 21 24 15 27; 20 26 17 29; 13 22 14 32 p,anovatab,stats=anova1(x),例 采用5种生长素处理豌豆，未处理为对照，待种子发芽后，分别每盆中移植4株，每组6盆，每盆一个处理，试验共有4组24盆，并按组排列于温室中，使同组各盆的环境条件一致。当各盆见第一朵花时记录4株豌豆的总节间数，结果见下表，试作方差分析。,方差分析表,x=60 62 61 60;
24、 65 65 68 65; 63 61 61 60; 64 67 63 61; 62 65 62 64; 61 62 62 65 p,table,stats=anova2(x),例 一些夏季害虫盛发期的早迟和春季温度高低有关。江苏武进连续9年测定3月下旬至4月中旬旬平均温度累积值(x)和水稻一代三化螟盛发期(y，以5月10日为0)的关系，得结果于表1。试计算其直线回归方程。,由观察值计算一级数据,因而有：,从而得到回归方程：,x=35.5 34.1 31.7 40.3 36.8 40.2 31.7 39.2 44.2 y=12 16 9 2 7 3 13 9 -1 p,s=polyfit(x,y,1),相关分析,x=35.5 34.1 31.7 40.3 36.8 40.2 31.7 39.2 44.2 y=12 16 9 2 7 3 13 9 -1 R,P,RLO,RUP=CORRCOEF(x,y),
