相关系数MATLAB提供了corrcoef函数-课件.ppt

1、标准方差与相关系数标准方差与相关系数1求标准方差求标准方差在在MATLAB中，提供了计算数据序列的标准方差的函数中，提供了计算数据序列的标准方差的函数std。对于向量对于向量X，std(X)返回一个标准方差。对于矩阵返回一个标准方差。对于矩阵A，std(A)返回一个行向量，它的各个元素便是矩阵返回一个行向量，它的各个元素便是矩阵A各列或各列或各行的标准方差。各行的标准方差。std函数的一般调用格式为：函数的一般调用格式为：Y=std(A,flag,dim)其中其中dim取取1或或2。当。当dim=1时，求各列元素的标准方差；当时，求各列元素的标准方差；当dim=2时，则求各行元素的标准方差。时

2、，则求各行元素的标准方差。flag取取0或或1，当，当flag=0时，按时，按1所列公式计算标准方差，当所列公式计算标准方差，当flag=1时，按时，按2所列公式计算标准方差。缺省所列公式计算标准方差。缺省flag=0，dim=1。例例6-7 对二维矩阵对二维矩阵x，从不同维方向求出其标准方差。，从不同维方向求出其标准方差。2相关系数相关系数MATLAB提供了提供了corrcoef函数，可以求出数函数，可以求出数据的相关系数矩阵。据的相关系数矩阵。corrcoef函数的调用格函数的调用格式为：式为：corrcoef(X)：返回从矩阵：返回从矩阵X形成的一个相关系形成的一个相关系数矩阵。此相关系

3、数矩阵的大小与矩阵数矩阵。此相关系数矩阵的大小与矩阵X一一样。它把矩阵样。它把矩阵X的每列作为一个变量，然后的每列作为一个变量，然后求它们的相关系数。求它们的相关系数。corrcoef(X,Y)：在这里，：在这里，X,Y是向量，它们与是向量，它们与corrcoef(X,Y)的作用一样。的作用一样。例例 生成满足正态分布的生成满足正态分布的100005随机矩阵，随机矩阵，然后求各列元素的均值和标准方差，再求然后求各列元素的均值和标准方差，再求这这5列随机数据的相关系数矩阵。列随机数据的相关系数矩阵。命令如下：命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=cor

4、rcoef(X)排序排序MATLAB中对向量中对向量X是排序函数是是排序函数是sort(X)，函数返，函数返回一个对回一个对X中的元素按升序排列的新向量。中的元素按升序排列的新向量。sort函数也可以对矩阵函数也可以对矩阵A的各列或各行重新排序，其的各列或各行重新排序，其调用格式为：调用格式为：Y,I=sort(A,dim)其中其中dim指明对指明对A的列还是行进行排序。若的列还是行进行排序。若dim=1，则按列排；若则按列排；若dim=2，则按行排。，则按行排。Y是排序后的矩是排序后的矩阵，而阵，而I记录记录Y中的元素在中的元素在A中位置。中位置。vubx vlb. .,minbAxt sR

5、xxCnT线性优化线性优化x=lp(C,A,b,vlb,vub) 例例 最小值线性优化f(x)=-5x1-4x2-6x3x1-x2+x3203x1+2x2+4x3423x1+2x230(0 x1, 0 x2,0 x3)%First, enter the coefficients:f = -5; -4; -6 ;A = 1 -1 1 3 2 4 3 2 0;b = 20; 42; 30;lb = 0,0,0; % x的最小值 0,0,0ub = inf,inf,inf;%Next, call a linear programming routine:x= lp(f,A,b,lb,ub)%Ente

6、ring xx = 0.0000 15.0000 3.0000例例 线性优化Min -400 x1-1000 x2-300 x3+200 x4 -2x2 + x3 + x4=0 2x1 +3x2 =16 3x1 +4x2 =0; x3=5c=-400,-1000,-300,200; %目标函数系数A=0 -2 1 1; 2 3 0 0; 3 4 0 0; %约束条件系数b=0; 16; 24;xLB=0,0,0,0; % x取值范围的最小值xUB=inf,inf,5,inf; % x取值范围的最大值x0=0,0,0,0; % x取迭代初始值nEq=1; % 约束条件中只有一个 = 号,其余为=

7、x=lp(c,A,b,xLB,xUB,x0,nEq)disp(最优值为: ,num2str(c*x)0 xg)(. .),(mintsRxxfn非线性优化非线性优化x=constr(f ,x0) fminbnd 计算下面函数在区间(0，1)内的最小值。x3exlogxxcosx)x( fx,fval,exitflag,output=fminbnd(x3+cos(x)+x*log(x)/exp(x),0,1) 在0，5上求下函数的最小值1) 3x()x( f3解：先自定义函数：在MATLAB编辑器中建立M文件为：function f = myfun(x)f = (x-3).2 - 1;保存为my

8、fun.m，然后在命令窗口键入命令：x=fminbnd(myfun,0,5)例例 最小值非线性优化Min f(x)=-x1x2x3, -x1-2x2-2x30,x1+2x2+2x372,初值: x = 10; 10; 10 x = 10; 10; 10%第一步：编写M文件 myfun.mfunction f,g=myfun(x)f=-x(1)*x(2)*x(3);g(1)=-x(1)-2*x(2)-2*x(3);g(2)=x(1)+2*x(2)+2*x(3)-72;%第二步：求解%在MATLAB工作窗中键入x0=10,10,10;x=constr(myfun,x0) %即可第一步：编写M文件

9、fxxgh.mfunction F,G=fxxgh(x)F=-x(1)*x(2);G(1)=(x(1)+x(2)*x(3)-120;第二步：求解在MATLAB工作窗中键入x=1,1,1; % x取迭代初始值options(13)=0; % 约束条件中有0个 = 号,其余为=XL=0,0,2; % x取值范围的最小值XU=inf;inf;inf; % x取值范围的最大值x,options=constr(fxxgh,x,options,XL,XU);options(8) %输出最小值x例例 非线性优化Min f(x)=-x1x2 (x1+ x2)x3=0; x3=2;无约束多元函数最小值 )x(

10、fminxx,x,xxn21多元函数最小值的标准形式为其中：x为向量，如使用fmins求其最小值 求222132131xxx10 xx4x2y的最小值点X=fminsearch(2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2, 0,0) 或在MATLAB编辑器中建立函数文件function f=myfun(x)f=2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2;保存为myfun.m，在命令窗口键入X=fminsearch (myfun, 0,0) 或 X=fminsearch(myfun, 0,0)利用函数fminunc求多变量无

11、约束函数最小值 当函数的阶数大于2时，使用fminunc比fminsearch更有效，但当所选函数高度不连续时，使用fminsearch效果较好 求222121xxx2x3)x( f的最小值 fun=3*x(1)2+2*x(1)*x(2)+x(2)2;x0=1 1;x,fval,exitflag,output,grad,hessian=fminunc(fun,x0)或用下面方法：fun=inline(3*x(1)2+2*x(1)*x(2)+x(2)2)x0=1 1x=fminunc(fun,x0)有约束的多元函数最小值 )x( fminx0)x(C0)x(CeqbxAbeqxAequbxlb非

12、线性有约束的多元函数的标准形式为：sub.to 其中：x、b、beq、lb、ub是向量，A、Aeq为矩阵，C(x)、Ceq(x)是返回向量的函数，f(x)为目标函数，f(x)、C(x)、Ceq(x)可以是非线性函数 fmincon 解：约束条件的标准形式为求下面问题在初始点（0，1）处的最优解21212221x5x2xxxx0 x) 1x(22106x3x221min s.t0 x) 1x(2216x3x22121212221x5x2xxxxmin s.t先在MATLAB编辑器中建立非线性约束函数文件：function c, ceq=mycon (x)c=(x(1)-1)2-x(2);ceq=

13、 ; %无等式约束然后，在命令窗口键入如下命令或建立M文件：fun=x(1)2+x(2)2-x(1)*x(2)-2*x(1)-5*x(2); %目标函数x0=0 1;A=-2 3; %线性不等式约束b=6;Aeq= ; %无线性等式约束beq= ;lb= ; %x没有下、上界ub= ;x,fval,exitflag,output,lambda,grad,hessian=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,mycon) 二次规划问题 xfxHx21minbxAbeqxAeqbuxbl二次规划问题（quadratic programming）的标准形式为：sub.to

14、 其中，H、A、Aeq为矩阵，f、b、beq、lb、ub、x为向量 quadprog 21212221x6x2xxxx21)x( fmin2xx212x2x213xx22121x0,x0 xfxHx21)x( f2111H62f21xxx求解下面二次规划问题sub.to 解：则， 在MATLAB中实现如下：H = 1 -1;-1 2 ;f = -2; -6 ;A = 1 1;-1 2; 2 1 ;b = 2; 2;3 ;lb = zeros(2,1) ;x,fval,exitflag,output,lambda = quadprog(H,f,A,b, , ,lb) “半无限”有约束的多元函数最

15、优解 x、b、beq、lb、ub都是向量；A、Aeq是矩阵；C(x)、Ceq(x)、是返回向量的函数，f(x)为目标函数；f(x)、C(x)、Ceq(x)是非线性函数；为半无限约束，通常是长度为2的向量 fseminf 先建立非线性约束和半无限约束函数文件，并保存为mycon.m：function C,Ceq,K1,K2,S = mycon(X,S)% 初始化样本间距：if isnan(S(1,1), S = 0.2 0; 0.2 0;end% 产生样本集：w1 = 1:S(1,1):100;w2 = 1:S(2,1):100;% 计算半无限约束：K1 = sin(w1*X(1).*cos(w

16、1*X(2) - 1/1000*(w1-50).2 -sin(w1*X(3)-X(3)-1;K2 = sin(w2*X(2).*cos(w2*X(1) - 1/1000*(w2-50).2 -sin(w2*X(3)-X(3)-1;% 无非线性约束：C = ; Ceq= ;% 绘制半无限约束图形plot(w1,K1,-,w2,K2,:),title(Semi-infinite constraints)然后在MATLAB命令窗口或编辑器中建立M文件：fun = sum(x-0.5).2);x0 = 0.5; 0.2; 0.3; % Starting guessx,fval = fseminf(fu

17、n,x0,2,mycon)x(f , )x(f , )x(f , )x(f , )x(f54321304x40 x48xx2)x(f212221122222x3x)x(f18x3x)x(f213214xx)x(f8xx)x(f215 求下列函数最大值的最小化问题其中：先建立目标函数文件，并保存为myfun.m：function f = myfun(x)f(1)= 2*x(1)2+x(2)2-48*x(1)-40*x(2)+304; f(2)= -x(1)2 - 3*x(2)2;f(3)= x(1) + 3*x(2) -18;f(4)= -x(1)- x(2);f(5)= x(1) + x(2)

18、 - 8;然后，在命令窗口键入命令：x0 = 0.1; 0.1; % 初始值x,fval = fminimax(myfun,x0) 求上述问题的绝对值的最大值最小化问题。目标函数为： | )x(f | , | )x(f | , | )x(f | , | )x(f | , | )x(f | 54321解：先建立目标函数文件（与上例相同）然后，在命令窗口或编辑器中建立M文件：x0 = 0.1; 0.1; % 初始点options = optimset(MinAbsMax,5); % 指定绝对值的最小化x,fval = fminimax(myfun,x0, , , , , , , ,options)

19、多目标规划问题 fgoalattain 在大量的应用领域中，人们经常面临用一个解析函数描述数据(通常是测量值)的任务。对这个问题有两种方法。在插值法里，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。曲线拟合或回归是人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。图1说明了这两种方法。连接数据点的实线描绘了线性内插，虚线是数据的最佳拟合。 曲 线 拟 合 00.20.40.60.81-2024681012xy=f(x)Second Order Curve Fittingx=0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1;y=-.447 1.978

20、3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; n=2; % polynomial orderp=polyfit(x, y, n)polyfit 的输出是一个多项式系数的行向量。其解是y = 9.8108x2 20.1293x0.0317。 为了将曲线拟合解与数据点比较，让我们把二者都绘成图。ezplot(-9.8108*x*x+20.1293*x-0.0317)xi=linspace(0, 1, 100); % x-axis data for plottingz=polyval(p, xi); 为了计算在xi数据点的多项式值，调用MATLAB的函数

21、polyval。plot(x, y, o , x, y, xi, z, : ) 画出了原始数据x和y，用o标出该数据点，在数据点之间，再用直线重画原始数据，并用点 : 线，画出多项式数据xi和z。xlabel( x ), ylabel( y=f(x) ), title( Second Order Curve Fitting )曲 线 拟 合 和 插 值 函 数polyfit(x, y, n)对描述n阶多项式y=f(x)的数据进行最小二乘曲线拟合interp1(x, y, xo)1维线性插值interp1(x, y, xo, spline )1维3次样条插值interp1(x, y, xo, c

22、ubic )1维3次插值interp2(x, y, Z, xi, yi)2维线性插值interp2(x, y, Z, xi, yi, cubic )2维3次插值interp2(x, y, Z, xi, yi, nearest )2维最近邻插值 x1=linspace(0, 2*pi, 60);x2=linspace(0, 2*pi, 6);plot(x1, sin(x1), x2, sin(x2), - ) xlabel( x ), ylabel( sin(x) ) title( Linear Interpolation )若不采用直线连接数据点，我们可采用某些更光滑的曲线来拟合数据点。最常用

23、的方法是用一个3阶多项式，即3次多项式，来对相继数据点之间的各段建模，每个3次多项式的头两个导数与该数据点相一致。这种类型的插值被称为3次样次样条条或简称为样条样条。函数interp1也能执行3次样条插值。 hours=1:12; % index for hour data was recordedtemps=5 8 9 15 25 29 31 30 22 25 27 24; % recorded temperaturesplot(hours, temps, hours, temps, + ) % view temperaturestitle( Temperature )xlabel( Hou

24、r ), ylabel( Degrees Celsius )为了说明一维插值，考虑下列问题，12小时内，一小时测量一次室外温度。数据存储在两个MATLAB变量中。 t=interp1(hours, temps, 9.3, spline ) % estimate temperature at hour=9.3 t=interp1(hours, temps, 4.7, spline ) % estimate temperature at hour=4.7 t=interp1(hours, temps, 3.2 6.5 7.1 11.7, spline ) 意，样条插值得到的结果，与上面所示的线性插

25、值的结果不同。因为插值是一个估计或猜测的过程，其意义在于，应用不同的估计规则导致不同的结果。样条插值是对数据进行平滑,也就是，给定一组数据，使用样条插值在更细的间隔求值。例如，hours=1:12; % index for hour data was recordedh=1:0.1:12; % estimate temperature every 1/10 hourtemps=5 8 9 15 25 29 31 30 22 25 27 24;t=interp1(hours, temps, h) ;plot(hours, temps, - , hours, temps, + , h, t) %

26、plot comparative resultstitle( Springfield Temperature )xlabel( Hour ), ylabel( Degrees Celsius )2 二维数据插值二维数据插值在在MATLAB中，提供了解决二维插值问题的函数中，提供了解决二维插值问题的函数interp2，其调用格式为：，其调用格式为：Z1=interp2(X,Y,Z,X1,Y1,method)其中其中X,Y是两个向量，分别描述两个参数的采样点，是两个向量，分别描述两个参数的采样点，Z是与参数采样点对应的函数值，是与参数采样点对应的函数值，X1,Y1是两个向是两个向量或标量，描述欲插

27、值的点。量或标量，描述欲插值的点。Z1是根据相应的插是根据相应的插值方法得到的插值结果。值方法得到的插值结果。 method的取值与一维插的取值与一维插值函数相同。值函数相同。X,Y,Z也可以是矩阵形式。也可以是矩阵形式。同样，同样，X1,Y1的取值范围不能超出的取值范围不能超出X,Y的给定范围，的给定范围，否则，会给出否则，会给出“NaN”错误。错误。例例 某实验对一根长某实验对一根长10米的钢轨进行热源的温度传播测试。用米的钢轨进行热源的温度传播测试。用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示测量时间表示测量时间0:30:60(秒秒)，用用T表示测试所得各点的温度表示测试所得各点的温度()。试用线性插值求出在。试用线性插值求出在一分钟内每隔一分钟内每隔20秒、钢轨每隔秒、钢轨每隔1米处的温度米处的温度TI。命令如下：命令如下：x=0:2.5:10;h=0:30:60;T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41;xi=0:10;hi=0:20:60;TI=interp2(x,h,T,xi,hi)