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为：*delta, y_fit+*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=[ ]';

x2=[ ]';

x=[x1 x2];

y=[ ]';

f=@(p,x)

*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).*x(:,2);

p0=[8 ]';

opt=optimset('TolFun',1e-3,'TolX',1e-3);%

[p R]=nlinfit(x,y,f,p0,opt)

例子例子例子例子例子例子例子例子例子例子例子例子例子例子例子例子

直线型例子

2.多项式型的一个例子

1900-2000年的总人口情况的曲线拟合

clear all;close all;

%cftool提供了可视化的曲线拟合！

t=[1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 2000]';

y=[ ]';

%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));