最小二乘拟合实验报告

实验名称： 最小二乘拟合

1 引言

在科学实验和生产实践中，经常要从一组实验数据(xi,yi)(i1,2,合的问题。

多项式的插值虽然在一定程度上解决了由函数表求函数近似表达式的问题，但用它来解决这里的问题，是有明显的缺陷的。首先，由实验提供的数据往往有测试误差。如果要求近似曲线y=φ(x)严格地通过所给的每个数据点(xi,yi)，就会使曲线保留原来的测试误差，因此当个别数据的误差较大的时候，插值的效果是不理想的。其次，当实验数据较多时，用插值法得到的近似表达式，明显缺乏实用价值。在实验中，我们常常用最小二乘法来解决这类问题。

定义i(xi)yi为拟合函数在xi处的残差。为了是近似曲线能尽量反映所给数据点的变化趋势，我们要求|i|尽可能小。在最小二乘法中，我们选取(x)，使得偏差平方和最小，即

,m)出发，寻求函

数y=f（x）的一个近似表达式y=φ(x)，称为经验公式，从几何上来看，这就是一个曲线拟

ii1m2[(xi)yi]i1m2min，这就是最小二乘法的原理。

2 实验目的和要求

运用matlab编写.m文件，要求用最小二乘法确定参数。

以下一组数据中x与y之间存在着yaebx的关系，利用最小二乘法确定式中的参数a和b，并计算相应的军方误差与最大偏差。数据如下： x y x y 1 0.8 11 5.46 2 2.38 12 6.53 3 3.07 13 10.9 4 1.84 14 16.5 5 2.02 15 22.5 6 1.94 16 35.7 7 2.22 17 50.6 8 2.77 18 61.6 9 4.02 19 81.8 10 4.76 3 算法原理与流程图

（1） 原理

最小二乘是要求对于给定数据列(xi,yi)(i1,2,{0(x),1(x),,m)，要求存在某个函数类

n(x)}(nm)中寻求一个函数：

*ann(x)，使得*(x)满足

**(x)a00(x)a1*1(x)

[(xi)yi]i*1n2min(x)[(xi)yi]。 i21n*根据以上条件可知，点(a0,a1*,*,an)是多元函数

S(a0,a1,*的极小点，从而a0,a1*,,an)[akk(xi)yi] ik210mn*满足方程组 ,anS0(k0,1,akm,n)

mmm0(xi)a1k(xi)1(xi)即a0k(xi)i1i1mank(xi)n(xi)i1k(xi)yi， i1记(h,g)h(xi)g(xi)，则上述方程组可表示成i1a0(k,0)a1(k,1)写成矩阵形式为

(0,0)(0,1)(1,0)(1,1)(n,0)(n,1)an(k,n)(k,f)，(k=0,1,…，n)

(0,n)(1,n)(n,n)a0(0,f)a1(1,f)，这个方程组成为法方程组，可以证明，当an(n,f)0(x),1(x),n(x)线性无关时，它有唯一解。

特别地，曲线拟合的一种常用情况为代数多项式，即取

0(x)1,1(x)x,n(x)x，则(j,k)(k,f)nxi1mjixkixiji1mk

xikyi (k=0,1,…，n) i1m故相应的法方程组变为

mmxii1mxnii1xii1mxii1m2xini1m1xi1mn1xii1mxi2ni1nimmxii1a0maxiyii1，这就是最小二乘法的原理。 manxnyiii1在解决本题时，为了简便起见，我们将指数转变成代数多项式去计算。

在yaebx两边取对数，得到lnylnabx，取y(1)lny,x(1)x，可见y(1),x(1)是呈线性关系的。这样我们可以方便地利用最小二乘法求取参数。 （2）流程图

输入 x i ,y i( i  1 ,2, , m ) 及m，n i=1，2，…，m 生成中间矩阵C ci11j=2，3，…，n+1 生成法方程组的系数矩T阵 A  C C 生成法方程组的右端向TY 量 b  C cijxi*cij1解法方程组得 ai(i0,1,,n)输出 a i( i  0,1, , n )

整体流程图

生成矩阵C流程图

4 程序代码及注释

%最小二乘拟合 %a为线性拟合中的常数，b为一次项系数 %t为均方误差，maxi为最大偏差 function [a,b,t,maxi]=polyfit(x0,y0,n) m=length(x0) ;p=length(y0); %x0和y0长度不等时，报错

if m~=p fprintf('Error! Please input again!\\n'); end %生成中间矩阵C for i=1:m C(i,1)=1; for j=2:n+1 C(i,j)=x0(i)*C(i,j-1); end end %生成系数矩阵A A=C'*C; %将题目中的y的每项求自然对数后得到y1 for i=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; for k=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=[0.8 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8]; >> [a b t max]=polyfit(x,y,1) Error! Please input again!

a = 0.7185 b = 0.2227 t = 31.6255 max = 22.0631 2、计算过程 （1）首先输入已知点 >> x=1:19; >> y=[0.8 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8]; （2）输出结果 >> [a b t max]=polyfit(x,y,1) a = 0.6814 b = 0.2306 t = 38.3255

max = 27.3047 其中，a，b，t，max分别为常数项，一次项系数，均方误差，最大偏差。 6 讨论与结论

1、 时间复杂度： >> tic;[a b t max]=polyfit(x,y,1);toc Elapsed time is 0.859861 seconds. 说明该算法具有一定的复杂性。 2、 直观展示： 输入以下命令： >> x=1:19; >> y=[0.8 2.38 3.07 1.84 2.02 1.94 2.22 2.77 4.02 4.76 5.46 6.53 10.9 16.5 22.5 35.7 50.6 61.6 81.8]; >> for i=1:19 y0(i)=log(y(i)); end >> x0=0:0.01:20; >> y1=0.6814*exp(0.2306*x0); >> plot(x,y,'+') >> hold on >> plot(x0,y1) >> plot(x,y0,'+') >> y2=log(0.6814)+0.2306*x0; >> hold on >> plot(x0,y2,'-') 可得如下图形：

参考文献

[1] 易大义，沈云宝，李有法. 计算方法(第2版)，浙江大学出版社. p.29-53. [2] 张琨 高思超 毕靖 编著 MATLAB2010从入门到精通 电子工业出版社