使用nlinfit、fminsearch在matlab中实现基于最小二乘法的
非线性参数拟合
（整理自网上资源）
最小二乘法在曲线拟合中比较普遍。拟合的模型主要有 1.直线型 2.多项式型 3.分数函数型 4.指数函数型 5.对数线性型 6.高斯函数型 ......
一般对于LS问题，通常利用反斜杠运算“\\”、fminsearch或优化工具箱提供的极小化函数求解。在Matlab中，曲线拟合工具箱也提供了曲线拟合的图形界面操作。在命令提示符后键入：cftool，即可根据数据，选择适当的拟合模型。
“\\”命令
1.假设要拟合的多项式是：y=a+b*x+c*x^2.首先建立设计矩阵X： X=[ones(size(x)) x x^2]; 执行： para=X\\y
para中包含了三个参数：para(1)=a;para(2)=b;para(3)=c; 这种方法对于系数是线性的模型也适应。 2.假设要拟合：y=a+b*exp(x)+cx*exp(x^2) 设计矩阵X为
X=[ones(size(x)) exp(x) x.*exp(x.^2)]; para=X\\y
3.多重回归(乘积回归)
设要拟合：y=a+b*x+c*t，其中x和t是预测变量，y是响应变量。设计矩阵为 X=[ones(size(x)) x t] %注意x,t大小相等！ para=X\\y
polyfit函数
polyfit函数不需要输入设计矩阵，在参数估计中，polyfit会根据输入的数据生成设计矩阵。
1.假设要拟合的多项式是：y=a+b*x+c*x^2 p=polyfit(x,y,2)
然后可以使用polyval在t处预测： y_hat=polyval(p,t)
polyfit函数可以给出置信区间。
[p S]=polyfit(x,y,2) %S中包含了标准差
[y_fit,delta] = polyval(p,t,S) %按照拟合模型在t处预测
在每个t处的95%CI为：(y_fit-1.96*delta, y_fit+1.96*delta)
2.指数模型也适应
假设要拟合：y = a+b*exp(x)+c*exp(x.?2) p=polyfit(x,log(y),2)
fminsearch函数
fminsearch是优化工具箱的极小化函数。LS问题的基本思想就是残差的平方和(一种范数，由此，LS产生了许多应用)最小，因此可以利用fminsearch函数进行曲线拟合。
假设要拟合：y = a+b*exp(x)+c*exp(x.?2)
首先建立函数，可以通过m文件或函数句柄建立： x=[......]'; y=[......]';
f=@(p,x) p(1)+p(2)*exp(x)+p(3)*exp(x.?2) %注意向量化:p(1)=a;p(2)=b;p(3)=c;
%可以根据需要选择是否优化参数 %opt=options()
p0=ones(3,1);%初值
para=fminsearch(@(p) (y-f(p,x)).^2,p0) %可以输出Hessian矩阵 res=y-f(para,x)%拟合残差
曲线拟合工具箱
提供了很多拟合函数，对大样本场合比较有效！ 非线性拟合nlinfit函数 clear all;
x1=[0.4292 0.4269 0.381 0.4015 0.4117 0.3017]';
x2=[0.00014 0.00059 0.0126 0.0061 0.00425 0.0443]'; x=[x1 x2];
y=[0.517 0.509 0.44 0.466 0.479 0.309]'; f=@(p,x)
2.350176*p(1)*(1-1/p(2))*(1-(1-x(:,1).^(1/p(2))).^p(2)).^2.*(x(:,1).^(-1/p(2))-1).^(-p(2)).*x(:,1).^(-1/p(2)-0.5).*x(:,2); p0=[8 0.5]';
opt=optimset('TolFun',1e-3,'TolX',1e-3);% [p R]=nlinfit(x,y,f,p0,opt)
2.多项式型的一个例子
1900-2000年的总人口情况的曲线拟合
clear all;close all;
?tool提供了可视化的曲线拟合！
t=[1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000]';
y=[75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 281.4220]'; %t太大，以t的幂作为基函数会导致设计矩阵尺度太差，列变量几乎线性相依。变换为[-1 1]上 s=(t-1950)/50; %plot(s,y,'ro'); %回归线：y=a+bx
mx=mean(s);my=mean(y); sx=std(s);sy=std(y); r=corr(s,y); b=r*sy/sx; a=my-b*mx; rline=a+b.*s; figure;
subplot(3,2,[1 2])
plot(s,y,'ro',s,rline,'k');% title('多项式拟合');
set(gca,'XTick',s,'XTickLabel',sprintf('%d|',t)); %hold on; n=4;
PreYear=[2010 2015 2020];%预测年份 tPreYear=(PreYear-1950)/50; Y=zeros(length(t),n); res=zeros(size(Y)); delta=zeros(size(Y));
PrePo=zeros(length(PreYear),n); Predelta=zeros(size(PrePo)); for i=1:n
[p S(i)]=polyfit(s,y,i);
[Y(:,i) delta(:,i)]=polyval(p,s,S(i));%拟合的Y
[PrePo(:,i) Predelta(:,i)]=polyval(p,tPreYear,S(i));%预测 res(:,i)=y-Y(:,i);%残差 end
% plot(s,Y); 09a自动添加不同颜色
% legend('data','regression line','1st poly','2nd poly','3rd poly','4th poly',2)
% plot(tPreYear,PrePo,'>'); % hold off
% plot(Y,res,'o');%残差图 r=corr(s,Y).^2 %R^2 %拟合误差估计CI
YearAdd=[t;PreYear']; tYearAdd=[s;tPreYear'];
CFtit={'一阶拟合','二阶拟合','三阶拟合','四阶拟合'};
for col=1:n
subplot(3,2,col+2);
plot(s,y,'ro',s,Y(:,col),'g-');%原始数据和拟合数据 legend('Original','Fitted',2); hold on;
plot(s,Y(:,col)+2*delta(:,col),'r:');?% CI plot(s,Y(:,col)-2*delta(:,col),'r:'); plot(tPreYear,PrePo(:,col),'>');%预测值
plot(tPreYear,PrePo(:,col)+2*Predelta(:,col));%预测95% CI plot(tPreYear,PrePo(:,col)-2*Predelta(:,col)); axis([-1.2 1.8 0 400]);
set(gca,'XTick',tYearAdd,'XTickLabel',sprintf('%d|',YearAdd)); title(CFtit{col}); hold off; end
figure;%残差图 for col=1:n
subplot(2,2,col);
plot(Y(:,i),res(:,i),'o'); end
一个非线性的应用例子(多元情况) 要拟合y=a*x1^n1+b*x2^n2+c*x3^n3
%注：只是作为应用，模型不一定正确！！！ %x2=x3!!!
y=[1080.94 1083.03 1162.80 1155.61 1092.82 1099.26 1161.06 1258.05 1299.03 1298.30 1440.22 1641.30 1672.21 1612.73 1658.64 1752.42 1837.99 2099.29 2675.47 2786.33 2881.07]';
x1=[1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2]'; x2=[1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.250 1.275 1.3 1.325 1.350 1.375 1.4 1.425 1.45 1.475 1.5]'; x3=[1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.250 1.275 1.3 1.325 1.350 1.375 1.4 1.425 1.45 1.475 1.5]'; x=[x1 x2 x3];
f=@(p,x) p(1)*x(:,1).^p(2)+p(3)*x(:,2).^p(4)+p(5)*x(:,3).^p(6); p0=ones(6,1);
p=fminsearch(@(p)sum(y-f(p,x)).^2,p0) res=y-f(p,x);
res2=res.^2 %失败的模型
Matlab 自定义函数 自定义函数的途径：
M文件函数(M file function) 在线函数(Inline Function) 匿名函数(Anonymous Function)
1.M文件函数 范例
function c=myadd(a,b)
%这里可以写函数的使用说明，前面以%开头 %在工作区中，help myadd将显示此处的说明 c=a+b;
%end %非必须的
第一行function告诉Matlab这是一个函数，a,b是输入，c是输出，myadd是函数名。以m文件定义的函数必须保存为函数名的形式，上例中，函数应保存为myadd.m。要使用myadd函数，该函数必须在Matlab的搜索路径中。 调用方式：
在Matlab命令符后输入 a=1;b=2;
c=myadd(a,b)
2.在线函数(Inline Function)
通常作为参数传递给另外一个函数。比如fminsearch，lsqcurvefit等函数需要以函数作为参数。
在线函数从字符串表达式创建函数，例如： f=inline('x.^2','x');
matlab非线性参数拟合估计_很好的参考材料
