淺論matlab多變量擬合

1、專業(yè)資料(兩個自變量，一個因變量，非線性擬合) 涉及三個變量：p, t和 乙要擬合出z = f(p,t)的關(guān)系式609.7387520.7536.598712013.572529.632550.9387518018.9787536.5987580.13875240207512538.2212590.92530022.05544.58104.772511.2z p t為了使得回歸分析的結(jié)果更加直觀，我調(diào)用 reglm，代碼如下(代碼中有調(diào)用方法和例子)。 首先寫一個M文件：function stats = reglm(y,X,model,varnames)%多重線性回歸分析或廣義線性回歸分析%r

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

3、這樣一組數(shù)據(jù)， (非線性的)。0.8word完美格式并以表格形式顯示在屏% reglm(y,X)，產(chǎn)生線性回歸分析的方差分析表和參數(shù)估計結(jié)果， 幕上.參% 數(shù)X是自變量觀測值矩陣，它是n行p列的矩陣.y是因變量觀測值向量，它是n行1列的列向量%，還返回一個包括了回歸分析的所有診斷統(tǒng)計量的結(jié)構(gòu)體變量% stats = reglm(y,X) stats.% stats = reglm(y,X,model)個字符串，%，用可選的model參數(shù)來控制回歸模型的類型.model是其可用的字符串如下'li near''in teracti on' 'quadrati

4、c' 'purequadratic'帶有常數(shù)項的線性模型(默認(rèn)情況) 帶有常數(shù)項、 帶有常數(shù)項、 帶有常數(shù)項、線性項和交叉項的模型線性項、交叉項和平方項的模型線性項和平方項的模型，用可選的varnames參數(shù)指定變量標(biāo)簽.它的每行的字符或每個元胞的字符串是一個變量的stats = reglm(y,X,model,var names) varn ames% 可以是字符矩陣或字符串元胞數(shù)組， 標(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 215136

5、.5 136.5 136.5 138.5 139.5 138.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')方差分析表%值%歸%方差來源5.0000殘差總計自由度值15.02775.000010.0000均方根誤差0.8839平方和均方3.00557.61220.02191.974217.00180.3948(Root MSE)0.6284判定系數(shù)(R-Square)%因變量均值(Dependent Me

6、an)R-Sq) 0.7678%4.7273調(diào)整的判定系數(shù)（Adj參數(shù)估計%變量估計值標(biāo)準(zhǔn)誤值%常數(shù)項30.94282011.11170.0154 C%X10.70400.64051.09920.3218%X2-0.848729.1537-0.02910.9779%X1*X2-0.00580.0044-1.31320.2461%X1*X10.00030.00030.83840.4400%X2*X20.00520.10550.04920.9626Copyright 2009 - 2010 xiezhh.$Revisio n: 1.0.0.0 $ $Date: 2009/12/22 21:41:0

7、0 $ if nargin < 2error('至少需要兩個輸入?yún)?shù)end');p = size(X,2); % X 的列數(shù), if nargin < 3 | isempty(model) model = 'li near：% modelend%生成變量標(biāo)簽 varnamesif n argin < 4 | isempty(var names)即變量個數(shù)參數(shù)的默認(rèn)值varnamel = strcat(X ,n um2str(1:p');varn ames = makevar names(var name1,model); %默認(rèn)的變量標(biāo)簽els

8、eif ischar(var names)varn ame1 = cellstr(var names);elseif iscell(var names)varn ame1 = varn ames(:);elseerror('varnames必須是字符矩陣或字符串元胞數(shù)組');endif size(var name1,1) = perror('變量標(biāo)簽數(shù)與X的列數(shù)不一致');elsevarn ames = makevar names(var name1,model); %指定的變量標(biāo)簽endendST = regstats(y,X,model); %調(diào)用regst

9、ats函數(shù)進(jìn)行線性回歸分析，返回結(jié)構(gòu)體變量STf = ST.fstat; % F檢驗相關(guān)結(jié)果t = ST.tstat; % t檢驗相關(guān)結(jié)果%顯示方差分析表fprin tf('n');fprintf(' 方差分析表');fprin tf('n');fprintf('%s%7s%15s%15s%15s%12s','方差來源','自由度','平方和','均方','F值','p');fprin tf('n');fmt = 

10、9;%s%13.4f%17.4f%17.4f%16.4f%12.4f;fprintf(fmt,'回歸',f.dfr,f.ssr,f.ssr/f.dfr,f.f,f.pval);fprin tf('n');fmt = '%s%13.4f%17.4f%17.4f;fprintf(fmt,'殘差',f.dfe,f.sse,f.sse/f.dfe);fprin tf('n');fmt = '%s%13.4f%17.4f;fprintf(fmt,'總計',f.dfe+f.dfr,f.sse+f.ssr);fp

11、rin tf('n');fprin tf('n');%顯示判定系數(shù)等統(tǒng)計量fmt = '%22s%15.4f%25s%10.4f;fprintf(fmt,'均方根誤差(Root MSE)',sqrt(ST.mse),'判定系數(shù)(R-Square)',ST.rsquare);fprin tf('n');fprintf(fmt,'因變量均值(Dependent Mean)',mean(y),'調(diào)整的判定系數(shù)(Adj R-Sq)',.ST.adjrsquare);fprin tf(

12、'n');fprin tf('n');%顯示參數(shù)估計及t檢驗相關(guān)結(jié)果fprintf(' 參數(shù)估計');fprin tf('n');fprintf('%8s%18s%15s%15s%12s','變量','估計值','標(biāo)準(zhǔn)誤','t 值','p 值');fprin tf('n');for i = 1:size(t.beta,1)if i = 1fmt = '%8s%20.4f%17.4f%17.4f%12.4fn：f

13、printf(fmt,'常數(shù)項',t.beta(i),t.se(i),t.t(i),t.pval(i);elsefmt = '%10s%20.4f%17.4f%17.4f%12.4fn'fprin tf(fmt,var namesi-1,t.beta(i),t.se(i),t.t(i),t.pval(i);endendif n argout = 1stats = ST; %返回一個包括了回歸分析的所有診斷統(tǒng)計量的結(jié)構(gòu)體變量end% 子函數(shù)function varn ames = makevar names(var name1,model)%生成指定模型的變量標(biāo)簽

14、p = size(var name1,1);varn ame2 =;for i = 1:p-1varn ame2 = var name2;strcat(var name1(i),'*',var name1(i+1:e nd);endvarn ame3 = strcat(var name1,'*',var namel);switch modelcase 'li near'varn ames = varn amel;case 'in teract ion'varn ames = var name1;var name2;case 

15、9;quadratic'varn ames = var name1;var name2;var name3;case 'purequadratic'varn ames = var name1;var name3;end調(diào)用自編的reglm函數(shù)進(jìn)行二次擬合，命令如下：>> 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.7725;>> p,t = meshgrid(

16、0.8 1 1.2,60:60:300);>> stats = reglm(z(:),p(:), t(:),'quadratic','p','t');>> pn ew,t new = meshgrid(li nspace(0.8,1.2,20),li nspace(60,300,20);>> pp = pn ew(:);>> tt = tn ew(:);>> zhat = on es(400,1) pp tt pp.*tt pp.A2 tt.A2*stats.beta;>>

17、mesh(p new,t new,reshape(zhat,20,20);>> hold on>> plot3(p(:),t(:),z(:),'k*')擬合結(jié)果：方差分析表方差來源自由度平萬和均萬F值p值回歸5.000011548.91472309.782993.47390.0000殘差9.0000222.394224.7105總計14.000011771.3089均方根誤差(Root MSE)4.9710判定系數(shù)(R-Square) 0.9811因變量均值(Dependent Mean)41.2168調(diào)整的判定系數(shù)(AdjR-Sq) 0.9706參數(shù)估

18、計-變量估計值標(biāo)準(zhǔn)誤t值 p 值常數(shù)項242.618869.04393.51400.0066p-513.7781137.3777-3.73990.0046t-0.36370.1212-3.00020.0150p*t0.60220.09266.50100.0001p*p272.262568.06773.99990.0031t*t-0.00030.0002-1.19460.26282.三個自變量，一個因變量clear,clcx1= 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15333.15 3

19、33.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15333.15 333.15 333.15 328.15 330.65 333.15 335.65 338.15 340.65 343.15 333.15 333.15333.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.311.35 1.39 1.43 1.23 1.23 1.23 1.231.23 1.2 1.2 1.2 1.2 1.

20、2 1.26 1.26 1.26 1.26 1.26 1.23 1.23 1.23 1.23 1.23 1.231.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 9080 80 80 80 80 80 80 80 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

21、48.5449.19 49.78 50.42 47.49 48.21 48.9 49.63 50.32 47.8 48.38 48.91 49.47 50.04 50.4950.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;ymi n=min( y);y=y_ymi n; fX=(b,x1,x2,x3)(b(1)+b(2)*x1+b(3)*x2+b(4)*x3+b(5)*x1.A2+b(6)*x2.A2+b(7)*x342 +b(8)*x1.*x2+b(9)*x

22、1.*x3+b(10)*x2.*x3+b(11)*x1.A3+b(12)*x243+b(13)*x3.A3)./(1 +b(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).A2+b (6)*X(:,2).A2+b (7)*X(:,3).A2+b(8)*X(:,1).*X(:,2)+b(9)*X(:,1).*X(:,3)+b(10)*X(:,2).*X(: ,3

23、)+b(11)*X(:,1).A3+b(12)*X(:,2).A3+b(13)*X(:,3).A3)./(1+b(14)*exp(b(15)*X(:,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

24、.90766627356,10.51891749 78167,-642.229950374061,0.00221335659769481,-244987.606559315,0.15540437371958 1,9.28886223888986e-05,-0.0142397533119651,13.4903417277274,0.0213803812532436 ,-0.00141251443766222,0.000377042917999337,-0.0845412180650883;b=bm;for l=1:10b=lsqcurvefit(fx2,b,X,y);b=nli nfit(X,y

25、,fx2,b);endbm1=mea n(x1);m2=mea n(x2);m3=mea n(x3);r1=ra nge(x1); r2=ran ge(x2);r3=ra nge(x3);ry=ra nge(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=len gth(y);str =nu m2str(1: n');figure(1),clfplot3(x1,x2,y,'o')stem3(x1,x2,y,'f

26、illed')text(x1,x2,y+.04*ry,str,'fo ntsize',12)pause(.0001)hold onx11,x22=meshgrid(x1a:r1/75:x1b,x2a:r2/75:x2b); y1=fx(bm,x11,x22,m3);surf(x11,x22,y1)axis(x1a x1b x2a x2b ya yb)alpha(.85)shad ing in terpaxis tightpause(1.0001)%clf% for 1=1:10%plot3(x1,x2,y,'o')% stem3(x1,x2,y,'

27、;filled')%text(x1,x2,y+.04*ry,str,'fo ntsize',12)%pause(.0001)% hold on%m3=x3a+l*r3/10;%y仁 fx(bm,x11,x22,m3);%surf(x11,x22,y1)% axis(x1a x1b x2a x2b ya yb)% alpha(.4)%shadi ng in terp% axis tight%pause(.5001)% endxlabel('X1'),ylabel('X2'),zlabel('Y')figure(2),clf

28、x11,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,'fo ntsize',12) pause(.0001)hold ony2=fx(bm,x11,m2,x33);surf(x11,x33,y2)axis(x1a x1b x3a x3b ya yb)alpha(.85)shad ing in terpaxis 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,'fo ntsize',12)%pause(.0001)% hold on%m2=x2a+(l-1)*r2/10;%y2=fx(bm,x11,m2,x33);%surf