浅论matlab多变量拟合

首先申明本人是土木专业的，因为有需要要用到matlab中的拟合用途，今天好好学习了一些关于matlab多变量拟合的东西，从网上下载了一些程序，也运行了一下，就举一些实例，附上源程序吧，主要是两个自变量和三个自变量，一个因变量的拟合。

让自己也更清楚，以后用起来也方便。

原理就是给出一个自变量和因变量的矩阵，然后给出一个自己认为的带有未知数的拟合方程，然后付一组初始值，根据matlab返回的初始值和误差在附一组初始值，知道最后的相关系数较大，也就是误差较小时，就能拟合的比较好，写出拟合后的方程了。

1.广义线性回归拟合和源码(两个自变量，一个因变量，非线性拟合）
【例】这里有这样一组数据，涉及三个变量：p，t 和z，要拟合出 z = f(p,t) 的关系式（非线性的）。

z p 0.8 1 1.2
t
60 9.73875 20.75 36.5987
120 13.5725 29.6325 50.93875
180 18.97875 36.59875 80.13875
240 2075125 38.22125 90.925
300 22.055 44.58 104.7725
为了使得回归分析的结果更加直观，我调用regstats函数，编写了一个更为实用的函数：reglm，代码如下（代码中有调用方法和例子）。

首先写一个M文件：
function stats = reglm(y,X,model,varnames)
% 多重线性回归分析或广义线性回归分析
%
% reglm(y,X)，产生线性回归分析的方差分析表和参数估计结果，并以表格形式显示在屏幕上. 参
% 数X是自变量观测值矩阵，它是n行p列的矩阵. y是因变量观测值向量，它是n行1列的列向量.
%
% stats = reglm(y,X)，还返回一个包括了回归分析的所有诊断统计量的结构体变量stats.
%
% stats = reglm(y,X,model)，用可选的model参数来控制回归模型的类型. model是一个字符串，
% 其可用的字符串如下
% 'linear' 带有常数项的线性模型（默认情况）
% 'interaction' 带有常数项、线性项和交叉项的模型
% 'quadratic' 带有常数项、线性项、交叉项和平方项的模型
% 'purequadratic' 带有常数项、线性项和平方项的模型
%
% stats = reglm(y,X,model,varnames)，用可选的varnames参数指定变量标签. varnames
% 可以是字符矩阵或字符串元胞数组，它的每行的字符或每个元胞的字符串是一个变量的标签，它的行
% 数或元胞数应与X的列数相同. 默认情况下，用X1,X2,…作为变量标签.
%
% 例:
% x = [215 250 180 250 180 215 180 215 250 215 215
% 136.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')
%
% ----------------------------------方差分析表
----------------------------------
% 方差来源自由度平方和均方 F
值 p值
% 回
归 5.0000 15.0277 3.0055 7.6122 0.0219
% 残差 5.0000 1.9742 0.3948
% 总计 10.0000 17.0018
%
% 均方根误差(Root MSE) 0.6284 判定系数
(R-Square) 0.8839
% 因变量均值(Dependent Mean) 4.7273 调整的判定系数(Adj
R-Sq) 0.7678
%
% -----------------------------------参数估计
-----------------------------------
% 变量估计值标准误 t值 p 值
% 常数项 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.0.0.0 $ $Date: 2009/12/22 21:41:00 $
if nargin < 2
error('至少需要两个输入参数');
end
p = size(X,2); % X的列数，即变量个数
if nargin < 3 || isempty(model)
model = 'linear'; % model参数的默认值
end
% 生成变量标签varnames
if nargin < 4 || isempty(varnames)
varname1 = strcat({'X'},num2str([1:p]'));
varnames = makevarnames(varname1,model); % 默认的变量标签
else
if ischar(varnames)
varname1 = cellstr(varnames);
elseif iscell(varnames)
varname1 = varnames(:);
else
error('varnames 必须是字符矩阵或字符串元胞数组');
end
if size(varname1,1) ~= p
error('变量标签数与X的列数不一致');
else
varnames = makevarnames(varname1,model); % 指定的变量标签
end
end
ST = regstats(y,X,model); % 调用regstats函数进行线性回归分析，返回结构体变量ST
f = ST.fstat; % F检验相关结果
t = ST.tstat; % t检验相关结果
% 显示方差分析表
fprintf('\n');
fprintf('------------------------------方差分析表
------------------------------');
fprintf('\n');
fprintf('%s%7s%15s%15s%15s%12s','方差来源','自由度','平方和','均方','F值','p值');
fprintf('\n');
fmt = '%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);
fprintf('\n');
fmt = '%s%13.4f%17.4f%17.4f';
fprintf(fmt,'残差',f.dfe,f.sse,f.sse/f.dfe);
fprintf('\n');
fmt = '%s%13.4f%17.4f';
fprintf(fmt,'总计',f.dfe+f.dfr,f.sse+f.ssr);
fprintf('\n');
fprintf('\n');
% 显示判定系数等统计量
fmt = '%22s%15.4f%25s%10.4f';
fprintf(fmt,'均方根误差(Root MSE)',sqrt(ST.mse),'判定系数
(R-Square)',ST.rsquare);
fprintf('\n');
fprintf(fmt,'因变量均值(Dependent Mean)',mean(y),'调整的判定系数(Adj R-Sq)',...
ST.adjrsquare);
fprintf('\n');
fprintf('\n');
% 显示参数估计及t检验相关结果
fprintf('-------------------------------参数估计
-------------------------------');
fprintf('\n');
fprintf('%8s%18s%15s%15s%12s','变量','估计值','标准误','t值','p值'); fprintf('\n');
for i = 1:size(t.beta,1)
if i == 1
fmt = '%8s%20.4f%17.4f%17.4f%12.4f\n';
fprintf(fmt,'常数项',t.beta(i),t.se(i),t.t(i),t.pval(i));
else
fmt = '%10s%20.4f%17.4f%17.4f%12.4f\n';
fprintf(fmt,varnames{i-1},t.beta(i),t.se(i),t.t(i),t.pval(i)); end
end
if nargout == 1
stats = ST; % 返回一个包括了回归分析的所有诊断统计量的结构体变量end
% -----------------子函数-----------------------
function varnames = makevarnames(varname1,model)
% 生成指定模型的变量标签
p = size(varname1,1);
varname2 = [];
for i = 1:p-1
varname2 = [varname2;strcat(varname1(i),'*',varname1(i+1:end))]; end
varname3 = strcat(varname1,'*',varname1);
switch model
case 'linear'
varnames = varname1;
case 'interaction'
varnames = [varname1;varname2];
case 'quadratic'
varnames = [varname1;varname2;varname3];
case 'purequadratic'
varnames = [varname1;varname3];
end
调用自编的reglm函数进行二次拟合，命令如下：
>> z = [9.73875 20.75 36.59875
13.5725 29.6325 50.93875
18.97875 36.59875 80.13875
20.75125 38.22125 90.925
22.055 44.58 104.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*')
拟合结果：
------------------------------------方差分析表
------------------------------------
方差来源自由度平方和均方 F
值 p值
回
归 5.0000 11548.9147 2309.7829 93.4739 0.0000 残差 9.0000 222.3942 24.7105
总计 14.0000 11771.3089
均方根误差(Root MSE) 4.9710 判定系数
(R-Square) 0.9811
因变量均值(Dependent Mean) 41.2168 调整的判定系数(Adj
R-Sq) 0.9706
-----------------------------------参数估计----------------------------------- 变量估计值标准误 t值 p值常数项 242.6188 69.0439 3.5140 0.0066 p -513.7781 137.3777 -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.2628
2.三个自变量，一个因变量
clear,clc
x1=[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 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.23 1.23 1.23 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 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 48.54
49.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.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)*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(:,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.8 8264106646625,934239.742471165,-25.5844409887743,-1301.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:10
b=lsqcurvefit(fx2,b,X,y);
b=nlinfit(X,y,fx2,b);
end
b
m1=mean(x1);m2=mean(x2);m3=mean(x3);
r1=range(x1); r2=range(x2);r3=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),clf
plot3(x1,x2,y,'o')
stem3(x1,x2,y,'filled')
text(x1,x2,y+.04*ry,str,'fontsize',12)
pause(.0001)
hold on
[x11,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)
shading interp
axis tight
pause(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 x2b ya yb])
% alpha(.4)
% shading interp
% axis tight
% pause(.5001)
% end
xlabel('X1'),ylabel('X2'),zlabel('Y')
figure(2),clf
[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,'fontsize',12)
pause(.0001)
hold on
y2=fx(bm,x11,m2,x33);
surf(x11,x33,y2)
axis([x1a x1b x3a x3b ya yb])
alpha(.85)
shading interp
axis tight
pause(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)
% end
xlabel('X1'),ylabel('X3'),zlabel('Y')
figure(3),clf
plot3(x2,x3,y,'o')
stem3(x2,x3,y,'filled')
text(x2,x3,y+.04*ry,str,'fontsize',12)
pause(.0001)
hold on
[x22,x33]=meshgrid(x2a:r2/75:x2b,x3a:r3/75:x3b);
y3=fx(bm,m1,x22,x33);
surf(x22,x33,y3)
axis([x2a x2b x3a x3b ya yb])
alpha(.85)
shading interp
axis tight
pause(5.0001)
%clf
% for l=1:10
% plot3(x2,x3,y,'o')
% stem3(x2,x3,y,'filled')
% text(x2,x3,y+.04*ry,str,'fontsize',12)
% pause(.0001)
% hold on
% m1=x1a+(l-1)*r1/10;
% y3=fx(bm,m1,x22,x33);
% surf(x22,x33,y3)
% axis([x2a x2b x3a x3b ya yb])
% alpha(.4)
% shading interp
% axis tight
% pause(.5001)
% end
xlabel('X2'),ylabel('X3'),zlabel('Y')
disp([' x1 x2 x3 y yh at'])
yhat=fx(b,x1,x2,x3);
[x1,x2,x3,y+ymin,yhat+ymin]
SSy=var(y)*(n-1)
RSS=(y-yhat)'*(y-yhat) Raqaure=(SSy-RSS)/SSy。