最小二乘拟合实验报告.docx
《最小二乘拟合实验报告.docx》由会员分享,可在线阅读,更多相关《最小二乘拟合实验报告.docx(12页珍藏版)》请在冰豆网上搜索。
最小二乘拟合实验报告
南昌工程学院
《计算方法》实验报告
课程名称计算方法
系院理学院
专业信息与计算科学
班级12级一班
学生姓名魏志辉
学号16
《最小二乘求解》
1引言
在科学实验和生产实践中,经常要从一组实验数据
出发,寻求函数y=f(x)的一个近似表达式y=φ(x),称为经验公式,从几何上来看,这就是一个曲线拟合的问题。
多项式的插值虽然在一定程度上解决了由函数表求函数近似表达式的问题,但用它来解决这里的问题,是有明显的缺陷的。
首先,由实验提供的数据往往有测试误差。
如果要求近似曲线y=φ(x)严格地通过所给的每个数据点
,就会使曲线保留原来的测试误差,因此当个别数据的误差较大的时候,插值的效果是不理想的。
其次,当实验数据较多时,用插值法得到的近似表达式,明显缺乏实用价值。
在实验中,我们常常用最小二乘法来解决这类问题。
定义
为拟合函数在
处的残差。
为了是近似曲线能尽量反映所给数据点的变化趋势,我们要求
尽可能小。
在最小二乘法中,我们选取
,使得偏差平方和最小,即
,这就是最小二乘法的原理。
2实验目的和要求
运用matlab编写.m文件,要求用最小二乘法确定参数。
以下一组数据中x与y之间存在着
的关系,利用最小二乘法确定式中的参数a和b,并计算相应的军方误差与最大偏差。
数据如下:
x
1
2
3
4
5
6
7
8
9
10
y
x
11
12
13
14
15
16
17
18
19
y
3算法原理与流程图
(1)原理
最小二乘是要求对于给定数据列
,要求存在某个函数类
中寻求一个函数:
,使得
满足
。
根据以上条件可知,点
是多元函数
的极小点,从而
满足方程组
即
,
记
,则上述方程组可表示成
,(k=0,1,…,n)
写成矩阵形式为
,这个方程组成为法方程组,可以证明,当
线性无关时,它有唯一解。
特别地,曲线拟合的一种常用情况为代数多项式,即取
,则
(k=0,1,…,n)
故相应的法方程组变为
,这就是最小二乘法的原理。
在解决本题时,为了简便起见,我们将指数转变成代数多项式去计算。
在
两边取对数,得到
,取
,可见
是呈线性关系的。
这样我们可以方便地利用最小二乘法求取参数。
(2)流程图
整体流程图生成矩阵C流程图
4程序代码及注释
%最小二乘拟合
%a为线性拟合中的常数,b为一次项系数
%t为均方误差,maxi为最大偏差
function[a,b,t,maxi]=polyfit(x0,y0,n)
m=length(x0);p=length(y0);
%x0和y0长度不等时,报错
ifm~=p
fprintf('Error!
Pleaseinputagain!
\n');
end
%生成中间矩阵C
fori=1:
m
C(i,1)=1;
forj=2:
n+1
C(i,j)=x0(i)*C(i,j-1);
end
end
%生成系数矩阵A
A=C'*C;
%将题目中的y的每项求自然对数后得到y1
fori=1:
m
y1(i)=log(y0(i));
end
%生成法方程组的右端向量B
B=C'*y1';
%求解拟合系数
X=A\B;
%题中,a=exp(X
(1)),b=X
(2)
a=exp(X
(1));
b=X
(2);
%先求偏差平方和,再求均方误差
sum=0;
fork=1:
m
y2(k)=a*exp(b*x0(k));
l(k)=y0(k)-y2(k);
sum=sum+l(k).^2;
end
t=sqrt(sum);
%最大偏差为偏差矩阵中绝对值最大的一项
maxi=max(max(abs(l)));
end
5算例分析
1、测试示例
>>x=1:
18;
>>y=[];
>>[abtmax]=polyfit(x,y,1)
Error!
Pleaseinputagain!
a=
b=
t=
max=
2、计算过程
(1)首先输入已知点
>>x=1:
19;
>>y=[];
(2)输出结果
>>[abtmax]=polyfit(x,y,1)
a=
b=
t=
max=
其中,a,b,t,max分别为常数项,一次项系数,均方误差,最大偏差。
6讨论与结论
1、时间复杂度:
>>tic;[abtmax]=polyfit(x,y,1);toc
Elapsedtimeisseconds.
说明该算法具有一定的复杂性。
2、直观展示:
输入以下命令:
>>x=1:
19;
>>y=[];
>>fori=1:
19
y0(i)=log(y(i));
end
>>x0=0:
:
20;
>>y1=*exp*x0);
>>plot(x,y,'+')
>>holdon
>>plot(x0,y1)
>>plot(x,y0,'+')
>>y2=log+*x0;
>>holdon
>>plot(x0,y2,'-')
可得如下图形:
参考文献
[1]易大义,沈云宝,李有法.计算方法(第2版),浙江大学出版社..
[2]张琨高思超毕靖编著MATLAB2010从入门到精通电子工业出版社