1、環(huán)境統(tǒng)計常見數(shù)據(jù)分析方法的MATLAB實現(xiàn)及運用第二講一、參數(shù)估計方法線性回歸非線性回歸網(wǎng)格搜索一、參數(shù)估計方法基于線性回歸/非線性回歸、網(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ù)是在置信度為100(1 - alpha)情況下得出的此時p應(yīng)該小于alpha，模型才成立。另外，假設(shè)回歸模型中沒有思索常數(shù)項，那么上述調(diào)用格式中的X為由np階

2、自變量組的觀測值構(gòu)成的矩陣每一列表示一個要素，假設(shè)回歸模型中包含常數(shù)項，那么X為由n(p+1)階矩陣，其第一列全部為1，后面p列由自變量組的觀測值構(gòu)成每列表示一個要素。一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索線性回歸-舉例一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索求解思緒一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索編程實現(xiàn)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;%構(gòu)造因變量自變量矩陣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一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索2、非線性回歸上述討論的線

4、性回歸中的“線性并非指y與x的關(guān)系，而是指y是系數(shù)b0、b1、b2等的線性函數(shù),在實踐科研任務(wù)中，y與參數(shù)之間的非線性關(guān)系更為常見。一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索非線性回歸-MATLAB函數(shù)一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索非線性回歸-舉例-nlinfit實驗實測環(huán)境要素和反響速度數(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 3

5、00 1208.55470 80 102.7511100 80 1200.056100 190 1014.3912285 300 1011.32一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索【求解】定義模型的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ū)間。M函數(shù)程序如下：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);%數(shù)組點運算一、MATLAB根本數(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

7、.3500 13.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ù)95%置信度下的估計區(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ū)間運轉(zhuǎn)結(jié)果：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本人練習(xí)二、顯著性檢驗基于方差分析非線性回歸-舉例2本人練習(xí)【求解】上述解析解含有余誤差函數(shù)，其手工計算普通要經(jīng)過查表的方法，而MATLAB中提供了余誤差函數(shù)

9、的求解函數(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ù)，進展參數(shù)估算。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、參數(shù)初始值Dx=nlinfit(t,C,c3fun39,Dx0);disp(估計出的縱向彌散系數(shù));Dx 二、顯著性檢驗基于方差分析非線性回歸-舉例2本人練習(xí)一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索-數(shù)學(xué)原理一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索-算法描畫一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索-運用舉例一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索-編程求解【求解】首先根據(jù)參數(shù)估計的根本思想構(gòu)建目的函數(shù)。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);一、參數(shù)估計方法基于線性回歸/非線性回歸、網(wǎng)格搜索3、網(wǎng)格搜索 然后，根據(jù)網(wǎng)格搜索素算法，編寫循環(huá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

12、，那么總節(jié)點數(shù)可如此計算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) ;%計算目的函數(shù)值 j=j+1; Z(j)=Zmin; %將目的值放到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

13、，這就是搜索到的最優(yōu)參數(shù)值ka=kka(ii)kd=kkd(ii)二、顯著性檢驗基于方差分析二、顯著性檢驗基于方差分析方差分析 在消費和科研中，不但影響某個事物的要素眾多，而且即使同一個要素在不同的程度下，影響也能夠不同。這些要素或同一要素下的不同程度有的影響大，有的影響小。 方差分析是充分利用現(xiàn)有觀測數(shù)據(jù)推斷某個要素或程度的影響能否顯著。 方差分析的根底是假設(shè)檢驗，這時假設(shè)H0為同一要素的不同程度觀測目的一樣，或者不同要素的影響觀測目的一樣。方差分析普通分為單要素方差分析和多要素方差分析。二、顯著性檢驗基于方差分析二、顯著性檢驗基于方差分析運用舉例二、顯著性檢驗基于方差分析編程求解【求解】編

14、制如下的簡單程序，可實現(xiàn)上述問題求解。x=455659;425263;465165;415763;465867;405158;p,tab,stats=anova1(x) %留意是anova1()而不是anoval()二、顯著性檢驗基于方差分析二、顯著性檢驗基于方差分析運用舉例二、顯著性檢驗基于方差分析編程求解【求解】該問題共兩個要素，每個要素又有4種程度，每個程度上又有5個反復(fù)。x=23252114;15201717;26211619;13162420;11222614; 12221923; 23151423; 14172323;P,tab, stats=anova2(x,5) %兩個要素相交

15、的單元內(nèi)有5個反復(fù)二、顯著性檢驗基于方差分析二、顯著性檢驗基于方差分析運用舉例二、顯著性檢驗基于方差分析編程求解【求解】這是一個多要素方差分析問題，可以如下編程處理。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;%調(diào)用方差分析計算時，計算一切2個程度交互作用

16、零假設(shè)的P值sstype=3;%默許的平方和計算類型gnames=cat;air;time;%用于表示三個影響要素P,tab,stats=anovan(X, group,model,sstype,gnames)二、顯著性檢驗基于方差分析三、趨勢面分析法污染空間分布三、趨勢面分析法污染空間分布趨勢面分析用某種方式的函數(shù)所代表的曲面來逼近環(huán)境要素的空間分布。環(huán)境要素在空間二維平面上的分布可用二元函數(shù)u=f(x,y)(趨勢面方程)近似表示，在空間三維的分布可用三元函數(shù)u=f(x,y,z) (趨勢面方程)近似表示。該函數(shù)從總體上反映環(huán)境要素空間區(qū)域性變化趨勢，稱為趨勢面部分；環(huán)境要素在空間分布的實測值

17、與該函數(shù)在對應(yīng)坐標處的對應(yīng)值之差，稱為偏向部分，偏向反映了部分的變化。例如在地質(zhì)數(shù)據(jù)分析中，用趨勢面方程來表示地質(zhì)特征的總的區(qū)域性變化規(guī)律，可以以為這是由大范圍的系統(tǒng)性要素引起的，用偏向部分反映部分性的變化特點可以以為是部分要素和隨機要素引起的，如地質(zhì)景象中的部分異常。趨勢面函數(shù)主要是多項式趨勢面，由于多項式實際上可以逼近恣意延續(xù)函數(shù)，故用多項式能較好地反映延續(xù)變化的分布趨勢。 三、趨勢面分析法污染空間分布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三、趨勢面分析法污染空間分布結(jié)果三、趨勢面分析法污染空間分布結(jié)果三、趨勢面分析法污染空間分布結(jié)果四 、聚類分析法環(huán)境樣本聚類四 、聚類分析法環(huán)境樣本聚類聚類分析是對一群不知道類別的察看對象按照彼此類似程度進展分類，到達“物以類聚的目的。聚類分析既可以對樣品進展聚類，也可以對變量目的進展聚類。從幾何角度講，聚類分析就是根據(jù)某種準那么將空間中某些比較接近的點歸為一類，而點之間的接近程度常用類似系數(shù)和間隔兩種量來表示。四 、聚類分析法環(huán)境樣本聚類四 、聚類分析法環(huán)境樣本聚類類似系數(shù)四 、聚類分析法環(huán)境樣本聚類間隔四 、聚類分析法環(huán)境樣本聚類聚類分析根本過程聚類分

25、析根本思緒是：開場先將n個樣本各自歸為一類，即n類然后取其中最類似者為一新類，此時總類數(shù)變?yōu)閚-1類，再計算新類與其它n2個類之間的類似性，選擇最相近者并為又一新類，此時總類數(shù)變?yōu)閚-2類依次類推，直到一切變量都歸為一類為止。該聚類過程可用聚類圖譜表示出來，并在合理選擇聚類間隔或類似系數(shù)后，得到最終的聚類類別。四 、聚類分析法環(huán)境樣本聚類四 、聚類分析法環(huán)境樣本聚類基于MATLAB的聚類分析1計算觀丈量樣本之間的間隔：Y = pdist(X)其中X為nm的矩陣，n為樣本數(shù)，m為目的變量數(shù)；前往的Y為有(n-1)n/2個匹配間隔的向量，這些間隔按照1,2、1,3、1,n、2,3、2,n、n-1,

26、n的順序陳列，Y也稱為類似矩陣。可以用squareform()將Y轉(zhuǎn)變?yōu)榉骄仃?，這樣矩陣中(i,j)位置的元素對應(yīng)樣本i和j之間的間隔?；蛘撸篩 = pdist (X,METRIC)其中METRIC為計算間隔時采用的方法，euclid表示歐氏間隔，seuclid 為規(guī)范化歐氏間隔， cityblock 表示布洛克間隔，mahal 表示馬氏間隔，minkowski 為明科夫斯基間隔?；蛘撸篩 = pdist (X, minkowski, p)表示運用明科夫斯基間隔計算X數(shù)據(jù)矩陣中樣本之間的間隔，p表示計算明科夫斯基間隔時取冪次。四 、聚類分析法環(huán)境樣本聚類2squareform()函數(shù)：Z =

27、squareform(Y)將pdist ()函數(shù)計算得到的Y轉(zhuǎn)變?yōu)榉骄仃嘮，這樣矩陣中(i,j)位置的元素對應(yīng)樣本i和j之間的間隔。3創(chuàng)建系統(tǒng)聚類樹函數(shù)：Z=linkage(Y)根據(jù)pdist ()函數(shù)計算得到的Y，運用最短間隔法快速創(chuàng)建一個系統(tǒng)聚類樹。或者：Z=linkage(Y,method)其中，method為聚類方法，single 最短間隔法，complete 最長間隔法，average 平均間隔法，centroid 重心間隔法，ward 平方和遞增法。4繪制聚類譜系圖：H=dendrogram(Z)生成由linkage()函數(shù)得到的系統(tǒng)聚類圖Z的冰柱圖。5計算Cophenetic相關(guān)

28、系數(shù)的函數(shù)：C = cophenet(Z,Y)前往Cophenetic相關(guān)系數(shù)C用以衡量linkage()函數(shù)得到的Z間隔信息和pdist()函數(shù)得到的Y間隔信息之間的擬合程度，該值越接近于1表示擬合程度越好。四 、聚類分析法環(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)境樣本聚類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)四 、聚類分析法環(huán)境樣本聚類五

31、、判別分析法環(huán)境樣本類別判別五、判別分析法環(huán)境樣本類別判別判別分析也屬于一種數(shù)值分類方法，但與聚類分析有明顯的差別。聚類分析前并不知道樣本所屬類別的特征，而在判別分析中用以建立判別函數(shù)的數(shù)據(jù)事先曾經(jīng)知道所屬類別。根據(jù)這些曾經(jīng)知道類別的數(shù)據(jù)建立判別函數(shù)，然后去判別未知類別的數(shù)據(jù)屬于哪一類，這就是判別分析。判別分析在環(huán)境科學(xué)、化學(xué)、地質(zhì)學(xué)、氣候等領(lǐng)域具有廣泛的運用，如根據(jù)環(huán)境樣品判別污染類型，根據(jù)巖石成分判別屬于哪一種巖石，根據(jù)化合物特性判別化合物類型，根據(jù)氣候信息判別近日天氣情況等等。五、判別分析法環(huán)境樣本類別判別五、判別分析法環(huán)境樣本類別判別五、判別分析法環(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 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