1、MatlabMatlab在统计分析中的应用在统计分析中的应用MATLAB概述与运算基础概述与运算基础nMATLAB语言：广泛应用于工程计算及数值分析领域的新型高级语言n1984年由美国 MathWorks 公司推向市场n已成为国际公认的最优秀的工程应用开发环境n功能强大、简单易学、编程效率高n在欧美各高等院校，MATLAB已成为线性代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等课程的基本教学工具，成为大学生、硕士生以及博士生必须掌握的基本技能。MATLAB：MATrix LABorotory一一. . MATLAB特点特点: :一一. MATLAB. MATLAB特点
2、特点: :1. 数值计算和符号计算功能n数值计算功能包括：矩阵运算、多项式和有理分式运算、数据统计分析、数值积分、优化处理等。n符号计算将得到问题的解析解。2.MATLAB语言n除了命令行的交互式操作以外，还可以程序方式工作。使用MATLAB很容易地实现C或FORTRAN语言的几乎全部功能。3.图形功能n两个层次的图形命令：一种是对图形句柄进行的低级图形命令，另一种是建立在低级图形命令之上的高级图形命令。利用MATLAB的高级图形命令可以轻而易举地绘制二维、三维乃至四维图形，并可进行图形和坐标的标识、视角和光照设计、色彩精细控制等等。4.应用工具箱n包括：基本部分和各种可选的工具箱。n基本部分
3、中有数百个内部函数。n工具箱分为两大类：功能性工具箱和学科性工具箱。a. 功能性工具箱：主要用来扩充符号计算功能、可视建模仿真功能及文字处理功能等。b. 学科性工具箱：专业性较强，如控制系统工具箱、信号处理工具箱、神经网络工具箱、最优化工具箱、金融工具箱等，用户可以直接利用这些工具箱进行相关领域的科学研究。n命令窗口的命令编辑区 用于输入命令和显示计算结果。n例：键入 A=1 2 3;4 5 6;7 8 9 输出 A = 1 2 3 4 5 6 7 8 9二二. MATLAB. MATLAB集成环境集成环境矩阵运算矩阵运算n【例】 求解线性方程组AX=B 1 1.5 2 9 7 0 3.6 0
4、.5 -4 4 其中A= 7 10 -3 22 33 3 7 8.5 21 6 3 8 0 90 -20 3 -4n B= 20 5 16 n 在MATLAB命令窗口输入命令： n 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;n b=3; -4; 20; 5; 16;n x=abn得到的结果是：n x =n 3.5653n -0.9255n -0.2695n 0.1435n 0.0101【例】 求方程 x4+7x3 +9x-20=0的全部根。 n命令窗口输入：n p=1,7,0,9,-20;
5、%建立多项式系数向量n x=roots(p) %求根n得到的结果是：n x =n -7.2254n -0.4286 + 1.5405in -0.4286 - 1.5405in 1.0826二维图形二维图形n一、一、 plot函数函数n函数格式：plot(x,y) 其中x和y为坐标向量n函数功能：以向量x、y为轴，绘制曲线。n【例例】 在区间0X2内，绘制正弦曲线Y=SIN（X）x=0:pi/100:2*pi;y=sin(x);plot(x,y)n【例例】同时绘制正、余弦两条曲线Y1=SIN(X)和Y2=COS(X)x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot
6、(x,y1,x,y2)n（一）线型与颜色（一）线型与颜色n格式：格式：plot(x,y1,cs,.)n其中c表示颜色， s表示线型一、一、 plot函数函数n【例】 用不同线型和颜色重新绘制上例图形 x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,go,x,y2,b-.)其中：参数go和b-.表示图形的颜色和线型。g：绿色 o：线型为圆圈b：蓝色 -.：表示图形线型为点划线n（二）图形标记（二）图形标记ntitle(加图形标题); nxlabel(加X轴标记); nylabel(加Y轴标记); ntext(X,Y,添加文本); 一、一、 plot函数
7、函数三维图形三维图形n一、一、 plot3函数函数函数格式：plot3(x1,y1,z1,c1,x2,y2,z2,c2,) 其中x1,y1,z1表示三维坐标向量；c1,c2表示线形或颜色。n函数功能：以向量x，y，z为坐标，绘制三维曲线。n【例例】 绘制三维螺旋曲线nt=0:pi/50:10*pi;nx=sin(t),y=cos(t);nplot3(x,y,t);ntitle(helix),text(0,0,0,origin);nxlabel(sin(t),ylabel(cos(t),zlabel(t);ngrid;n二、二、mesh函数函数n绘制三维网格图三维网格图。函数格式：mesh(x,
8、y,z,c)n其中：x，y控制X和Y轴坐标 矩阵z是由(x，y)求得Z轴坐标 (x,y,z)组成三维空间的网格点 c用于控制网格点颜色n【例】 绘制三维网格曲面图nx=0:0.15:2*pi;ny=0:0.15:2*pi;nz=sin(y)*cos(x); %矩阵相乘nmesh(x,y,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); % 计算
9、函数值，zz也是21x21的矩阵 mesh(xx, yy, zz); % 画出立体网状图 n三、三、surf函数函数n绘制三维曲面图三维曲面图，各线条之间的补面用颜色填充。surf函数和mesh函数的调用格式一致。n函数格式: surf (x,y,z)n其中x，y控制X和Y轴坐标，矩阵z是由x，y求得的曲面上Z轴坐标。n【例例】 绘制三维曲面图nx=0:0.15:2*pi;ny=0:0.15:2*pi;nz=sin(y)*cos(x); %矩阵相乘nsurf(x,y,z);nxlabel(x-axis),ylabel(y-axis),zlabel(z-axis);ntitle(3-D surf
10、);【例例】剔透玲珑球剔透玲珑球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);hold off %画外球面hidden off %产生透视效果axis off %不显示坐标轴【例例】卫星返回地球的运动轨线示意。R0=1;%以地球半径为一个单位以地球半径为一个单位a=12*R0;b=9*R0;T0=2*pi; %T0是轨道周期是轨道周期T=5*T0;dt=pi/100;t
11、=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,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 %画地球画地球
12、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 %设视角、画运动轨线设视角、画运动轨线 n五、等高线图五、等高线图n函数contour3n【例例】 多峰函数peaks的等高线图nx,y,z=peaks(30);ncontour3(x,y,z,16);nxlabel(x-axis),ylabel(y-axis),zlabel(z-axis);ntitle(contour3 of pe
13、aks)4.6 动画设计动画设计动画功能函数：getframe、moviein和movien【例例】 播放一个不断变化的眼球程序。nm=moviein(20); %建立一个20个列向量组成的矩阵nfor j=1:20n plot(fft(eye(j+10) %绘制出每一幅眼球图并保存到m矩阵中n m(:,j)=getframe;nendnmovie(m,10);%以每秒10幅的速度播放画面【例例 】求下列三阶线性代数方程组的近似解5426255452321321321xxxxxxxxxMATLAB程序为：A=2 -5 4;1 5 -2;-1 2 4；b=5;6;5；x=Ab1323211123
14、45111xxx543321112345111yyyn解法解法1：分别解方程组 (1)Ax=b1；(2)Ay=b2nA=1 -1 1;5 -4 3;2 1 1;nb1=2;-3;1;nb2=3;4;-5;nx=Ab1nx =n -3.8000n 1.4000n 7.2000y=Ab2 -3.6000 -2.2000 4.4000得两个线性代数方程组的解： (1) x1= -3.8， x2= 1.4， x3= 7.2； (2) y1= -3.6， y2=2.2， y3= 4.4n解法解法2：将两个方程组连在一起求解：Az=bnb=2 3;-3 4;1 -5nz=Abnz =n -3.8000 -
15、3.6000n 1.4000 -2.2000n 7.2000 4.4000一、 基本统计处理n1、查取最大值、查取最大值nMAX函数的命令格式：nY,I= max (X)：将max(X)返回矩阵X的各列中的最大元素值及其该元素的位置赋予行向量Y与I；当X为向量时，则Y与I为单变量。nY,I=max(X, ,DIM)：按数组X的第DIM维的方向查取其最大的元素值及其该元素的位置赋予向量Y与I。n【例例】查找下面数列x的最大值。nx=3 5 9 6 1 8 % 产生数列xnx = 3 5 9 6 1 8ny=max(x) % 查出数列x中的最大值赋予yny = 9ny,l=max(x) % 查出数
16、列x中的最大值及其该元素的位置赋予y,lny = 9nl = 3n【例例】分别查找下面34的二维数组x中各列和各行元素中的最大值。nx=1 8 4 2;9 6 2 5;3 6 7 1 % 产生二维数组xnx = 1 8 4 2n 9 6 2 5n 3 6 7 1ny=max(x) % 查出二维数组x中各列元素的最大值产生赋予行向量yny = 9 8 7 5n2、查取最小值、查取最小值 MIN函数n3、求中值、求中值nY=median(X)：将median(X)返回矩阵X各列元素的中值赋予行向量Y。若X为向量，则Y为单变量。nY=median(X,DIM)：按数组X的第DIM维方向的元素求其中值
17、赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。n【例例】分别求下面数列x1与x2的中值。nx1=9 -2 5 7 12; % 奇数个元素ny1=median(x)ny1 =n 7nx2=9 -2 5 6 7 12; % 偶数个元素ny2=median(x)ny2 =n 6.5000n4、求和、求和n命令格式：nY=sum(X)：将sum(X)返回矩阵X各列元素之和赋予行向量Y；若X为向量，则Y为单变量。Y=sum(X,DIM)：按数组X的第DIM维的方向的元素求其和赋予Y。若DIM=1，为按列操作；若DIM=2，为按
18、行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。n例：nx=4 5 6;1 4 8nx =n 4 5 6n 1 4 8ny=sum(x,1)ny =n 5 9 14ny=sum(x,2)ny =n 15n 13n5、求平均值、求平均值nMEAN函数调用的命令格式：nY= mean(X)：将mean (X)返回矩阵X各列元素之的平均值赋予行向量Y。若X为向量，则Y为单变量。nY= mean(X,DIM)：按数组X的第DIM维的方向的元素求其平均值赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。n例：n
19、x=4 5 6;1 4 8;ny1= mean(x,1)ny1 =n 2.5000 4.5000 7.0000ny2= mean(x,2)ny2 =n 5.0000n 4.3333n6、求积、求积n命令格式：nY= prod(X)：将prod(X)返回矩阵X各列元素之积赋予行向量Y。若X为向量，则Y为单变量。nY= prod(X,DIM)：按数组X的第DIM维的方向的元素求其积赋予向量Y。若DIM=1，为按列操作；若DIM=2，为按行操作。若X为二维数组，Y为一个向量；若X为一维数组，则Y为单变量。n例：nx=4 5 6;1 4 8;ny1= prod(x,1)ny1 =n 4 20 48ny
20、2= prod(x,2)ny2 =n 120n 327、 求累计和、累积积、标准方差与升序排序求累计和、累积积、标准方差与升序排序n累计和函数CUMSUMn累积积CUMPROD n标准方差STD微积分n求极限nsyms x; %定义符号变量nf=(x*(exp(sin(x)+1)-2*(exp(tan(x)-1)/sin(x)3; %确定符号表达式nw=limit(f) %求函数的极限xeextgxxx3sin0sin) 1(2) 1(lim解方程n解方程函数的格式：solve(expr1,expr2,.,exprN,var1,var2,.varN) n或 solve(expr1,expr2,
21、.,exprN) n功能：求解代数方程组expr1,expr2,.,exprN的根，未知数为var1,var2,.varN。n说明：若不指明符号表达式expr1,expr2,.,exprN的值，系统默认为0。如给出一个表达式x2-3*x-8,则系统将按x2-3*x-8=0进行运算；n【例】解代数方程：a*x2-b*x-6=0nsyms a b xnsolve(a*x2-b*x-6) nans =n 1/2/a*(b+(b2+24*a)(1/2)n 1/2/a*(b-(b2+24*a)(1/2)n即该方程有两个根： x1=1/2/a*(b+(b2+24*a)(1/2)；x2=1/2/a*(b-(
22、b2+24*a)(1/2) 例 研究矮壮素使玉米矮化的效果，在抽穗期测定喷矮壮素小区8株、对照区玉米9株，其观察值如下表：y y1 1( (喷施矮壮素喷施矮壮素) ) 160160160160200200160160200200170170150150210210y y2 2( (对照对照) )170170270270180180250250270270290290270270230230170170Matlab实现实现t检验检验 从理论上判断，喷施矮壮素只可能矮化无从理论上判断，喷施矮壮素只可能矮化无效而不可能促进植物长高，因此假设效而不可能促进植物长高，因此假设H H0 0：喷施：喷施矮壮
23、素的株高与未喷的相同或更高，即矮壮素的株高与未喷的相同或更高，即H H0 0： H HA A： _2_1)(_2_1yysyyt) 1() 1()()(21_2222_1121212nnyyyySSSSse2212_2_1nsnsseeyy184005 .37873 .233,3 .17621_2_1SSSScmycmy)(688.18)9181(17.147917.1479875 .378718400_2_12cmssyye05. 3688.183 .2333 .176t 按按=7+8=15=7+8=15，查，查t t 表得一尾表得一尾t t0.050.05=1.753(=1.753(一尾测
24、验一尾测验t t0.050.05等于两尾测验的等于两尾测验的t t0.100.10) )，现实得，现实得t t=-3.05=-3.05- - t t0.050.05=-1.753=-1.753，故故P P0.050.05。推断：否定推断：否定H H0 0： 1 12 2，接受，接受H HA A： 1 12 2，即认为玉米喷施矮壮素后，其株，即认为玉米喷施矮壮素后，其株高显著地矮于对照。高显著地矮于对照。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
25、.05,left)n例例 选生长期、发育进度、植株大小和其它方面皆比较一致的两株番茄构成一组，共得7组，每组中一株接种A处理病毒，另一株接种B处理病毒，以研究不同处理方法的纯化的病毒效果。组别组别y y1 1(A(A法法) )y y2 2(B(B法法) )d d1 110102525-15-152 2131312121 13 38 81414-6-64 43 31515-12-125 520202727-7-76 620202020 0 07 76 61818-12-12假设：两种处理对纯化病毒无不同效果，即：假设：两种处理对纯化病毒无不同效果，即： ；对；对测验计算：测验计算： 3 . 77
26、/517/)12(.1)15(_d) 1()(2_nnddsd_ddsdt4286.2277/)51()12(.1)15(2222dSS13. 33323. 23 . 73323. 2674286.227_tsd查附表查附表4 4， y1=10 13 8 3 20 20 6;y2=25 12 14 15 27 20 18;h,p,ci,stat=ttest(y1,y2,0.05,both)例例 以以A A、B B、C C、D4D4种药剂处理水稻种子，种药剂处理水稻种子，其中其中A A为对照，每处理各得为对照，每处理各得4 4个苗高观察值个苗高观察值(cm)(cm)，试做方差分析。，试做方差分析
27、。药剂药剂苗高观察值苗高观察值总和总和T Ti i平均数平均数 A A181821212020131372721818B B202024242626222292922323C C101015151717141456561414D D28282727292932321161162929T T=336=336=21=21iyy变异来源变异来源DFDFSSSSMSMSF F显著显著F F值值药剂处理间药剂处理间3 3504504168.00168.0020.520.56 6F F0.050.05=3.49=3.49药剂处理内药剂处理内( (误差误差) )121298988.178.17F F0.01
28、0.01=5.95=5.95总变异总变异1515602602水稻药剂处理苗高方差分析表水稻药剂处理苗高方差分析表x=18 20 10 28; 21 24 15 27; 20 26 17 29; 13 22 14 32p,anovatab,stats=anova1(x)例例 采用采用5 5种生长素处理豌豆，未处理为对照，待种子种生长素处理豌豆，未处理为对照，待种子发芽后，分别每盆中移植发芽后，分别每盆中移植4 4株，每组株，每组6 6盆，每盆一个处盆，每盆一个处理，试验共有理，试验共有4 4组组2424盆，并按组排列于温室中，使同盆，并按组排列于温室中，使同组各盆的环境条件一致。当各盆见第一朵花
29、时记录组各盆的环境条件一致。当各盆见第一朵花时记录4 4株豌豆的总节间数，结果见下表，试作方差分析。株豌豆的总节间数，结果见下表，试作方差分析。处理（处理（A A）组（组（B B）总和总和T Ti.i.平均平均I IIIIIIIIIIIIVIV未处理未处理(CK)(CK)606062626161606024324360.860.8赤霉素赤霉素656565656868656526326365.865.8动力精动力精636361616161606024524561.361.3吲哚乙酸吲哚乙酸硫酸腺嘌呤硫酸腺嘌呤马来酸马来酸64646262616167676565626263636262626261
30、616464656525525525325325025063.863.863.363.362.562.5总和总和T.T.j j375375382382377377375375T T15091509变异来源变异来源DFDFSSSSMSMSF FF F0.050.05F F0.010.01组间组间3 35.455.451.821.821 1处理间处理间5 565.8765.8713.1713.174.564.562.902.904.564.56误差误差151543.3043.302.892.89总变异总变异2323114.62114.62方差分析表方差分析表x=60 62 61 60; 65 65
31、 68 65; 63 61 61 60; 64 67 63 61; 62 65 62 64; 61 62 62 65p,table,stats=anova2(x)例例 一些夏季害虫盛发期的早迟和春季温度高低有关。江苏武进连续9年测定3月下旬至4月中旬旬平均温度累积值(x)和水稻一代三化螟盛发期(y，以5月10日为0)的关系，得结果于表1。试计算其直线回归方程。x累积温累积温35.534.131.740.336.840.231.739.244.2y盛发期盛发期12169273139-1由观察值计算一级数据7 .3332 .441 .345 .359xn49.125172 .441 .345 .3
32、52222x70) 1(1612y因而有：因而有：0778.379/7 .333/nxx778. 79/70/nyy度旬天).( /0996. 16356.144/0444.159/xSSSPb)(5485.48)0778.370996. 1(7778. 7天xbyaxy0996. 15845.48从而得到回归方程：从而得到回归方程：x=35.5 34.1 31.7 40.3 36.8 40.2 31.7 39.2 44.2y=12 16 9 2 7 3 13 9 -1p,s=polyfit(x,y,1)相关分析x=35.5 34.1 31.7 40.3 36.8 40.2 31.7 39.2 44.2y=12 16 9 2 7 3 13 9 -1R,P,RLO,RUP=CORRCOEF(x,y)
