2环境统计常见数据分析方法的MATLAB实现及应用

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,环境统计常见数据分析方法的MATLAB实现及应用,第二讲,1,一、参数估计方法,线性回归非线性回归网格搜索,2,一、参数估计方法基于线性回归/非线性回归、网格搜索,1、线性回归,MATLAB中调用函数：,b = regress(y,X),或b,bint,r,rint,stats = regress(y,X,alpha),其中b为估计的系数，bint为b的估计区间；r为回归残差，rint为r的估计区间；向量stats给出依次给出了R,2,统计量、F值以及P值；上述参数是在置信度为100(1 - alpha)情况下得出的（此时p应该小于alpha，模型才成立）。,另外，如果回归模型中没有考虑常数项，则上述调用格式中的,X,为由np阶自变量组的观测值构成的矩阵（每一列表示一个因素），如果回归模型中包含常数项，则,X,为由n(p+1)阶矩阵，其第一列全部为1，后面p列由自变量组的观测值构成（每列表示一个因素）。,3,一、参数估计方法基于线性回归/非线性回归、网格搜索,线性回归,-,举例,4,一、参数估计方法基于线性回归/非线性回归、网格搜索,求解思路,5,一、参数估计方法基于线性回归/非线性回归、网格搜索,编程实现,M=10000000;u=0.5;A=20;xx=500;%给出已知条件,t=1803004806609001140156018002100240030003600;,C=141504506246565783933022121476932;,y=log(C.*sqrt(t);x1=1./t;x2=t;X=ones(size(t,1),1),x1,x2;%构造因变量自变量矩阵,b012,bint,r,rint,stats=regress(y,X,0.05) %多元线性回归,T=xx/u;,B=b012(3)*(-1) %观察两种途径求得的B是否相等？,B=(-1)*b012(2)/T2 %观察两种途径求得的B是否相等？,A0=exp(b012(1)-2*B*T);,disp(由B算Dx，);Dx=u2/(4*B),disp(由A0算Dx，);Dx=(M/(A0*A*sqrt(4*pi)2,6,一、参数估计方法基于线性回归/非线性回归、网格搜索,2、非线性回归,上述讨论的线性回归中的“线性”并非指,y,与,x,的关系，而是指,y,是系数b0、b1、b2等的线性函数,在实际科研工作中，,y,与参数之间的非线性关系更为常见。,7,一、参数估计方法基于线性回归/非线性回归、网格搜索,非线性回归-MATLAB函数,8,一、参数估计方法基于线性回归/非线性回归、网格搜索,非线性回归-举例-nlinfit,实验实测环境因素和反应速度数值,序号,x,1,x,2,x,3,yrate,序号,x,1,x,2,x,3,yrate,1,470 300 10,8.55,7,100 80 65,2.54,2,285 80 10,3.79,8,470 190 65,4.35,3,470 300 120,4.82,9,100 300 54,13,4,470 80 120,0.02,10,100 300 120,8.5,5,470 80 10,2.75,11,100 80 120,0.05,6,100 190 10,14.39,12,285 300 10,11.32,9,一、参数估计方法基于线性回归/非线性回归、网格搜索,【求解】定义模型的M函数，并给出参数初始值,beta,0=,b,10,b,20,b,30,b,40,b,50，然后调用nlinfit()函数得到估计的参数,beta,、回归残差,r,、雅可比矩阵,J,。,利用以上输出结果以及函数nlparci()得到非线性模型估计参数的95%置信度下的置信区间。,调用nlpredci()函数得到非线性模型响应值（因变量）的置信区间。,M函数程序如下：,function yrate=c2fun213(b,x),x1=x(:,1);x2=x(:,2);x3=x(:,3);,yrate=(b(1)*x2-x3/b(5)./(1+b(2)*x1+b(3)*x2+b(4)*x3);%数组点运算,10,一、MATLAB基本数学运算,X=,470 300 10,285 80 10,470 300 120,470 80 120,470 80 10,100 190 10,100 80 65,470 190 65,100 300 54,100 300 120,100 80 120,285 300 10,;%定义自变量x,y=8.5500 3.7900 4.8200 0.0200 2.7500 14.3900 2.5400 4.3500 13.0000 8.5000 0.0500 11.3200;%定义因变量y,beta0=1 0.5 0.2 0.1 2;%给出参数初始值,beta,r,J=nlinfit(X,y,c2fun213,beta0) %调用函数求取参数,betaci=nlparci(beta, r, J) %,求参数95%置信度下的估计区间,xinput=470 300 10;285 80 10;470 300 120;470 80 120;470 80 10;%,给出自变量一些值,ypred, yci=nlpredci(c2fun213,xinput,beta, r,J) %得到因变量的估计区间,运行结果,：,beta =1.3871 0.0689 0.0455 0.1220 1.0874,betaci=-0.7541 3.5282；-0.0377 0.1755；-0.0318 0.1228；-0.0602 0.3042；-0.6126 2.7873,ypred = 8.4315 3.9904 4.9571 0.0118 2.6603,yci = 0.2459 0.2219 0.1644 0.1667 0.1419,11,二、显著性检验基于方差分析,非线性回归-举例2（自己练习）,12,二、显著性检验基于方差分析,非线性回归-举例2（自己练习）,【求解】上述解析解含有余误差函数，其手工计算一般要通过查表的方法，而MATLAB中提供了余误差函数的求解函数erfc()，可以直接实现其求解。,%首先编制描述解析解模型的函数,function C=c3fun39(Dx,t),c0=350; %mg/L,x=1000;%m,u=0.6;%m/s,C=(c0/2)*(erfc(x-u*t)./(2*sqrt(Dx*t)+exp(u*x/Dx)*erfc(x+u*t)./(2*sqrt(Dx*t);,%,然后调用主要函数，进行参数估算,。,t=60*3 914 2124 2935 3744 5056 60;,C=0.000.05 6.00 80.01130.95210.31280.20313.59330.27341.11345.43349.00;,Dx0=50;%,给出参数初始值,Dx=nlinfit(t,C,c3fun39,Dx0);disp(估计出的纵向弥散系数);Dx,13,二、显著性检验基于方差分析,非线性回归-举例2（自己练习）,14,一、参数估计方法基于线性回归/非线性回归、网格搜索,3、网格搜索-数学原理,15,一、参数估计方法基于线性回归/非线性回归、网格搜索,3、网格搜索-算法描述,16,一、参数估计方法基于线性回归/非线性回归、网格搜索,3、网格搜索-应用举例,17,一、参数估计方法基于线性回归/非线性回归、网格搜索,3、网格搜索-编程求解,【求解】首先根据参数估计的基本思想构建目标函数。,function Zmin=c3fun317(kd,ka),os=8.32;%mg/L,l0=23;%mg/L,o0=8.2;%mg/L,u=4.2*24;%km/d,x=0 9293855;%自变量观测值,DO=8.28.07.36.47.1;%因变量观测值,O=os-(os-o0)*exp(-ka*x./u)+kd*l0./(ka-kd)*(exp(-ka*x./u)-exp(-kd*x./u);,Zmin=sum(O-DO).2);,18,一、参数估计方法基于线性回归/非线性回归、网格搜索,3、网格搜索,然后，根据网格搜索素算法，编写循环进行网格搜索,j=0; kamin=2; kamax=5; kdmin= 0.1;kdmax=1.5;tka=0.05;tkd=0.01;,%取ka,kd的步长分别为0.01和0.005，则总节点数可如此计算,N=(kamax-kamin)/tka*(kdmax-kdmin)/tkd;,n=1:1:N;,kka(n)=0; kkd(n)=0;Z(n)=0; %首先对矩阵进行占位,for i1=kdmin:tkd:kdmax,for i2=kamin:tka:kamax,Zmin=c3fun317(i1,i2) ;%计算目标函数值,j=j+1; Z(j)=Zmin; %将目标值放到Z中,kka(j)=i2;kkd(j)=i1;%将相应参数置于kka,kkd中,end,end,Zmin,ii=min(Z);%返回矩阵Z中的最小值Zmin和对应的位置ii,%从位置ii处提取出矩阵MI中的参数k，这就是搜索到的最优参数值,ka=kka(ii),kd=kkd(ii),19,二、显著性检验,基于方差分析,20,二、显著性检验基于方差分析,方差分析,在生产和科研中，不但影响某个事物的因素众多，而且即使同一个因素在不同的水平下，影响也可能不同。这些因素或同一因素下的不同水平有的影响大，有的影响小。,方差分析是充分利用现有观测数据推断某个因素或水平的影响是否显著。,方差分析的基础是假设检验，这时假设,H,0为同一因素的不同水平观测指标相同，或者不同因素的影响观测指标相同。方差分析一般分为单因素方差分析和多因素方差分析。,21,二、显著性检验基于方差分析,22,二、显著性检验基于方差分析,应用举例,23,二、显著性检验基于方差分析,编程求解,【求解】,编制如下的简单程序，可实现上述问题求解。,x=455659;425263;465165;415763;465867;405158;,p,tab,stats=,anova,1(x) %注意是,anova,1()而不是,anoval,(),24,二、显著性检验基于方差分析,25,二、显著性检验基于方差分析,应用举例,26,二、显著性检验基于方差分析,编程求解,【求解】该问题共两个因素，每个因素又有4种水平，每个水平上又有5个重复。,x=23252114;15201717;26211619;13162420;11222614; 12221923; 23151423; 14172323;,P,tab, stats=anova2(x,5) %两个因素相交的单元内有5个重复,27,二、显著性检验基于方差分析,28,二、显著性检验基于方差分析,应用举例,29,二、显著性检验基于方差分析,编程求解,【求解】这是一个多因素方差分析问题，可以如下编程解决。,X=33 62 37 63 58 75 63 80;,group=cat20;cat20;cat20;cat20;cat40;cat40;cat40;cat40,air200;air200;air400;air400;air200;air200;air400;air400,time1;time2;time1;time2;time1;time2;time1;time2;,model=2;%调用方差分析计算时，计算所有2个水平交互作用零假设的P值,sstype=3;%默认的平方和计算类型,gnames=cat;air;time;%用于表示三个影响因素,P,tab,stats=anovan(X, group,model,sstype,gnames),30,二、显著性检验基于方差分析,31,三、趋势面分析法,污染空间分布,32,三、趋势面分析法,污染空间分布,趋势面分析用某种形式的函数所代表的曲面来逼近环境要素的空间分布。,环境要素在空间二维平面上的分布可用二元函数,u=f(x,y),(趋势面方程)近似表示，在空间三维的分布可用三元函数,u=f(x,y,z),(趋势面方程)近似表示。,该函数从总体上反映环境要素空间区域性变化趋势，称为趋势面部分；环境要素在空间分布的实测值与该函数在对应坐标处的对应值之差，称为偏差部分，偏差反映了局部的变化。,例如在地质数据分析中，用趋势面方程来表示地质特征的总的区域性变化规律，可以认为这是由大范围的系统性因素引起的，用偏差部分反映局部性的变化特点可以认为是局部因素和随机因素引起的，如地质现象中的局部异常。,趋势面函数主要是多项式趋势面，因为多项式理论上可以逼近任意连续函数，故用多项式能较好地反映连续变化的分布趋势。,33,三、趋势面分析法,污染空间分布,1、,一次趋势面模型,34,三、趋势面分析法,污染空间分布,2、,二次趋势面模型,35,三、趋势面分析法,污染空间分布,3、趋势面拟合程度的检验,36,三、趋势面分析法,污染空间分布,举例,37,三、趋势面分析法,污染空间分布,38,三、趋势面分析法,污染空间分布,MATLAB程序。,%-cfun151,clear;clc;,X=2 2 2 4 5 5 7 7 8 10 11 12 12 12 15;,Y=3 10 13 1 8 14 3 6 11 8 13 3 6 10 13;,U=1.9000 2.3000 1.1000 2.6000 2.2000 1.8000 3.5000 3.1000 1.3000 1.2000 1.4000 1.7000 1.8000 1.2000 1.0000;,alpha=0.01;,disp(一次趋势面拟合);,X0=ones(length(X),1);,X1=X;X2=Y;XX1=X0,X1,X2;yy=U;,A1,bint1,r1,rint1,stats1 = regress(yy,XX1,alpha),UU=A1(1)+A1(2)*X1+A1(3)*X2;,R1=sum(UU-mean(U).2)/sum(U-mean(U).2),xxx=1:1:15;yyy=1:1:15;XXX,YYY=meshgrid(xxx,yyy);,UUU=A1(1)+A1(2)*XXX+A1(3)*YYY;,figure(1);c1,h1=contour(XXX,YYY,UUU,8);clabel(c1,h1);title(一次趋势面拟合);,xlabel(X/km);ylabel(y/km);,hold on ;plot(X,Y,bp);for ii=1:1:length(X);text(X(ii),Y(ii),num2str(U(ii);end,39,三、趋势面分析法,污染空间分布,disp(二次趋势面拟合);,X3=X.*X;X4=X.*Y;X5=Y.*Y;XX2=X0,X1,X2,X3,X4,X5;,A2,bint2,r2,rint2,stats2 = regress(yy,XX2,alpha),UU=A2(1)+A2(2)*X1+A2(3)*X2+A2(4)*X3+A2(5)*X4+A2(6)*X5;,R2=sum(UU-mean(U).2)/sum(U-mean(U).2),xxx=1:1:15;yyy=1:1:15;XXX,YYY=meshgrid(xxx,yyy);,UUU=A2(1)+A2(2)*XXX+A2(3)*YYY+A2(4)*XXX.2+A2(5)*XXX.*YYY+A2(6)*YYY.2;,figure(2);c2,h2=contour(XXX,YYY,UUU,8,b-.);clabel(c2,h2);title(二次趋势面拟合);,xlabel(X/km);ylabel(y/km);,hold on ;plot(X,Y,bp);for ii=1:1:length(X);text(X(ii),Y(ii),num2str(U(ii);end,40,三、趋势面分析法,污染空间分布,disp(三次趋势面拟合);,X6=X.*X.*X;X7=X.*X.*Y;X8=X.*Y.*Y;X9=Y.*Y.*Y;,XX3=X0,X1,X2,X3,X4,X5,X6,X7,X8,X9;,A3,bint3,r3,rint3,stats3 = regress(yy,XX3,alpha),UU=A3(1)+A3(2)*X1+A3(3)*X2+A3(4)*X3+A3(5)*X4+A3(6)*X5+A3(7)*X6+A3(8)*X7+A3(9)*X8+A3(10)*X9;,R3=sum(UU-mean(U).2)/sum(U-mean(U).2),xxx=1:1:15;yyy=1:1:15;,XXX,YYY=meshgrid(xxx,yyy);,UUU=A3(1)+A3(2)*XXX+A3(3)*YYY+A3(4)*XXX.*XXX+A3(5)*XXX.*YYY+A3(6)*YYY.*YYY+A3(7).*XXX.3+A3(8).*XXX.*XXX.*YYY+A3(9).*XXX.*YYY.*YYY+A3(10).*YYY.*YYY.*YYY;,figure(3);c3,h3=contour(XXX,YYY,UUU,8,b:);clabel(c3,h3);,title(三次趋势面拟合);xlabel(X/km);ylabel(y/km);,hold on ;plot(X,Y,bp);for ii=1:1:length(X);text(X(ii),Y(ii),num2str(U(ii);end,41,三、趋势面分析法,污染空间分布,结果,42,三、趋势面分析法,污染空间分布,结果,43,三、趋势面分析法,污染空间分布,结果,44,四 、聚类分析法环境样本聚类,45,四 、聚类分析法环境样本聚类,聚类分析是对一群不知道类别的观察对象按照彼此相似程度进行分类，达到“物以类聚”的目的。,聚类分析既可以对样品进行聚类，也可以对变量（指标）进行聚类。,从几何角度讲，聚类分析就是根据某种准则将空间中某些比较接近的点归为一类，而点之间的接近程度常用相似系数和距离两种量来表示。,46,四 、聚类分析法环境样本聚类,47,四 、聚类分析法环境样本聚类,相似系数,48,四 、聚类分析法环境样本聚类,距离,49,四 、聚类分析法环境样本聚类,聚类分析基本过程,聚类分析基本思路是：开始先将,n,个样本各自归为一类，即,n,类,然后取其中最相似者为一新类，此时总类数变为,n,-1类，再计算新类与其它,n,2个类之间的相似性，选择最相近者并为又一新类，此时总类数变为,n,-2类,依次类推，直到所有变量都归为一类为止。该聚类过程可用聚类图谱表示出来，并在合理选择聚类距离或相似系数后，得到最终的聚类类别。,50,四 、聚类分析法环境样本聚类,51,四 、聚类分析法环境样本聚类,基于MATLAB的聚类分析,（1）计算观测量（样本）之间的距离：,Y = pdist(X),其中X为nm的矩阵，n为样本数，m为指标（变量）数；返回的Y为有(n-1)n/2个匹配距离的向量，这些距离按照（1,2）、（1,3）、（1,n）、（2,3）、（2,n）、（n-1,n）的顺序排列，Y也称为相似矩阵。可以用squareform()将Y转变为方矩阵，这样矩阵中(i,j)位置的元素对应样本i和j之间的距离。,或者：,Y = pdist (X,METRIC),其中METRIC为计算距离时采用的方法，euclid表示欧氏距离，seuclid 为标准化欧氏距离， cityblock 表示布洛克距离，mahal 表示马氏距离，minkowski 为明科夫斯基距离。,或者：,Y = pdist (X, minkowski, p),表示使用明科夫斯基距离计算X数据矩阵中样本之间的距离，p表示计算明科夫斯基距离时取幂次。,52,四 、聚类分析法环境样本聚类,（2）squareform()函数：,Z = squareform(Y),将,pdist (,)函数计算得到的Y转变为方矩阵Z，这样矩阵中(i,j)位置的元素对应样本i和j之间的距离。,（3）创建系统聚类树函数：,Z=linkage(Y),根据,pdist (,)函数计算得到的Y，使用最短距离法快速创建一个系统聚类树。,或者：,Z=linkage(Y,method),其中，method为聚类方法，single 最短距离法，complete 最长距离法，average 平均距离法，centroid 重心距离法，ward 平方和递增法。,（4）绘制聚类谱系图：,H=dendrogram(Z),生成由linkage()函数得到的系统聚类图Z的冰柱图。,（5）,计算Cophenetic相关系数的函数：,C = cophenet(Z,Y),返回Cophenetic相关系数C用以衡量linkage()函数得到的Z距离信息和pdist()函数得到的Y距离信息之间的拟合程度，该值越接近于1表示拟合程度越好。,53,四 、聚类分析法环境样本聚类,（6）聚类分析函数：,T = cluster(Z, CUTOFF),根据linkage()函数得到的Z来创建聚CUTOFF个类别。,（7）系统聚类分析函数：,T = clusterdata(X, CUTOFF),根据数据矩阵X创建分类，当0=1，CUTOFF可以解释为系统聚类树中分类的最大个数。,54,四 、聚类分析法环境样本聚类,举例,55,四 、聚类分析法环境样本聚类,MATLAB程序如下 ：,clear all; clc;,X =3.6600 2.5400 2.2100; 3.3400 2.2700 2.1200;,3.2900 5.7100 1.9000; 6.6400 1.3000 1.9000;,3.8900 1.3100 1.5200; 8.6500 1.0700 3.5000;,4.5500 6.1600 4.2500; 4.7500 5.6000 2.7500;,5.8900 1.3900 1.2300; 4.0500 3.4500 2.5100;,12.5300 3.2800 1.4800; 3.0200 1.5800 1.4300;,0.6400 1.1000 1.0400; 3.6600 1.3200 1.1700;,3.1700 2.8000 1.1500; 3.8400 1.0800 1.0100;,3.9600 1.3600 1.0900; 3.4200 1.6800 1.2500;,3.6600 0.8900 1.1000; 1.1800 0.7800 1.2400;,Y=pdist(X);,Z1=squareform(Y);,Z2=linkage(Y);,H=dendrogram(Z2);,T=cluster(Z2,10);,N,M=size(X);,NN=1:1:N;,TT=NN,T,C=cophenet(Z2,Y),56,四 、聚类分析法环境样本聚类,57,五、判别分析法,环境样本类别判断,58,五、判别分析法,环境样本类别判断,判别分析也属于一种数值分类方法，但与聚类分析有明显的差别。,聚类分析前并不知道样本所属类别的特征，而在判别分析中用以建立判别函数的数据事先已经知道所属类别。,根据这些已经知道类别的数据建立判别函数，然后去判断未知类别的数据属于哪一类，这就是判别分析。,判别分析在环境科学、化学、地质学、气象等领域具有广泛的应用，如根据环境样品判别污染类型，根据岩石成分判别属于哪一种岩石，根据化合物特性判别化合物类型，根据气象信息判断近日天气状况等等。,59,五、判别分析法,环境样本类别判断,60,五、判别分析法,环境样本类别判断,61,五、判别分析法,环境样本类别判断,举例,62,五、判别分析法,环境样本类别判断,%-c6fun6_5,clear all;clc;,TRAINING=0.0560 0.0840 0.0310 0.0380 0.0081 0.0220;,0.0400 0.0550 0.1000 0.1100 0.0220 0.0073 ;,0.0500 0.0740 0.0410 0.4800 0.0071 0.0200 ;,0.0450 0.0500 0.1100 0.1000 0.0250 0.0063 ;,0.0380 0.1300 0.0790 0.1700 0.0580 0.0430 ;,0.0300 0.1100 0.0700 0.1600 0.0500 0.0460 ;,0.0340 0.0950 0.0580 0.1600 0.2000 0.0290 ;,0.0300 0.0900 0.0680 0.1800 0.2200 0.0390 ;,0.0840 0.0660 0.0290 0.3200 0.0120 0.0410 ;,0.0850 0.0760 0.0190 0.3000 0.0100 0.0400 ;,0.0640 0.0720 0.1000 0.2100 0.0280 1.3800 ;,0.0540 0.0650 0.0220 0.2800 0.0210 0.0400 ;,0.0450 9.0000 0.0720 0.2000 0.0350 0.0320 ;,0.0480 0.0890 0.0620 0.2600 0.0380 0.0360 ;,0.0690 0.0870 0.0270 0.0500 0.0890 0.0210 ;,SAMPLE=0.0397 0.0667 0.0893 0.1233 0.0823 0.0142 ;,0.0570 0.0983 0.0533 0.2433 0.0353 0.0397 ;,0.0570 0.0983 0.0533 0.2433 0.0353 0.0397 ;,GROUP=1 1 1 1 2 2 1 1 2 2 2 2 2 2 1;,CLASS = classify(SAMPLE,TRAINING,GROUP),63,五、判别分析法,环境样本类别判断,64,思考与练习,65,问题描述,66,问题求解,【求解】首先编制非线性模型的M函数。,function Zmin=c3fun311(k),D=0.1 0.15 0.30.55 0.8 11.2 1.6;%自变量管道直径观测值,pdata=30 55 82153 294452 592 930;%因变量价格观测值,k1=k(1);k2=k(2);k3=k(3); p=k1+k2*D.k3; %因变量价格预测值,Zmin=pdata-p; %目标函数因变量观测值和预测值的差,调用lsqnonlin ()主函数，开始估算,。,k0=1 1 1;k=lsqnonlin (c3fun311,k0);%拟合出参数k,disp(估计出的参数: ); k1=k(1),k2=k(2),k3=k(3),67,要求用不同的估计方法,68,2、,69,用非线性回归,用非线性回归,70,71,简述聚类分析主要应用在哪些方面？,结合讲义实例掌握聚类分析的MATLAB实现函数。,请问判别分析的主要应用在哪些方面？,结合讲义实例掌握判别分析的MATLAB实现函数。,72,