2環(huán)境統(tǒng)計常見數(shù)據(jù)分析方法的MATLAB實現(xiàn)及應(yīng)用

1、環(huán)境統(tǒng)計常見數(shù)據(jù)分析方環(huán)境統(tǒng)計常見數(shù)據(jù)分析方法的法的MATLAB實現(xiàn)及應(yīng)用實現(xiàn)及應(yīng)用第二講第二講一、參數(shù)估計方法一、參數(shù)估計方法線性回歸線性回歸非線性回歸非線性回歸網(wǎng)格搜索網(wǎng)格搜索一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索1、線性回歸、線性回歸MATLAB中調(diào)用函數(shù)：b = regress(y,X) 或或b,bint,r,rint,stats = regress(y,X,alpha)其中b為估計的系數(shù)，bint為b的估計區(qū)間；r為回歸殘差，rint為r的估計區(qū)間；向量stats給出依次給出了R2統(tǒng)計量、F值以及P值；上述參數(shù)是在置信度為1

2、00(1 - alpha)情況下得出的（此時p應(yīng)該小于alpha，模型才成立）。另外，如果回歸模型中沒有考慮常數(shù)項，則上述調(diào)用格式中的X為由np階自變量組的觀測值構(gòu)成的矩陣（每一列表示一個因素），如果回歸模型中包含常數(shù)項，則X為由n(p+1)階矩陣，其第一列全部為1，后面p列由自變量組的觀測值構(gòu)成（每列表示一個因素）。一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索l線性回歸線性回歸-舉例舉例一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索p求解思路求解思路一、參數(shù)估計方法一、參數(shù)估計方法基于線

3、性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索p編程實現(xiàn)編程實現(xiàn)M=10000000;u=0.5;A=20;xx=500;%給出已知條件t=1803004806609001140156018002100240030003600;C=14 1504506246565783933022121476932;y=log(C.*sqrt(t);x1=1./t;x2=t;X=ones(size(t,1),1),x1,x2;%構(gòu)造因變量自變量矩陣b012,bint,r,rint,stats=regress(y,X,0.05) %多元線性回歸T=xx/u;B=b012(3)*(-1) %觀察兩種

4、途徑求得的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一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索2、非線性回歸、非線性回歸上述討論的線性回歸中的上述討論的線性回歸中的“線性線性”并非指并非指y與與x的關(guān)系，而是指的關(guān)系，而是指y是系數(shù)是系數(shù)b0、b1、b2等的線性函數(shù)等的線性函數(shù),在實際科研工作中，在實際科研工作中，y與參數(shù)之間的非線性關(guān)

5、系更為常與參數(shù)之間的非線性關(guān)系更為常見。見。一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索p非線性回歸非線性回歸-MATLAB函數(shù)函數(shù)一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索p非線性回歸非線性回歸-舉例舉例-nlinfit已已知知某某種種化化學(xué)學(xué)物物質(zhì)質(zhì)在在環(huán)環(huán)境境中中的的降降解解速速度度和和多多個個環(huán)環(huán)境境因因素素或或者者條條件件有有關(guān)關(guān)，它它們們之之間間的的定定量量關(guān)關(guān)系系可可用用非非線線性性模模型型表表示示為為：34231253211/xbxbxbbxxbyrate。其其中中 b

6、1b5表表示示要要估估計計的的未未知知參參數(shù)數(shù)；x1x5表表示示環(huán)環(huán)境境因因素素；yrate 表表示示降降解解速速度度。自自變變量量和和因因變變量量觀觀測測值值如如表表。 實實驗驗實實測測環(huán)環(huán)境境因因素素和和反反應(yīng)應(yīng)速速度度數(shù)數(shù)值值序號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 300 1208.55470 80 102.7511100 80 1200.056100 1

7、90 1014.3912285 300 1011.32一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索p【求解】定義模型的M函數(shù)，并給出參數(shù)初始值beta0=b10, b20, b30, b40, b50，然后調(diào)用nlinfit()函數(shù)得到估計的參數(shù)beta、回歸殘差r、雅可比矩陣J。 利用以上輸出結(jié)果以及函數(shù)nlparci()得到非線性模型估計參數(shù)的95%置信度下的置信區(qū)間。 調(diào)用nlpredci()函數(shù)得到非線性模型響應(yīng)值（因變量）的置信區(qū)間。pM函數(shù)程序如下：function yrate=c2fun213(b,x)x1=x(:,1);x2

8、=x(:,2);x3=x(:,3);yrate=(b(1)*x2-x3/b(5)./(1+b(2)*x1+b(3)*x2+b(4)*x3);%數(shù)組點運算一、一、MATLAB基本數(shù)學(xué)運算基本數(shù)學(xué)運算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.3500 13.

9、0000 8.5000 0.0500 11.3200;%定義因變量ybeta0=1 0.5 0.2 0.1 2;%給出參數(shù)初始值beta,r,J=nlinfit(X,y,c2fun213,beta0) %調(diào)用函數(shù)求取參數(shù)betaci=nlparci(beta, r, J) %求參數(shù)求參數(shù)95%置信度下的估計區(qū)間置信度下的估計區(qū)間xinput=470 300 10;285 80 10;470 300 120;470 80 120;470 80 10;%給出自變量一些值給出自變量一些值ypred, yci=nlpredci(c2fun213,xinput,beta, r,J) %得到因變量的估計區(qū)

10、間運行結(jié)果：beta =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 二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p非線性回歸-舉例2（自己練習(xí)）二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p非線性回歸-舉例2（自己練習(xí)）【求解】上述

11、解析解含有余誤差函數(shù)，其手工計算一般要通過查表的方法，而MATLAB中提供了余誤差函數(shù)的求解函數(shù)erfc()，可以直接實現(xiàn)其求解。%首先編制描述解析解模型的函數(shù)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);%然后調(diào)用主要函數(shù)，進(jìn)行參數(shù)估算然后調(diào)用主要函數(shù)，進(jìn)行參數(shù)估算。t=60*3 9 14 2124 2935 3744 5056 60;C=0.000.05 6.00 80.01

12、130.95210.31280.20313.59330.27341.11345.43349.00;Dx0=50;%給出參數(shù)初始值Dx=nlinfit(t,C,c3fun39,Dx0);disp(估計出的縱向彌散系數(shù));Dx 二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p非線性回歸-舉例2（自己練習(xí)）一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索、網(wǎng)格搜索-數(shù)學(xué)原理數(shù)學(xué)原理一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索、網(wǎng)格搜索-算法描述算法描述一、參數(shù)估計方法

13、一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索、網(wǎng)格搜索-應(yīng)用舉例應(yīng)用舉例一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索、網(wǎng)格搜索-編程求解編程求解【求解求解】首先根據(jù)參數(shù)估計的基本思想構(gòu)建目標(biāo)函數(shù)。首先根據(jù)參數(shù)估計的基本思想構(gòu)建目標(biāo)函數(shù)。function Zmin=c3fun317(kd,ka)os=8.32;%mg/Ll0=23;%mg/Lo0=8.2;%mg/Lu=4.2*24;%km/dx=0 9293855;%自變量觀測值自變量觀測值DO=8.28.07.36.47.1;%

14、因變量觀測值因變量觀測值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);一、參數(shù)估計方法一、參數(shù)估計方法基于線性回基于線性回歸歸/非線性回歸、網(wǎng)格搜索非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索、網(wǎng)格搜索 然后，根據(jù)網(wǎng)格搜索素算法，編寫循環(huán)進(jìn)行網(wǎng)格搜索j=0; kamin=2; kamax=5; kdmin= 0.1;kdmax=1.5;tka=0.05;tkd=0.01;%取ka,kd的步長分別為0.01和0.005，則總節(jié)點數(shù)可如此計算N=(kamax-kamin)/tka

15、*(kdmax-kdmin)/tkd;n=1:1:N;kka(n)=0; kkd(n)=0;Z(n)=0; %首先對矩陣進(jìn)行占位for i1=kdmin:tkd:kdmax for i2=kamin:tka:kamax Zmin=c3fun317(i1,i2) ;%計算目標(biāo)函數(shù)值 j=j+1; Z(j)=Zmin; %將目標(biāo)值放到Z中 kka(j)=i2;kkd(j)=i1;%將相應(yīng)參數(shù)置于kka,kkd中 endendZmin,ii=min(Z);%返回矩陣Z中的最小值Zmin和對應(yīng)的位置ii%從位置ii處提取出矩陣MI中的參數(shù)k，這就是搜索到的最優(yōu)參數(shù)值ka=kka(ii)kd=kkd(i

16、i)二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p方差分析方差分析 在生產(chǎn)和科研中，不但影響某個事物的因素眾多，而且即使同一個因素在不同的水平下，影響也可能不同。這些因素或同一因素下的不同水平有的影響大，有的影響小。 方差分析是充分利用現(xiàn)有觀測數(shù)據(jù)推斷某個因素或水平的影響是否顯著。 方差分析的基礎(chǔ)是假設(shè)檢驗，這時假設(shè)H0為同一因素的不同水平觀測指標(biāo)相同，或者不同因素的影響觀測指標(biāo)相同。方差分析一般分為單因素方差分析和多因素方差分析。二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析

17、p應(yīng)用舉例二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p編程求解【求解】編制如下的簡單程序，可實現(xiàn)上述問題求解。x=45 5659;42 5263;46 5165;41 5763;46 5867;40 5158;p,tab,stats=anova1(x) %注意是anova1()而不是anoval()二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p應(yīng)用舉例二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p編程求解【求解】該問題共兩個因素，每個因素又有4種水平，每個水平上又有5個重復(fù)。x=23252114;15201717;2

18、6211619;13162420;11222614; 12221923; 23151423; 14172323;P,tab, stats=anova2(x,5) %兩個因素相交的單元內(nèi)有5個重復(fù)二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p應(yīng)用舉例二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析p編程求解【求解】這是一個多因素方差分析問題，可以如下編程解決。X=33 62 37 63 58 75 63 80;group=cat20;cat20;cat20;cat20;cat40;cat40;cat40;cat40, air200;

19、air200;air400;air400;air200;air200;air400;air400, time1;time2;time1;time2;time1;time2;time1;time2;model=2;%調(diào)用方差分析計算時，計算所有2個水平交互作用零假設(shè)的P值sstype=3;%默認(rèn)的平方和計算類型gnames=cat;air;time;%用于表示三個影響因素P,tab,stats=anovan(X, group,model,sstype,gnames)二、顯著性檢驗二、顯著性檢驗基于方差分析基于方差分析三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布三、趨勢面分析法三、趨勢

20、面分析法污染空間分布污染空間分布p趨勢面分析用某種形式的函數(shù)所代表的曲面來逼近環(huán)境要趨勢面分析用某種形式的函數(shù)所代表的曲面來逼近環(huán)境要素的空間分布。素的空間分布。p環(huán)境要素在空間二維平面上的分布可用二元函數(shù)環(huán)境要素在空間二維平面上的分布可用二元函數(shù)u=f(x,y)(趨勢面方程趨勢面方程)近似表示，在空間三維的分布可近似表示，在空間三維的分布可用三元函數(shù)用三元函數(shù)u=f(x,y,z) (趨勢面方程趨勢面方程)近似表示。近似表示。p該函數(shù)從總體上反映環(huán)境要素空間區(qū)域性變化趨勢，稱為該函數(shù)從總體上反映環(huán)境要素空間區(qū)域性變化趨勢，稱為趨勢面部分；環(huán)境要素在空間分布的實測值與該函數(shù)在對趨勢面部分；環(huán)境要

21、素在空間分布的實測值與該函數(shù)在對應(yīng)坐標(biāo)處的對應(yīng)值之差，稱為偏差部分，偏差反映了局部應(yīng)坐標(biāo)處的對應(yīng)值之差，稱為偏差部分，偏差反映了局部的變化。的變化。p例如在地質(zhì)數(shù)據(jù)分析中，用趨勢面方程來表示地質(zhì)特征的總的區(qū)域性例如在地質(zhì)數(shù)據(jù)分析中，用趨勢面方程來表示地質(zhì)特征的總的區(qū)域性變化規(guī)律，可以認(rèn)為這是由大范圍的系統(tǒng)性因素引起的，用偏差部分變化規(guī)律，可以認(rèn)為這是由大范圍的系統(tǒng)性因素引起的，用偏差部分反映局部性的變化特點可以認(rèn)為是局部因素和隨機(jī)因素引起的，如地反映局部性的變化特點可以認(rèn)為是局部因素和隨機(jī)因素引起的，如地質(zhì)現(xiàn)象中的局部異常。質(zhì)現(xiàn)象中的局部異常。p趨勢面函數(shù)主要是多項式趨勢面，因為多項式理論上可

22、以逼近任意連趨勢面函數(shù)主要是多項式趨勢面，因為多項式理論上可以逼近任意連續(xù)函數(shù)，故用多項式能較好地反映連續(xù)變化的分布趨勢。續(xù)函數(shù)，故用多項式能較好地反映連續(xù)變化的分布趨勢。 三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布1、一次趨勢面模型一次趨勢面模型 三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布2、二次趨勢面模型二次趨勢面模型 三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布3、趨勢面擬合程度的檢驗 三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布p舉例三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布三、趨勢面分析法三、趨勢面分析法污染空間分布污

23、染空間分布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=

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

25、 ;plot(X,Y,bp);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

26、=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三、趨勢面分析法三、趨

27、勢面分析法污染空間分布污染空間分布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=

28、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,

29、bp);for ii=1:1:length(X);text(X(ii),Y(ii),num2str(U(ii);end三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布p結(jié)果三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布p結(jié)果三、趨勢面分析法三、趨勢面分析法污染空間分布污染空間分布p結(jié)果四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p聚類分析是對一群不知道類別的觀察對象按照彼聚類分析是對一群不知道類別的觀察對象按照彼此相似程度進(jìn)行分類，達(dá)到此相似程度進(jìn)行分類，達(dá)到“物以類聚物以類聚”的目的。的目的。p聚類分析既可以對

30、樣品進(jìn)行聚類，也可以對變量聚類分析既可以對樣品進(jìn)行聚類，也可以對變量（指標(biāo)）進(jìn)行聚類。（指標(biāo)）進(jìn)行聚類。p從幾何角度講，聚類分析就是根據(jù)某種準(zhǔn)則將空從幾何角度講，聚類分析就是根據(jù)某種準(zhǔn)則將空間中某些比較接近的點歸為一類，而點之間的接間中某些比較接近的點歸為一類，而點之間的接近程度常用相似系數(shù)和距離兩種量來表示。近程度常用相似系數(shù)和距離兩種量來表示。四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p相似系數(shù)四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p距離四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p聚類分析基本過程

31、聚類分析基本過程聚類分析基本思路是：開始先將聚類分析基本思路是：開始先將n n個樣本各自歸為一個樣本各自歸為一類，即類，即n n類類然后取其中最相似者為一新類，此時總類數(shù)變?yōu)槿缓笕∑渲凶钕嗨普邽橐恍骂悾藭r總類數(shù)變?yōu)閚 n-1-1類，再計算新類與其它類，再計算新類與其它n n2 2個類之間的相似性，個類之間的相似性，選擇最相近者并為又一新類，此時總類數(shù)變?yōu)檫x擇最相近者并為又一新類，此時總類數(shù)變?yōu)閚 n-2-2類類依次類推，直到所有變量都?xì)w為一類為止。該聚類依次類推，直到所有變量都?xì)w為一類為止。該聚類過程可用聚類圖譜表示出來，并在合理選擇聚類過程可用聚類圖譜表示出來，并在合理選擇聚類距離或相似系

32、數(shù)后，得到最終的聚類類別。距離或相似系數(shù)后，得到最終的聚類類別。四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p基于MATLAB的聚類分析（1）計算觀測量（樣本）之間的距離：Y = pdist(X)其中X為nm的矩陣，n為樣本數(shù)，m為指標(biāo)（變量）數(shù)；返回的Y為有(n-1)n/2個匹配距離的向量，這些距離按照（1,2）、（1,3）、（1,n）、（2,3）、（2,n）、（n-1,n）的順序排列，Y也稱為相似矩陣?？梢杂胹quareform()將Y轉(zhuǎn)變?yōu)榉骄仃嚕@樣矩陣中(i,j)位置的元素對應(yīng)樣本i和j之間的距離?；蛘撸篩 = pdis

33、t (X,METRIC)其中METRIC為計算距離時采用的方法，euclid表示歐氏距離，seuclid 為標(biāo)準(zhǔn)化歐氏距離， cityblock 表示布洛克距離，mahal 表示馬氏距離，minkowski 為明科夫斯基距離?；蛘撸篩 = pdist (X, minkowski, p)表示使用明科夫斯基距離計算X數(shù)據(jù)矩陣中樣本之間的距離，p表示計算明科夫斯基距離時取冪次。四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類（2）squareform()函數(shù)：Z = squareform(Y)將pdist ()函數(shù)計算得到的Y轉(zhuǎn)變?yōu)榉骄仃嘮，這樣矩陣中(i,j)位置的元素對應(yīng)樣本i和j之間的距

34、離。（3）創(chuàng)建系統(tǒng)聚類樹函數(shù)：Z=linkage(Y)根據(jù)pdist ()函數(shù)計算得到的Y，使用最短距離法快速創(chuàng)建一個系統(tǒng)聚類樹?；蛘撸篫=linkage(Y,method)其中，method為聚類方法，single 最短距離法，complete 最長距離法，average 平均距離法，centroid 重心距離法，ward 平方和遞增法。（4）繪制聚類譜系圖：H=dendrogram(Z)生成由linkage()函數(shù)得到的系統(tǒng)聚類圖Z的冰柱圖。（5）計算Cophenetic相關(guān)系數(shù)的函數(shù)：C = cophenet(Z,Y)返回Cophenetic相關(guān)系數(shù)C用以衡量linkage()函數(shù)得到的

35、Z距離信息和pdist()函數(shù)得到的Y距離信息之間的擬合程度，該值越接近于1表示擬合程度越好。四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類（6）聚類分析函數(shù)：T = cluster(Z, CUTOFF)根據(jù)linkage()函數(shù)得到的Z來創(chuàng)建聚CUTOFF個類別。（7）系統(tǒng)聚類分析函數(shù)：T = clusterdata(X, CUTOFF)根據(jù)數(shù)據(jù)矩陣X創(chuàng)建分類，當(dāng)0 CUTOFF=1，CUTOFF可以解釋為系統(tǒng)聚類樹中分類的最大個數(shù)。四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類p舉例四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類MATLAB程序如下 ：clear a

36、ll; 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.

37、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)四四 、聚類分析法、聚類分析法環(huán)境樣本聚類環(huán)境樣本聚類五、

38、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷五、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷p判別分析也屬于一種數(shù)值分類方法，但與聚類分析有明顯判別分析也屬于一種數(shù)值分類方法，但與聚類分析有明顯的差別。的差別。p聚類分析前并不知道樣本所屬類別的特征，而在判別分析聚類分析前并不知道樣本所屬類別的特征，而在判別分析中用以建立判別函數(shù)的數(shù)據(jù)事先已經(jīng)知道所屬類別。中用以建立判別函數(shù)的數(shù)據(jù)事先已經(jīng)知道所屬類別。p根據(jù)這些已經(jīng)知道類別的數(shù)據(jù)建立判別函數(shù)，然后去判斷根據(jù)這些已經(jīng)知道類別的數(shù)據(jù)建立判別函數(shù)，然后去判斷未知類別的數(shù)據(jù)屬于哪一類，這就是判別分析。未知類別的數(shù)據(jù)屬于哪一類，這

39、就是判別分析。p判別分析在環(huán)境科學(xué)、化學(xué)、地質(zhì)學(xué)、氣象等領(lǐng)域具有廣判別分析在環(huán)境科學(xué)、化學(xué)、地質(zhì)學(xué)、氣象等領(lǐng)域具有廣泛的應(yīng)用，如根據(jù)環(huán)境樣品判別污染類型，根據(jù)巖石成分泛的應(yīng)用，如根據(jù)環(huán)境樣品判別污染類型，根據(jù)巖石成分判別屬于哪一種巖石，根據(jù)化合物特性判別化合物類型，判別屬于哪一種巖石，根據(jù)化合物特性判別化合物類型，根據(jù)氣象信息判斷近日天氣狀況等等。根據(jù)氣象信息判斷近日天氣狀況等等。五、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷五、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷五、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷p舉例五、判別分析法五、判別分析法環(huán)境樣本類別判斷環(huán)境樣本類別判斷%-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