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

1、环境统计常见数据分析方法的MATLAB实现及运用第二讲一、参数估计方法线性回归非线性回归网格搜索一、参数估计方法基于线性回归/非线性回归、网格搜索1、线性回归MATLAB中调用函数：b = regress(y,X) 或b,bint,r,rint,stats = regress(y,X,alpha)其中b为估计的系数，bint为b的估计区间；r为回归残差，rint为r的估计区间；向量stats给出依次给出了R2统计量、F值以及P值；上述参数是在置信度为100(1 - alpha)情况下得出的此时p应该小于alpha，模型才成立。另外，假设回归模型中没有思索常数项，那么上述调用格式中的X为由np阶

2、自变量组的观测值构成的矩阵每一列表示一个要素，假设回归模型中包含常数项，那么X为由n(p+1)阶矩阵，其第一列全部为1，后面p列由自变量组的观测值构成每列表示一个要素。一、参数估计方法基于线性回归/非线性回归、网格搜索线性回归-举例一、参数估计方法基于线性回归/非线性回归、网格搜索求解思绪一、参数估计方法基于线性回归/非线性回归、网格搜索编程实现M=10000000;u=0.5;A=20;xx=500;%给出知条件t=1803004806609001140156018002100240030003600;C=141504506246565783933022121476932;y=log(C.*

3、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一、参数估计方法基于线性回归/非线性回归、网格搜索2、非线性回归上述讨论的线

4、性回归中的“线性并非指y与x的关系，而是指y是系数b0、b1、b2等的线性函数,在实践科研任务中，y与参数之间的非线性关系更为常见。一、参数估计方法基于线性回归/非线性回归、网格搜索非线性回归-MATLAB函数一、参数估计方法基于线性回归/非线性回归、网格搜索非线性回归-举例-nlinfit实验实测环境要素和反响速度数值序号x1 x2 x3yrate序号x1 x2 x3yrate1470 300 108.557100 80 652.542285 80 103.798470 190 654.353470 300 1204.829100 300 54134470 80 1200.0210100 3

5、00 1208.55470 80 102.7511100 80 1200.056100 190 1014.3912285 300 1011.32一、参数估计方法基于线性回归/非线性回归、网格搜索【求解】定义模型的M函数，并给出参数初始值beta0=b10, b20, b30, b40, b50，然后调用nlinfit()函数得到估计的参数beta、回归残差r、雅可比矩阵J。 利用以上输出结果以及函数nlparci()得到非线性模型估计参数的95%置信度下的置信区间。 调用nlpredci()函数得到非线性模型呼应值因变量的置信区间。M函数程序如下：function yrate=c2fun213

6、(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);%数组点运算一、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;%定义自变量xy=8.5500 3.7900 4.8200 0.0200 2.7500 14.3900 2.5400 4

7、.3500 13.0000 8.5000 0.0500 11.3200;%定义因变量ybeta0=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

8、=1.3871 0.0689 0.0455 0.1220 1.0874betaci=-0.7541 3.5282；-0.0377 0.1755；-0.0318 0.1228；-0.0602 0.3042；-0.6126 2.7873ypred = 8.4315 3.9904 4.9571 0.0118 2.6603yci = 0.2459 0.2219 0.1644 0.1667 0.1419 二、显著性检验基于方差分析非线性回归-举例2本人练习二、显著性检验基于方差分析非线性回归-举例2本人练习【求解】上述解析解含有余误差函数，其手工计算普通要经过查表的方法，而MATLAB中提供了余误差函数

9、的求解函数erfc()，可以直接实现其求解。%首先编制描画解析解模型的函数function C=c3fun39(Dx,t)c0=350; %mg/Lx=1000;%mu=0.6;%m/sC=(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;%给出

10、参数初始值Dx=nlinfit(t,C,c3fun39,Dx0);disp(估计出的纵向弥散系数);Dx 二、显著性检验基于方差分析非线性回归-举例2本人练习一、参数估计方法基于线性回归/非线性回归、网格搜索3、网格搜索-数学原理一、参数估计方法基于线性回归/非线性回归、网格搜索3、网格搜索-算法描画一、参数估计方法基于线性回归/非线性回归、网格搜索3、网格搜索-运用举例一、参数估计方法基于线性回归/非线性回归、网格搜索3、网格搜索-编程求解【求解】首先根据参数估计的根本思想构建目的函数。function Zmin=c3fun317(kd,ka)os=8.32;%mg/Ll0=23;%mg/L

11、o0=8.2;%mg/Lu=4.2*24;%km/dx=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);一、参数估计方法基于线性回归/非线性回归、网格搜索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

12、，那么总节点数可如此计算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中 endendZmin,ii=min(Z);%前往矩阵Z中的最小值Zmin和对应的位置ii%从位置ii处提取出矩阵MI中的参数k

13、，这就是搜索到的最优参数值ka=kka(ii)kd=kkd(ii)二、显著性检验基于方差分析二、显著性检验基于方差分析方差分析 在消费和科研中，不但影响某个事物的要素众多，而且即使同一个要素在不同的程度下，影响也能够不同。这些要素或同一要素下的不同程度有的影响大，有的影响小。 方差分析是充分利用现有观测数据推断某个要素或程度的影响能否显著。 方差分析的根底是假设检验，这时假设H0为同一要素的不同程度观测目的一样，或者不同要素的影响观测目的一样。方差分析普通分为单要素方差分析和多要素方差分析。二、显著性检验基于方差分析二、显著性检验基于方差分析运用举例二、显著性检验基于方差分析编程求解【求解】编

14、制如下的简单程序，可实现上述问题求解。x=455659;425263;465165;415763;465867;405158;p,tab,stats=anova1(x) %留意是anova1()而不是anoval()二、显著性检验基于方差分析二、显著性检验基于方差分析运用举例二、显著性检验基于方差分析编程求解【求解】该问题共两个要素，每个要素又有4种程度，每个程度上又有5个反复。x=23252114;15201717;26211619;13162420;11222614; 12221923; 23151423; 14172323;P,tab, stats=anova2(x,5) %两个要素相交

15、的单元内有5个反复二、显著性检验基于方差分析二、显著性检验基于方差分析运用举例二、显著性检验基于方差分析编程求解【求解】这是一个多要素方差分析问题，可以如下编程处理。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个程度交互作用

16、零假设的P值sstype=3;%默许的平方和计算类型gnames=cat;air;time;%用于表示三个影响要素P,tab,stats=anovan(X, group,model,sstype,gnames)二、显著性检验基于方差分析三、趋势面分析法污染空间分布三、趋势面分析法污染空间分布趋势面分析用某种方式的函数所代表的曲面来逼近环境要素的空间分布。环境要素在空间二维平面上的分布可用二元函数u=f(x,y)(趋势面方程)近似表示，在空间三维的分布可用三元函数u=f(x,y,z) (趋势面方程)近似表示。该函数从总体上反映环境要素空间区域性变化趋势，称为趋势面部分；环境要素在空间分布的实测值

17、与该函数在对应坐标处的对应值之差，称为偏向部分，偏向反映了部分的变化。例如在地质数据分析中，用趋势面方程来表示地质特征的总的区域性变化规律，可以以为这是由大范围的系统性要素引起的，用偏向部分反映部分性的变化特点可以以为是部分要素和随机要素引起的，如地质景象中的部分异常。趋势面函数主要是多项式趋势面，由于多项式实际上可以逼近恣意延续函数，故用多项式能较好地反映延续变化的分布趋势。 三、趋势面分析法污染空间分布1、一次趋势面模型 三、趋势面分析法污染空间分布2、二次趋势面模型 三、趋势面分析法污染空间分布3、趋势面拟合程度的检验 三、趋势面分析法污染空间分布举例三、趋势面分析法污染空间分布三、趋势

18、面分析法污染空间分布MATLAB程序。%-cfun151clear;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,

19、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,b

20、p);for ii=1:1:length(X);text(X(ii),Y(ii),num2str(U(ii);end三、趋势面分析法污染空间分布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

21、=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三、趋势面分析法污染空间分布disp(三次趋势面拟合);X6=X.*X.*

22、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(

23、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

24、),num2str(U(ii);end三、趋势面分析法污染空间分布结果三、趋势面分析法污染空间分布结果三、趋势面分析法污染空间分布结果四 、聚类分析法环境样本聚类四 、聚类分析法环境样本聚类聚类分析是对一群不知道类别的察看对象按照彼此类似程度进展分类，到达“物以类聚的目的。聚类分析既可以对样品进展聚类，也可以对变量目的进展聚类。从几何角度讲，聚类分析就是根据某种准那么将空间中某些比较接近的点归为一类，而点之间的接近程度常用类似系数和间隔两种量来表示。四 、聚类分析法环境样本聚类四 、聚类分析法环境样本聚类类似系数四 、聚类分析法环境样本聚类间隔四 、聚类分析法环境样本聚类聚类分析根本过程聚类分

25、析根本思绪是：开场先将n个样本各自归为一类，即n类然后取其中最类似者为一新类，此时总类数变为n-1类，再计算新类与其它n2个类之间的类似性，选择最相近者并为又一新类，此时总类数变为n-2类依次类推，直到一切变量都归为一类为止。该聚类过程可用聚类图谱表示出来，并在合理选择聚类间隔或类似系数后，得到最终的聚类类别。四 、聚类分析法环境样本聚类四 、聚类分析法环境样本聚类基于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,

26、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表示计算明科夫斯基间隔时取幂次。四 、聚类分析法环境样本聚类2squareform()函数：Z =

27、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相关

28、系数的函数：C = cophenet(Z,Y)前往Cophenetic相关系数C用以衡量linkage()函数得到的Z间隔信息和pdist()函数得到的Y间隔信息之间的拟合程度，该值越接近于1表示拟合程度越好。四 、聚类分析法环境样本聚类6聚类分析函数：T = cluster(Z, CUTOFF)根据linkage()函数得到的Z来创建聚CUTOFF个类别。7系统聚类分析函数：T = clusterdata(X, CUTOFF)根据数据矩阵X创建分类，当0 CUTOFF=1，CUTOFF可以解释为系统聚类树中分类的最大个数。四 、聚类分析法环境样本聚类举例四 、聚类分析法环境样本聚类MATLA

29、B程序如下 ：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.1

30、000 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,TC=cophenet(Z2,Y)四 、聚类分析法环境样本聚类五

31、、判别分析法环境样本类别判别五、判别分析法环境样本类别判别判别分析也属于一种数值分类方法，但与聚类分析有明显的差别。聚类分析前并不知道样本所属类别的特征，而在判别分析中用以建立判别函数的数据事先曾经知道所属类别。根据这些曾经知道类别的数据建立判别函数，然后去判别未知类别的数据属于哪一类，这就是判别分析。判别分析在环境科学、化学、地质学、气候等领域具有广泛的运用，如根据环境样品判别污染类型，根据岩石成分判别属于哪一种岩石，根据化合物特性判别化合物类型，根据气候信息判别近日天气情况等等。五、判别分析法环境样本类别判别五、判别分析法环境样本类别判别五、判别分析法环境样本类别判别举例五、判别分析法环境样本类别判别%-c6fun6_5clear 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