淺析matlab多變量擬合

1、 8/8 首先申明本人是土木專(zhuān)業(yè)的，因?yàn)橛行枰玫絤atlab中的擬合用途，今天好好學(xué)習(xí)了一些關(guān)于matlab多變量擬合的東西，從網(wǎng)上下載了一些程序，也運(yùn)行了一下，就舉一些實(shí)例，附上源程序吧，主要是兩個(gè)自變量和三個(gè)自變量，一個(gè)因變量的擬合。讓自己也更清楚，以后用起來(lái)也方便。 原理就是給出一個(gè)自變量和因變量的矩陣，然后給出一個(gè)自己認(rèn)為的帶有未知數(shù)的擬合方程，然后付一組初始值，根據(jù)matlab返回的初始值和誤差在附一組初始值，知道最后的相關(guān)系數(shù)較大，也就是誤差較小時(shí)，就能擬合的比較好，寫(xiě)出擬合后的方程了。1.廣義線性回歸擬合和源碼(兩個(gè)自變量，一個(gè)因變量，非線性擬合）【例】這里有這樣一組數(shù)據(jù)，涉

2、及三個(gè)變量：p，t 和z，要擬合出 z = f(p,t) 的關(guān)系式（非線性的）。z p 0.8 1 1.2t609.73875 20.75 36.5987120 13.572529.6325 50.93875180 18.9787536.59875 80.138752402075125 38.22125 90.925300 22.05544.58104.7725為了使得回歸分析的結(jié)果更加直觀，我調(diào)用regstats函數(shù)，編寫(xiě)了一個(gè)更為實(shí)用的函數(shù)：reglm，代碼如下（代碼中有調(diào)用方法和例子）。首先寫(xiě)一個(gè)M文件：function stats = reglm(y,X,model,varnames)

3、% 多重線性回歸分析或廣義線性回歸分析% reglm(y,X)，產(chǎn)生線性回歸分析的方差分析表和參數(shù)估計(jì)結(jié)果，并以表格形式顯示在屏幕上. 參% 數(shù)X是自變量觀測(cè)值矩陣，它是n行p列的矩陣. y是因變量觀測(cè)值向量，它是n行1列的列向量.% stats = reglm(y,X)，還返回一個(gè)包括了回歸分析的所有診斷統(tǒng)計(jì)量的結(jié)構(gòu)體變量stats.% stats = reglm(y,X,model)，用可選的model參數(shù)來(lái)控制回歸模型的類(lèi)型. model是一個(gè)字符串，% 其可用的字符串如下% linear 帶有常數(shù)項(xiàng)的線性模型（默認(rèn)情況）% interaction 帶有常數(shù)項(xiàng)、線性項(xiàng)和交叉項(xiàng)的模型% q

4、uadratic 帶有常數(shù)項(xiàng)、線性項(xiàng)、交叉項(xiàng)和平方項(xiàng)的模型% purequadratic 帶有常數(shù)項(xiàng)、線性項(xiàng)和平方項(xiàng)的模型% stats = reglm(y,X,model,varnames)，用可選的varnames參數(shù)指定變量標(biāo)簽. varnames% 可以是字符矩陣或字符串元胞數(shù)組，它的每行的字符或每個(gè)元胞的字符串是一個(gè)變量的標(biāo)簽，它的行% 數(shù)或元胞數(shù)應(yīng)與X的列數(shù)相同. 默認(rèn)情況下，用X1,X2,作為變量標(biāo)簽.% 例:% x = 215 250 180 250 180 215 180 215 250 215 215% 136.5 136.5 136.5 138.5 139.5 138.5

5、 140.5 140.5 140.5 138.5 138.5;% y = 6.2 7.5 4.8 5.1 4.6 4.6 2.8 3.1 4.3 4.9 4.1;% reglm(y,x,quadratic)% 方差分析表% 方差來(lái)源 自由度 平方和 均方 F值 p值% 回歸 5.0000 15.0277 3.0055 7.6122 0.0219% 殘差 5.0000 1.9742 0.3948% 總計(jì) 10.0000 17.0018% 均方根誤差(Root MSE) 0.6284 判定系數(shù)(R-Square) 0.8839% 因變量均值(Dependent Mean) 4.7273 調(diào)整的判定

6、系數(shù)(Adj R-Sq) 0.7678% 參數(shù)估計(jì)% 變量 估計(jì)值 標(biāo)準(zhǔn)誤 t值 p值% 常數(shù)項(xiàng) 30.9428 2011.1117 0.0154 0.9883% X1 0.7040 0.6405 1.0992 0.3218% X2 -0.8487 29.1537 -0.0291 0.9779% X1*X2 -0.0058 0.0044 -1.3132 0.2461% X1*X1 0.0003 0.0003 0.8384 0.4400% X2*X2 0.0052 0.1055 0.0492 0.9626% Copyright 2009 - 2010 xiezhh.% $Revision: 1.

7、0.0.0 $ $Date: 2009/12/22 21:41:00 $if nargin 2 error(至少需要兩個(gè)輸入?yún)?shù));endp = size(X,2); % X的列數(shù)，即變量個(gè)數(shù)if nargin 3 | isempty(model) model = linear; % model參數(shù)的默認(rèn)值end% 生成變量標(biāo)簽varnamesif nargin z = 9.73875 20.75 36.5987513.5725 29.6325 50.9387518.97875 36.59875 80.1387520.75125 38.22125 90.92522.055 44.58 104.

8、7725; p,t = meshgrid(0.8 1 1.2,60:60:300); stats = reglm(z(:),p(:), t(:),quadratic,p,t); pnew,tnew = meshgrid(linspace(0.8,1.2,20),linspace(60,300,20); pp = pnew(:); tt = tnew(:); zhat = ones(400,1) pp tt pp.*tt pp.2 tt.2*stats.beta; mesh(pnew,tnew,reshape(zhat,20,20); hold on plot3(p(:),t(:),z(:),k

9、*)擬合結(jié)果：方差分析表方差來(lái)源 自由度 平方和 均方 F值 p值回歸 5.0000 11548.9147 2309.7829 93.4739 0.0000殘差 9.0000 222.3942 24.7105總計(jì) 14.0000 11771.3089 均方根誤差(Root MSE) 4.9710 判定系數(shù)(R-Square) 0.9811因變量均值(Dependent Mean) 41.2168 調(diào)整的判定系數(shù)(Adj R-Sq) 0.9706參數(shù)估計(jì) 變量 估計(jì)值 標(biāo)準(zhǔn)誤 t值 p值 常數(shù)項(xiàng) 242.6188 69.0439 3.5140 0.0066 p -513.7781 137.377

10、7 -3.7399 0.0046 t -0.3637 0.1212 -3.0002 0.0150 p*t 0.6022 0.0926 6.5010 0.0001 p*p 272.2625 68.0677 3.9999 0.0031 t*t -0.0003 0.0002 -1.1946 0.26282.三個(gè)自變量，一個(gè)因變量clear,clcx1=333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15

11、333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 328.15 330.65 333.15 335.65 338.15 340.65 343.15 333.15 333.15 333.15 323.15 325.65 345.65 348.15;x2=1.19 1.206 1.228 1.23 1.252 1.27 1.277 1.31 1.35 1.39 1.43 1.23 1.23 1.23 1.23 1.23 1.2 1.2 1.2 1.2 1.2 1.26 1.26 1.26 1.26 1.26 1.231.23 1.23

12、1.23 1.23 1.23 1.23 1.15 1.47 1.51 1.23 1.23 1.23 1.23;x3=80 80 80 80 80 80 80 80 80 80 80 77 78 79 80 81 67 68 69 70 71 86 87 88 89 90 80 80 80 80 80 80 8080 80 80 80 80 80 80;y=59.49 55.16 50.18 49.78 45.75 42.96 41.96 37.87 33.96 30.83 28.29 47.92 48.54 49.19 49.78 50.42 47.49 48.21 48.9 49.63 50

13、.32 47.8 48.38 48.91 49.47 50.04 50.49 50.14 49.79 49.45 49.13 48.81 48.5 74.13 26.18 24.39 51.22 50.85 48.21 47.92;X=x1,x2,x3;ymin=min(y);y=y-ymin;fx=(b,x1,x2,x3)(b(1)+b(2)*x1+b(3)*x2+b(4)*x3+b(5)*x1.2+b(6)*x2.2+b(7)*x3.2+b(8)*x1.*x2+b(9)*x1.*x3+b(10)*x2.*x3+b(11)*x1.3+b(12)*x2.3+b(13)*x3.3)./(1+b(

14、14)*exp(b(15)*x1+b(16)*x2+b(17)*x3+b(18)*x1.*x2+b(19)*x1.*x3+b(20)*x2.*x3);fx2=(b,X,y)(b(1)+b(2)*X(:,1)+b(3)*X(:,2)+b(4)*X(:,3)+b(5)*X(:,1).2+b(6)*X(:,2).2+b(7)*X(:,3).2+b(8)*X(:,1).*X(:,2)+b(9)*X(:,1).*X(:,3)+b(10)*X(:,2).*X(:,3)+b(11)*X(:,1).3+b(12)*X(:,2).3+b(13)*X(:,3).3)./(1+b(14)*exp(b(15)*X(:

15、,1)+b(16)*X(:,2)+b(17)*X(:,3)+b(18)*X(:,1).*X(:,2)+b(19)*X(:,1).*X(:,3)+b(20)*X(:,2).*X(:,3);bm=105091.513651451,1328.10332025611,-711027.452435498,-1213.61405762992,-1.88264106646625,934239.742471165,-25.5844409887743,-1301.90766627356,10.5189174978167,-642.229950374061,0.00221335659769481,-244987.

16、606559315,0.155404373719581,9.28886223888986e-05,-0.0142397533119651,13.4903417277274,0.0213803812532436,-0.00141251443766222,0.000377042917999337,-0.0845412180650883;b=bm;for l=1:10 b=lsqcurvefit(fx2,b,X,y); b=nlinfit(X,y,fx2,b);endbm1=mean(x1);m2=mean(x2);m3=mean(x3);r1=range(x1); r2=range(x2);r3=

17、range(x3);ry=range(y);x1a=min(x1);x1b=max(x1);x2a=min(x2);x2b=max(x2);x3a=min(x3);x3b=max(x3);ya=min(y);yb=max(y);n=length(y);str=num2str(1:n);figure(1),clfplot3(x1,x2,y,o)stem3(x1,x2,y,filled)text(x1,x2,y+.04*ry,str,fontsize,12)pause(.0001)hold onx11,x22=meshgrid(x1a:r1/75:x1b,x2a:r2/75:x2b);y1=fx(

18、bm,x11,x22,m3);surf(x11,x22,y1)axis(x1a x1b x2a x2b ya yb)alpha(.85)shading interpaxis tightpause(1.0001)%clf% for l=1:10% plot3(x1,x2,y,o)% stem3(x1,x2,y,filled)% text(x1,x2,y+.04*ry,str,fontsize,12)% pause(.0001)% hold on% m3=x3a+l*r3/10;% y1=fx(bm,x11,x22,m3);% surf(x11,x22,y1)% axis(x1a x1b x2a

19、x2b ya yb)% alpha(.4)% shading interp% axis tight% pause(.5001)% endxlabel(X1),ylabel(X2),zlabel(Y)figure(2),clfx11,x33=meshgrid(x1a:r1/75:x1b,x3a:r3/75:x3b);plot3(x1,x3,y,o)stem3(x1,x3,y,filled)text(x1,x3,y+.04*ry,str,fontsize,12)pause(.0001)hold ony2=fx(bm,x11,m2,x33);surf(x11,x33,y2)axis(x1a x1b

20、x3a x3b ya yb)alpha(.85)shading interpaxis tightpause(5.0001)%clf% for l=1:10% plot3(x1,x3,y,o)% stem3(x1,x3,y,filled)% text(x1,x3,y+.04*ry,str,fontsize,12)% pause(.0001)% hold on% m2=x2a+(l-1)*r2/10;% y2=fx(bm,x11,m2,x33);% surf(x11,x33,y2)% axis(x1a x1b x3a x3b ya yb)% alpha(.4)% shading interp% axis tight% pause(.5001)% endxlabel(X1),ylabel(X3),zlabel(Y)figure(3),c