2.1 线性最小二乘法
曲线拟合问题的提法
是,已知一组(二维)数据,即平面上的
n
个点
(
x
i
,
y
i
)
,
i
=
1,2,
L
,
n
,
x
i
互不相同,寻求一个函数(曲线)
y
=
f
(
x
)
,使
f
(
x
)
在某种准则下 与所有数据点最为接近,即曲线拟合得最好。
线性最小二乘法
是解决曲线拟合最常用的方法,基本思路是,令
![](https://i-blog.csdnimg.cn/blog_migrate/c3579feb97f03844f304f21f0c946b77.png)
![](https://i-blog.csdnimg.cn/blog_migrate/7e4885c852cd2e8e14f55173ec521557.png)
2.2 最小二乘法的 Matlab 实现
2.2.1
解方程组方法
在上面的记号下,
![](https://i-blog.csdnimg.cn/blog_migrate/bdd5962c403aa5c0de8a723e2419b5ac.png)
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x.^2];
ab=r\y
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
2.2.2
多项式拟合方法
如果取
{
r
1
(
x
),
L
,
r
m
+
1
(
x
)}
=
{1,
x
,
L
,
x
m
}
,即用
m
次多项式拟合给定数据,
Matlab 中有现成的函数
a=polyfit(x0,y0,m)
其中输入参数
x0,y0
为要拟合的数据,
m
为拟合多项式的次数,输出参数
a
为拟合多项
式
y=a
m
x
m
+
…+a
1
x+a
0
系数 a=[
a
m
,
…, a
1
, a
0
]。
多项式在 x 处的值 y 可用下面的函数计算
y=polyval(a,x)。
例 5 某乡镇企业 1990-1996 年的生产利润如表 5
![](https://i-blog.csdnimg.cn/blog_migrate/5e51c7a9d4b1366647ea2d6aca790431.png)
试预测 1997 年和 1998 年的利润。
作已知数据的的散点图,
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
plot(x0,y0,'*')
发现该乡镇企业的年生产利润几乎直线上升。因此,我们可以用
y
=
a
1
x
+
a
0
作为
拟合函数来预测该乡镇企业未来的年利润。编写程序如下:
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
a=polyfit(x0,y0,1)
y97=polyval(a,1997)
y98=polyval(a,1998)
2.3 最小二乘优化
在无约束最优化问题中,有些重要的特殊情形,比如目标函数由若干个函数的平方和构成。这类函数一般可以写成:
![](https://i-blog.csdnimg.cn/blog_migrate/2adea1420803ad8c1d47227b12b99407.png)
lsqlin 函数
Matlab 中的函数为:
x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x.^2];
ab=lsqlin(r,y)
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
lsqcurvefit 函数
解 该问题即解最优化问题:
![](https://i-blog.csdnimg.cn/blog_migrate/3e9f33ffd9627c382b735ab9f17547ce.png)
function f=fun1(x,tdata);
f=x(1)+x(2)*exp(-0.02*x(3)*tdata); %其中 x(1)=a,x(2)=b,x(3)=k
调用函数 lsqcurvefit,编写程序如下:
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
x0=[0.2 0.05 0.05];
x=lsqcurvefit(@fun1,x0,td,cd)
lsqnonlin 函数
function f=fun2(x);
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
f=x(1)+x(2)*exp(-0.02*x(3)*td)-cd;
调用函数 lsqnonlin,编写程序如下:
x0=[0.2 0.05 0.05]; %初始值是任意取的
x=lsqnonlin(@fun2,x0)
lsqnonneg 函数
编写程序如下:
c=[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
d=[0.8587;0.1781;0.0747;0.8405];
x=lsqnonneg(c,d)
曲线拟合的用户图形界面求法
Matlab
工具箱提供了命令
cftool
,该命令给出了一维数据拟合的交互式环境。具体
执行步骤如下:
(
1
)把数据导入到工作空间;
(
2
)运行
cftool
,打开用户图形界面窗口;
(
3
)对数据进行预处理;
(
4
)选择适当的模型进行拟合;
(
5
)生成一些相关的统计量,并进行预测。
可以通过帮助(运行
doc cftool
)熟悉该命令的使用细节。
2.4 曲线拟合与函数逼近
前面讲的曲线拟合是已知一组离散数据{(
x
i ,
y
i),
i =1,L,
n},选择一个较简单的函数
f (
x),如多项式,在一定准则如最小二乘准则下,最接近这些数据。
![](https://i-blog.csdnimg.cn/blog_migrate/fc9113e86fdcb9349241065f5a21ef7c.png)
![](https://i-blog.csdnimg.cn/blog_migrate/5b2190d1952c37a34ca7d51676a4fedd.png)
解
编写程序如下:
syms x
base=[1,x^2,x^4];
y1=base.'*base
y2=cos(x)*base.'
r1=int(y1,-pi/2,pi/2)
r2=int(y2,-pi/2,pi/2)
a=r1\r2
xishu1=double(a)
digits(8),xishu2=vpa(a)
求得 xishu1=0.9996 -0.4964 0.0372,即所求的最佳平方逼近多项式为
![](https://i-blog.csdnimg.cn/blog_migrate/1a97d3c7ebebce06276d2f20dfd4171f.png)
2.5 黄河小浪底调水调沙问题
clc,clear
load data.txt %data.txt 按照原始数据格式把水流量和排沙量排成 4 行,12 列
liu=data([1,3],:);
liu=liu';liu=liu(:);
sha=data([2,4],:);
sha=sha';sha=sha(:);
y=sha.*liu;y=y';
i=1:24;
t=(12*i-4)*3600;
t1=t(1);t2=t(end);
pp=csape(t,y);
xsh=pp.coefs %求得插值多项式的系数矩阵,每一行是一个区间上多项式的系数。
TL=quadl(@(tt)ppval(pp,tt),t1,t2)
也可以利用
3
次
B
样条函数进行插值,求得总的排沙量也为
1.844
×
10
9
t
,,计算
的
Matlab
程序如下:
clc,clear
load data.txt %data.txt 按照原始数据格式把水流量和排沙量排成 4 行,12 列
liu=data([1,3],:);
liu=liu';liu=liu(:);
sha=data([2,4],:);
sha=sha';sha=sha(:);
y=sha.*liu;y=y';
i=1:24;
t=(12*i-4)*3600;
t1=t(1);t2=t(end);
pp=spapi(4,t,y) %三次 B 样条
pp2=fn2fm(pp,'pp') %把 B 样条函数转化为 pp 格式
TL=quadl(@(tt)fnval(pp,tt),t1,t2)
对于问题(
2),研究排沙量与水量的关系,从试验数据可以看出,开始排沙量是随
着水流量的增加而增长,而后是随着水流量的减少而减少。显然,变化规律并非是线性
的关系,为此,把问题分为两部分,从开始水流量增加到最大值
2720m
3
/s
(即增长的过
程)为第一阶段,从水流量的最大值到结束为第二阶段,分别来研究水流量与排沙量的
关系。
![](https://i-blog.csdnimg.cn/blog_migrate/dbb3307b3e1b07d6d9e29741b64e5db2.png)
load data.txt
liu=data([1,3],:);
liu=liu';liu=liu(:);
sha=data([2,4],:);
sha=sha';sha=sha(:);
y=sha.*liu;
subplot(1,2,1), plot(liu(1:11),y(1:11),'*')
subplot(1,2,2), plot(liu(12:24),y(12:24),'*')
从散点图可以看出,第一阶段基本上是线性关系,第二阶段准备依次用二次、三次、
四次曲线来拟合,看哪一个模型的剩余标准差小就选取哪一个模型。最后求得第一阶段
排沙量
y
与水流量
v
之间的预测模型为
clc, clear
load data.txt %data.txt 按照原始数据格式把水流量和排沙量排成 4 行,12 列
liu=data([1,3],:); liu=liu'; liu=liu(:);
sha=data([2,4],:); sha=sha'; sha=sha(:);
y=sha.*liu;
%以下是第一阶段的拟合
format long e
nihe1_1=polyfit(liu(1:11),y(1:11),1) %拟合一次多项式,系数排列从高次幂到低次幂
nihe1_2=polyfit(liu(1:11),y(1:11),2)
yhat1_1=polyval(nihe1_1,liu(1:11)); %求预测值
yhat1_2=polyval(nihe1_2,liu(1:11));
%以下求误差平凡和与剩余标准差
cha1_1=sum((y(1:11)-yhat1_1).^2); rmse1_1=sqrt(cha1_1/9)
cha1_2=sum((y(1:11)-yhat1_2).^2); rmse1_2=sqrt(cha1_2/8)
%以下是第二阶段的拟合
for j=1:3
str1=char(['nihe2_' int2str(j) '=polyfit(liu(12:24),y(12:24),' int2str(j+1) ')']);
eval(str1)
str2=char(['yhat2_' int2str(j) '=polyval(nihe2_' int2str(j) ',liu(12:24));']);
eval(str2)
str3=char(['cha2_' int2str(j) '=sum((y(12:24)-yhat2_' int2str(j) ').^2);'...
'rmse2_' int2str(j) '=sqrt(cha2_' int2str(j) '/(11-j))']);
eval(str3)
end
format
更多美赛资料,请点击下方名片获取: